% ### ODErkEX1.m ### 09.19.14 % numerically solve the Logistic equation % P'(t)= k*P(1-P/L) % Program calculates in four ways: 1&2. solve explicitly using Euler and RK4, % 3. solve via ode45 (via external function call) and 4. actual solution clear; figure(1); clf % ************************ % User Inputs k= 1; % intrinsic growth rate (const.) L= 5; % carrying capacity (const.) Pinit= L/15; % initial condition at tI(1) tI= [0 10]; % time boundaries stepsize= 0.05; % for RK4 % ************************ % Runge-Kutta 4 (and Euler's method) m=1; % counter for j=tI(1):stepsize:tI(2) if j == tI(1) P(m)= Pinit; Peuler(m)= Pinit; else P0= k*P(m-1)*(1-P(m-1)/L); % deriv. at y=y0 (last meas.) P1= k*(P(m-1)+(stepsize/2)*P0)*(1-(P(m-1)+(stepsize/2)*P0)/L); P2= k*(P(m-1)+(stepsize/2)*P1)*(1-(P(m-1)+(stepsize/2)*P1)/L); P3= k*(P(m-1)+(stepsize)*P2)*(1-(P(m-1)+(stepsize)*P2)/L); P(m)= P(m-1) + (stepsize/6)*(P0+ 2*P1+ 2*P2+ P3); % RK4 solution Peuler(m)= P(m-1) + stepsize* P0; % also store away Euler's method value end t(m)= j; % keep track of t for plotting m= m+ 1; % increment counter end % visualize plot(t,Peuler,'g+'); grid on; hold on plot(t,P,'o--'); xlabel('t'); ylabel('P(t)'); % ************************ % can also solve via Matlab's builtin ode45, but need to use an external % function to define the ODE [tM,PM] = ode45(@(t,P) logistic(t,P,k,L),tI, Pinit); plot(tM,PM,'kx','LineWidth',2); % ************************ % also plot analytic solution A= (L-Pinit)/Pinit; sol= L./(1+A*exp(-k*t)); plot(t,sol, 'r-') legend('Euler','RK4','ode45 solution','Exact solution','Location','SouthEast') title('Solution to logisitc equation using various numerical methods')