Math 448/548 - Example 5.1 - Tree Problem

This Matlab code illustrates the Eigenvalue Method for studying stability of steady states for CONTINUOUS-TIME dynamical system described in the tree problem.

Contents

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

Stability Analysis via Eigenvalue Method

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