This Matlab code illustrates the Eigenvalue Method for studying stability of steady states for CONTINUOUS-TIME dynamical system described in the tree problem.
Declare the variables
clear all, close all, clc syms x1 x2
Write down the RHS for the governing system of differential equations (assuming the rate of decay due to competition is 1/2 the rate due to overcrowding)
f1 = 0.1*x1*(1-(1/10000)*x1- 0.5*(1/10000)*x2); f2 = 0.25*x2*(1-(1/6000)*x2- 0.5*(1/6000)*x1);
Find the equilibrium state(s)
[x1steady,x2steady] = solve(f1,f2);
N=length(x1steady);
fprintf('The equilibrium points are\n')
disp([x1steady x2steady])
The equilibrium points are [ 0, 0] [ 0, 6000] [ 10000, 0] [ 28000/3, 4000/3]
The following code displays the direction field (f1,f2)
M=15; % number of samples points x1min=0; x1max=11000; % domain specification x2min=0; x2max=8000; ezplot(f1,[x1min x1max x2min x2max]); hold on ezplot(f2,[x1min x1max x2min x2max]); hold on [X1,X2]=meshgrid(x1min:(x1max-x1min)/M:x1max,x2min:(x2max-x2min)/M:x2max); dX1 = 0.1*X1.*(1-(1/10000)*X1-0.5*(1/10000)*X2); % x1-component dX2 = 0.25*X2.*(1-(1/6000)*X2-0.5*(1/6000)*X1); % x2-component quiver(X1,X2,dX1,dX2); % matlab routine axis([x1min x1max x2min x2max]); title('Direction field (the vectors are rescaled!)'); xlabel('Hardwoods'); ylabel('Softwoods'); hold on
First, compute the Jacobian matrix
DF = [diff(f1,x1), diff(f1,x2); diff(f2,x1), diff(f2,x2)];
and then, for each equilibrium point, compute its eigenvalues:
for i=1:N x1num = double(x1steady(i)); x2num = double(x2steady(i)); A=subs(DF,[x1,x2],[x1num,x2num]) lambda = eig(A); fprintf('The eigenvalues for the equilibrium (') fprintf('%1.0f, %1.0f',x1num, x2num); fprintf(') are ') fprintf('%1.2f %1.2f\n',lambda(1), lambda(2)); plot(x1num,x2num,'ro','MarkerSize',10, 'MarkerFaceColor','g'); end
A = 0.1000 0 0 0.2500 The eigenvalues for the equilibrium (0, 0) are 0.10 0.25 A = 0.0700 0 -0.1250 -0.2500 The eigenvalues for the equilibrium (0, 6000) are -0.25 0.07 A = -0.1000 -0.0500 0 0.0417 The eigenvalues for the equilibrium (10000, 0) are -0.10 0.04 A = -0.0933 -0.0467 -0.0278 -0.0556 The eigenvalues for the equilibrium (9333, 1333) are -0.12 -0.03
Conclusion: Based on the eigenvalues obtained above, one can decide which equilibria are stable (sink point), semi-stable (saddle point) or unstable (source point). You can visualize this by plotting a few trajectories of the dynamical system, e.g. using pplane7.m