Programming Project: Solving Systems of Differential Equations Numerically

Solving Systems of Differential Equations
Numerically
~

Programming Project for Math 15

% Complete this program to perform the explicit midpoint rule on a % differential
equation system. We are given a differential equation % of the form dydt=f(t,y),
where y=y(t) outputs a 2×1 column vector for % each (scalar) input t. Here, we
assume the first component of the output % represents the number of rabbits
at time t, and the second component of % the output represents the number of
wolves at time t. Once completed, the % code should output 179 rabbits and 15
wolves at the end. Please read the % entire code before making changes.
clear all % This command erases past variables & other data format compact %
This suppresses extra lines in the output
f = @(t,y) [.6y(1)-.013y(1)y(2);.01y(1)y(2)-.82y(2)]; % f is a function of a scalar
t and a 2×1 vector y % dydt=f(t,y) is the differential equation system to solve
% this is an example of a Lotka-Volterra predator-prey model
a = 0; % lower bound for t b = 10; % upper bound for t N = 40; % panel count
h = (b-a)N; % step size y0 = [200;20]; % initial values of y
t = linspace(a,b,N+1); % creates a row vector t of N+1 equally-spaced %
elements from a to b, including endpoints. % We think of t as housing the time
values over % which we are interested in obtaining a solution.
w = zeros(2,N+1); % creates a 2x(N+1) matrix w, which we will fill in with %
the approximate solution. % w(1,i) is our approximation to the first component
of % the solution at time t(i); w(2,i) is our approximation % to the second
component of the solution at time t(i).
w(,1) = y0; % sets the first column of w to be the initial values
%————you only need to modify code below this line———— for i=1N
w(,i+1)=w(,i); % (modify this) end %————you only need to modify code
above this line————
hold on % allows us to superimpose plots plot(w(1,),w(2,)) % plots the rabbit
vs wolf population counts for all t u = gradient(w(1,)); % calculates differences
between successive iterates v = gradient(w(2,)); quiver(w(1,),w(2,),u,v); % gives
directional arrows on the graph, so that % we can see how the populations

1

evolve over % time title(‘Phase Plane of Rabbits vs. Wolves’) xlabel(‘Rabbits’);
ylabel(‘Wolves’)
disp(sprintf(‘n ~~ ~’)); disp(sprintf([‘final rabbit count %1.0f’],w(1,end)));
~
disp(sprintf([‘final wolf count %1.0f’],w(2,end)));