How to use Euler's method to solve the logistic grown model?

Hi all,

I need help solving the logistic growth model (an ODE) using Euler's Method in MATLAB.

The function is: (dy/dx) = r*y*(1-(y/K)) where r is the growth rate and K is the carrying capacity.

I have solved this out by hand but I am having a difficult time implementing it as a function.

I'm meant to write a function with two inputs (a vector time, t, and an initial y, y0) and this function is meant to output a vector of solutions to the ode for each time t and plot the results.

I know there is mean to be some h as a time step to evaluate a new tangent line at the function and the smaller that step is, the more accurate the answer is, but I'm having a hard time writing the code.

Can anyone help?

[James Tursa]

What code have you written so far? Do you know what Euler's Method is? Do you know how to write a loop? Do you know how to plot?

[Star Strider]

‘Yes’ to the last two: How to solve for a dependent system of difference (recursive) equations?

[John D'Errico]

So make an effort! You should know how Euler's method works. Surely it has been discussed.

[Erin W]

@ James This is what I have and I know it's not right but it's what I've got:

% logistic growth model: (dy/dt) = ry(1-(y/K))
% r = growth rate, K = carrying cap
% function DE = odeEuler(t,y0)
DE = @(t,y0) (r.*y.*(1-(y/K)));
% this function will take in a vector t for time 
% and initial condition y0
% It will output a vector of approx sol'n to logistic eq'n at t
% inital values
r = 1;
K = 10;
t = linspace(0,10,11);
h = 0:t:0.5 %time-step determined by t where t evenly spaced over [0,10]
y = zeros(1,length(t));
for i = 2:length(t)
  y(i) = y(i-1) + h.*DE((i-1),y(i-1));
  DE = [t(i) y(i)];
end

I know how to write loops but I'm fairly new (I just started a couple of weeks ago) and I know how to plot. I also know what Euler's Method is, I'm just finding it challenging to code.

[Erin W]

@ Star Strider Nice to see you again. Last time I had a better grasp on the code. My teacher provided this support code but it seems more confusing than other similar codes I've seen.

function yE=eulerExample(h)
% Euler example uses Euler's Method and the Improved Euler's Method to
% numerically solve dy/dt=-y; y(0)=2
% Inputs: time step h (take h=0.1);
% Outputs:  numerical solution using Euler's method, yE
t=0:h:0.5;
y=[2, 1.80967, 1.6374, 1.48163, 1.34064, 1.21306];
yE=[2, 1.8, 1.62, 1.458, 1.3122, 1.118090];
yIE=[2, 1.81, 1.63805, 1.482435, 1.341603, 1.2141515];
DelE=abs(y-yE);
DelIE=abs(y-yIE);
% Plot solutions
figure
plot(t,y,t,yE,t,yIE)
xlabel('time t')
ylabel('Solutions')
legend('True Soln', 'Euler Approx', 'Improved Euler Approx')
% Numerical solutions are discrete, so they are properly plotted as points
figure
plot(t,y,'Marker','o')
hold on
plot(t,yE,'Marker','x', 'Color','green')
plot(t,yIE,'Marker','*', 'Color','red')
xlabel('time t')
ylabel('Solutions')
legend('True Soln', 'Euler Approx', 'Improved Euler Approx')
% Plot error
figure
plot(t,DelE,'Marker', 'x', 'Color','green')
hold on
plot(t,DelIE,'Marker', '*', 'Color','red')
xlabel('time t')
ylabel('Error')
legend('Euler Error', 'Improved Euler Error')
hold off

I'm totally ok with the plotting parts but I'm unsure as to how she found her y yE and yIE.

[Star Strider]

How were you told to implement the Euler method?

[Erin W]

"We will write a function in MATLAB to solve this equation numerically using Euler’s method.

Open a new file in MATLAB and declare a function odeEuler. There are 2 inputs to this function: a vector t of evenly spaced points over the time interval of integration, [0, 10] and an initial value y0. Note that your time step will be determined by the spacing in t. There is one output: the function should return a vector of the approximate solution to the logistic equation at the time points t."

[Star Strider]

I should have been clearer. Sorry.

What algorithm for the Euler method were you told to use, or was discussed in class?

Answer 1

[James Tursa]

One step of Euler's Method is simply this:

(value at new time) = (value at old time) + (derivative at old time) * time_step

So to put this in a loop, the outline of your program would be as follows assuming y is a scalar:

t = your time vector
y0 = your initial y value
dt = the time step (you write code here to calculate this from the t vector based on your instructions)
y = zeros(size(t));  % allocate a vector to hold the output values, same size as t
y(1) = y0; % the initial value
for k=2:numel(t)  % loop to find all the subsequent values y(2) through y(end)
    y(k) = y(k-1) + (you write code for derivative dy/dt here) * dt;  % The Euler step
end

So, complete the code as indicated above, and then add code to wrap all of this inside a function per the instructions and you will be done. Your function will return the y vector. (Note that t and y0 will be coming in as function arguments, not set directly in the function itself)

[Erin W]

I'm not sure what I've messed up but I'm only getting the value 1 returned.

(value at new time) = (value at old time) + (derivative at old
% time)*time_step
function DE = odeEuler(t,y0)
% initial conditions 
r = 1;
K = 10;
dt = 0:t:10; % time step
y = zeros(size(t)); % zero vector to fill in for y
y(1) = y0; % initial y
for k=2:numel(t)
  y(k) = y(k-1) + ((K*y0)/(y0+(K-y0)*e^(-r*t)))*dt;
end
DE = y;

I'm going to try to debug it because your response was incredibly helpful but if you see something really wrong, I'd love to know... :(

But thank you so much. And, when it works, when I go to plot would it be something like: plot(t,y) right?

[Erin W]

I just changed the for k=2:numel(t) index k to i because of the capital K for carrying cap.

I now have:

for i=2:numel(t)
    y(i) = y(i-1) + ((K*y0)/(y0+(K-y0)*e^(-r*t)))*dt;
end

But it's still yielding 1. :/

[Erin W]

Sorry, one more thing, I know it is supposed to look something like this, which I calculated by hand, only it's supposed to be a lot more refined.

c = [2.31969,4.50853,6.90568,8.58486,9.42826,9.78178,9.91860,9.969810,9.98891,9.99592]
t = 1:1:length(c)
plot(t,c)

And it should be approaching K at 10, which is good, but my loop must be wrong because I'm not getting a vector returned.

[James Tursa]

Problem with this line:

dt = 0:t:10; % time step

dt is the distance between the current time and the next time. Since your assignment states that these distances are constant ("... evenly space points ..."), you can simply calculate one single scalar dt value and use that at each iteration. So ask yourself, if you are given a vector of times t that is evenly spaced, how would you calculate the distance between any two of these values that are next to each other? Put that calculation on the rhs of this line. (Hint: this is not some complex formula ... it is really really simple!)

And a problem with this line:

y(i) = y(i-1) + ((K*y0)/(y0+(K-y0)*e^(-r*t)))*dt;

Your derivative code uses the y0 value. This is not what your derivative equation looks like in your original post. If you look above, the derivative equation has y on the rhs, not y0. In other words, you should be using the current value of y, not the initial value of y. The current value of y for the rhs is simply y(i-1). So replace y0 with y(i-1) on the rhs.

Another problem, and it may be a bit confusing to you because of nomenclature, is the t that you are using on the rhs. Again, the derivative equation written above implicitly assumes that this t is the current time value for the rhs, not the entire t vector coming into the function. So the confusing part is that t is used differently ... in the function argument t is a vector of times, but in the derivative equation t is a single time value (from the vector t). So, bottom line is the time value for the derivative equation for the rhs is simply going to be t(i-1) since that is the current time of the rhs stuff. So replace the t on the rhs of this line with t(i-1).

And yes, just plot(t,y) will work as long as t and y are the names of the variables holding the data in the current workspace.

[Erin W]

I got it now! Thank you so much! :)