Question

Now assume that the harvesting is not done at a constant rate, but rather at rates that vary at different times of the year. This can be modeled by ??/?? = .25? (1 − ?/4 ) − ?(1 + sin(?)). This equation cannot be solved by any technique we have learned. In fact, it cannot be solved analytically, but it can still be analyzed graphically.

8. Let c=0.16. Use MATLAB to graph a slopefield and approximate solutions for several different values of p(0), and interpret what you see. Turn in the graphs together with your analysis. Note: 0<p<5, and 0<t<50 is a reasonable viewing window to start with, though you may want to change it as you proceed.

Answer 1

%%Matlab code for ode solution and direction field
clear all
close all
%function to be solved
f = @(t,p) 0.25.*p.*(1-(p/4))-0.16.*(1+sin(t));
t_ini=0;
t_end=50;
grid on
%Plotting Direction field
figure(1)
dirrfield(f,0:2:50,0:.25:5)
hold on
%plot for different initial velocity
for v0=0:1:5
[ts,ys] = ode45(f,[t_ini,t_end],v0); plot(ts,ys,'b','Linewidth',2)
end

hold off
title('phase potrait plot')
xlabel('t')
ylabel('p(t) ')

%%Matlab function for direction filed
function dirrfield(f,tval,yval)
% dirfield(f, t1:dt:t2, y1:dy:y2)
%
%   plot direction field for first order ODE y' = f(t,y)
%   using t-values from t1 to t2 with spacing of dt
%   using y-values from y1 to t2 with spacing of dy
%
%   f is an @ function, or an inline function,
%     or the name of an m-file with quotes.
%
% Example: y' = -y^2 + t
%   Show direction field for t in [-1,3], y in [-2,2], use
%   spacing of .2 for both t and y:
%
%   f = @(t,y) -y^2+t
%   dirfield(f, -1:.2:3, -2:.2:2)

[tm,ym]=meshgrid(tval,yval);
dt = tval(2) - tval(1);
dy = yval(2) - yval(1);
fv = vectorize(f);
if isa(f,'function_handle')
fv = eval(fv);
end
yp=feval(fv,tm,ym);
s = 1./max(1/dt,abs(yp)./dy)*0.75;
h = ishold;
quiver(tval,yval,s,s.*yp,0,'r'); hold on;
%quiver(tval,yval,-s,-s.*yp,0,'r');
if h
hold on
else
hold off
end
axis([tval(1)-dt/2,tval(end)+dt/2,yval(1)-dy/2,yval(end)+dy/2])
end