Math 448/548 - Example 6.2 - Simulation of Whale Problem

This Matlab code illustrates the Euler method applied to the Whale Problem problem.

clear all, close all, clc

alpha = 10^(-7);

M=10; % number of sample points
x1min=0; x1max=80000; % domain specification
x2min=0; x2max=400000;
[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 are rescaled!)'); hold on
xlabel('Blue Whales'); ylabel('Fin Whales');

Choose initial condition and terminal time T, plus the time step h

x=[5000;70000];
T = 200; % terminal time

h=1; % time step (the smaller, the better)
N = floor(T/h); % numer of iterations

This is the iteration and the resultion approximate solution in the phase portrait.

fprintf('    Blue Whales   Fin Whales\n\n')
fprintf(' Initial population\n      %5.2f        %5.2f\n', x(1), x(2))

for n=1:N
    rhs = [0.05*x(1)*(1-x(1)/150000)-alpha*x(1)*x(2);
           0.08*x(2)*(1-x(2)/400000)-alpha*x(1)*x(2)];

    xnew = x + h*rhs;
    plot([x(1),xnew(1)], [x(2),xnew(2)],'ro',...
        'MarkerFaceColor','k','MarkerSize',2), hold on
    x = xnew;
end
xlabel('Blue whales'); ylabel('Fin Whales');
fprintf('\n Final population \n     %5.2f        %5.2f\n', x(1), x(2))
    Blue Whales   Fin Whales

 Initial population
      5000.00        70000.00

 Final population 
     28220.23        386275.74

One can plot also the dynamics of each individual species as a function of time.

t=0:h:T;
x=[5000;70000]; X=x;
for n=1:N
    rhs = [0.05*x(1)*(1-x(1)/150000)-alpha*x(1)*x(2);
           0.08*x(2)*(1-x(2)/400000)-alpha*x(1)*x(2)];
    x = x + h*rhs;
    X=[X x];
end

figure
subplot(2,1,1)
plot(t,X(1,:),'bo','MarkerFaceColor','k','MarkerSize',2)
title('Blue whale dynamics vs time');
xlabel('t'); ylabel('Blue Whales');
subplot(2,1,2)
plot(t,X(2,:),'go', 'MarkerFaceColor','k','MarkerSize',2)
title('Fin whale dynamics vs time');
xlabel('t'); ylabel('Fin Whales');