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');