This Matlab code illustrates the solution of steady state analysis for the dynamical system describing the blue and fin whale populations.
clear all, close all, clc
syms x1 x2
Choose a numerical value for the rate of competition between species (this is the most uncertain assumption):
alpha=10^(-7);
The dynamical system is given by dx1/dt=f1; dx2/dt=f2 where
f1=0.05*x1*(1-x1/150000) - alpha*x1*x2; f2=0.08*x2*(1-x2/400000) - alpha*x1*x2;
The steady states are found by solving the system f1=0, f2=0:
[x1steady,x2steady]=solve(f1,f2);
disp('The equilibrium points are')
disp([x1steady x2steady])
The equilibrium points are [ 0, 0] [ 150000, 0] [ 0, 400000] [ 600000/17, 6500000/17]
The following code displays the direction field f1,f2
M=10; % number of samples points x1min=0; x1max=900000; % domain specification x2min=0; x2max=600000; [X1,X2]=meshgrid(x1min:(x1max-x1min)/M:x1max,x2min:(x2max-x2min)/M:x2max); dX1=0.05*X1.*(1-X1/150000) - alpha*X1.*X2; % x1-component dX2=0.08*X2.*(1-X2/400000) - alpha*X1.*X2; % x2-component quiver(X1,X2,dX1,dX2); % matlab routine axis([x1min x1max x2min x2max]); title('Direction field (the vectors may be rescaled!)'); hold on xlabel('Blue Whales'); ylabel('Fin Whales');
TO PLOT THE PHASE PORTRAIT (including DIRECTION FIELD, THE STEADY STATES (EQUILIBRIUM) POPULATIONS AND A FEW SOLUTION CURVES PLEASE USE THE pplane7.m INTERACTIVE PROGRAM (available for download on the E-Companion site -> Doc Sharing, or at <http://math.rice.edu/~dfield).>
Here we can only easily plot the nullclines
ezplot(f1,[0 900000 0 600000]), hold on
ezplot(f2,[0 900000 0 600000])
Conclusion: There is a stable equlibrium in which both species will coexist (x1>0,x2>0). The stability is a consequence of the phase portrait (done separately with pplane7.m) or via the stability analysis (next chapter in the book). Any initial condition will evolve to this stable equlibrium.
The main question is, will there exist a stable coexistence equilibrium even when the value of alpha is so uncertain?
syms alpha
f1=.05*x1*(1-x1/150000)-alpha*x1*x2;
f2=.08*x2*(1-x2/400000) - alpha*x1*x2;
[x1steady,x2steady] = solve(f1,f2)
x1steady = 0 0 150000 150000*(-1+8000000*alpha)/(-1+15000000000000*alpha^2) x2steady = 0 400000 0 400000*(-1+1875000*alpha)/(-1+15000000000000*alpha^2)
A run of the previous lines shows that there are 4 equilibrium points, only the fourth indicating nonzero steady coexistence levels.
x1alpha=x1steady(4); x2alpha=x2steady(4); pretty(x1alpha), pretty(x2alpha)
-1 + 8000000 alpha 150000 -------------------------- 2 -1 + 15000000000000 alpha -1 + 1875000 alpha 400000 -------------------------- 2 -1 + 15000000000000 alpha
The values for alpha for which at least one species becomes extinct:
alpha1 = solve(x1alpha); alpha2 = solve(x2alpha); format short e double(alpha1), double(alpha2)
ans = 1.2500e-007 ans = 5.3333e-007
Using a graphical approach we study the range of values for alpha for which there exists a coexistence equilibrium (x1a>0,x2a>0)
figure ezplot(x1alpha,[0 8*10^(-7)]), hold on grid on ezplot(x2alpha,[0 8*10^(-7)]), hold on title('Level of coexisting populations vs parameter \alpha');
To be able to manipulate data and illustrate the results graphically, it is always best to work with numerical values:
a = linspace(0,9*10^(-7)); % generates a row of values between 0 and 9*10^(-7) x1a=subs(x1alpha,alpha,a); x2a=subs(x2alpha,alpha,a); ind=find(x1a>0 & x2a>0); %finds the indices for which both x1 and x2 are positive plot(a(ind),x1a(ind),'bo'); hold on plot(a(ind), x2a(ind),'ro');
From the graph above, it is clear that there will exist a coexistence equilibrium state when alpha is between 0 and 1.25*10^(-7), or greater than 5.3*10^(-7).
Here we perform a sensitivity of the equlibirum to the value \alpha=10^(-7)
format bank
Sx1alpha = diff(x1alpha, alpha)*(alpha/x1alpha);
Sx2alpha = diff(x2alpha, alpha)*(alpha/x2alpha);
Sx1a = subs(Sx1alpha, alpha, 10^(-7))
Sx2a = subs(Sx2alpha, alpha, 10^(-7))
Sx1a = -3.65 Sx2a = 0.12
This means that the blue whale population is more sensitive to the parameter alpha, but the conclusion still remains true (the equlibrium in the first quadrant will remain stable, even when alpha is not certain)