Math 448/548 - Example 4.2 - Whale Problem

This Matlab code illustrates the solution of steady state analysis for the dynamical system describing the blue and fin whale populations.

Contents

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;

Steady State Analysis

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.

Sensitivity analysis to the parameter *alpha*

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)