The Euler Method

Often it is not possible or desirable to solve a differential equation,

$$\frac{dP}{dt} = f(P),$$

analytically, and one turns to numerical or computational methods. A numerical method seeks to approximate the solution to the equation at discrete times. Time is subdivided into intervals of length $\Delta t$, so that $t_n = n \Delta t$, and then the method approximates the solution at those times, $P_n \approx P(t_n)$. One of the oldest ideas for doing this is the Euler method. Since $\frac{dP}{dt} \approx \frac{\Delta P}{\Delta t}$ one may write

$$\frac{\Delta P}{\Delta t} = \frac{P(t_n+\Delta t)-P(t_n)}{\D... ...frac{P_{n+1} - P_n} {\Delta t} \approx \frac{dP}{dt} = f(P_n).$$

Solving this expression for $P_{n+1}$ you end up with a discrete equation which predicts a future value of $P$, $P_{n+1}$, in terms of a past value:

$$P_{n+1} = P_n + \Delta t \ f(P_n).$$

This can be used to approximate solutions to the differential equation.

For example, suppose you want to solve the exponential growth equation for aphids, but include a constant `harvesting' rate due to predation by ladybird beetles, as in

$$\frac{dP}{dt} = r P - h.$$

Using the Euler approximation, with $f(P)=rP$, gives an approximate discrete rule :

$$P_{n+1} = P_n + r \Delta t \ P_n - \Delta t \ h.$$

A MATLAB script implementing the Euler approximation would look like

 %  variables for the discretization:

       tmax=10;             % set the time to finish solving
       N=100;               % set number of time intervals
       dt=tmax/N;           % determine delta t
       t=linspace(0,tmax,N+1);     % not required for the DE, but useful for
                                   % plotting

 %  set parameters for the ODE:
       h=.2;                % harvesting rate
       r=.1;                % growth rate
       p(1)=3;              % set the initial population
 
       p=zeros(1,N+1);      % define an array to hold the solution

       for n=1:N            % begin loop  for Euler solution

          p(n+1)=(1+r*dt)*p(n)-h*dt;

       end                  % end loop  for Euler solution
       
       plot(t,p)            % plot the solution

Save this script in an m-file and run it to see how solutions look. The differential equation has a steady state at $\bar{P}=\frac{h}{R}$; can you confirm this numerically? Is the steady state stable or unstable? In this particular case we know the analytic solution; as we calculated in class

$$P=\frac{h}{r} + \left( P_0-\frac{h}{r} \right) e^{rt}.$$

Having run the script above, now try

       pactual=h/R+(p(1)-h/r)*exp(r*t);
       hold on, plot(t,pactual,'g'), hold off

which should plot the actual solution on the same axes as your last numericl solution. What do you think?

	EXERCISE 1: Apply the Euler method to solve the logistic equation $$\frac{dP}{dt} = r P \left( 1 - \frac{P}{K} \right), \quad P(0) = P_0.$$ with $r=.5$, $K=10$ and $P_0=.1$, over a time interval from 0 to 50. Try changing the number of steps ($N$); start with $N=200$ and then decrease to $N=10$. What do you observe? Do you have an explanation? Compare your discretization to the discrete logistic equation and think about stability.

As is illustrated in the previous exercise, it is possible for the Euler method (and, in fact, for any numerical approach) to go wrong, particularly when $\Delta t$ becomes large. In addition, the behavior of dynamics calculated using the Euler approximation generally `lag' actual system dynamics, as we will see when we compare Euler solutions to the analytic solution of the logistic equation (in the next EXERCISE). However, MATLAB has several built-in solvers for differential equations which avoid most numerical problems and are highly accurate; we will explore using these in the next sections.