Question

In the study of population dynamics one of the most famous models for a growing but bounded population is the logistic equation dP/dt = P(a − bP), where a and b are positive constants. First find the general solution, then given a particular solution with p(0)=1. Then plot the solution in MATLAB with two sets of (a,b) values that you choose.

Answer 1

clear all
close all
%exact soultion
syms p(t) a b
eqn = diff(p,t) == p.*(a-b.*p);
cond = p(0) == 1;
Psol(t) = dsolve(eqn,cond);
fprintf('General Exact solution for given ode P(t)=')
disp(Psol)

%exact soultion
syms p(t)
eqn = diff(p,t) == p.*(10-5.*p);
cond = p(0) == 1;
Psol(t) = dsolve(eqn,cond);
fprintf('Exact solution for given ode for a=10; b=5 P(t)=')
disp(Psol)

%function for which Euler method have to calculate
a=10; b=5;
f=@(t,P) P.*(a-b*P);
P_in=1; t_in=0; t_end=2;
h=0.01;
[t_result,P_result] = euler_method(f,P_in,t_in,t_end,h);
figure(1)
hold on
plot(t_result,P_result)
plot(t_result,Psol(t_result))
xlabel('time')
ylabel('P(t)')
title('P(t) vs. t plot a=10 b=5')
legend('Euler Solution','Exact Solution')
box on

%exact soultion
syms p(t)
eqn = diff(p,t) == p.*(15-3.*p);
cond = p(0) == 1;
Psol(t) = dsolve(eqn,cond);
fprintf('Exact solution for given ode for a=15; b=3 P(t)=')
disp(Psol)

%function for which Euler method have to calculate
a=15; b=3;
f=@(t,P) P.*(a-b*P);
P_in=1; t_in=0; t_end=2;
h=0.01;
[t_result,P_result] = euler_method(f,P_in,t_in,t_end,h);
figure(2)
hold on
plot(t_result,P_result)
plot(t_result,Psol(t_result))
xlabel('time')
ylabel('P(t)')
title('P(t) vs. t plot a=15 b=3')
legend('Euler Solution','Exact Solution')
box on

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%Matlab code for Euler's forward
function [t_result,y1_result] = euler_method(f,y0,t0,tend,h)
%function for Euler equation solution
  
    %all step size
    N=(tend-t0)/h;
    %Initial values

    %t end values
    tn=t0:h:tend;
    % Euler steps
    y1_result(1)=y0;
    t_result(1)=t0;
    fprintf('\n\tFor step size %2.2f\n',h)
    for i=1:length(tn)-1
        t_result(i+1)= t_result(i)+h;
        y1_result(i+1)=y1_result(i)+h*double(f(t_result(i),y1_result(i)));
        %fprintf('At x=%2.2f y(%2.2f)=%f\n',t_result(i+1),t_result(i+1),y1_result(i+1))
    end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%