Euler Approximations
of The Initial Value Problem:

x[t0] = x0
dx = f[t,x] dt

by
K. D.  Stroyan
University of Iowa

copyright 2007 by Keith Stroyan - All rights reserved.

Special Mathematica Functions used in this NoteBook

eulerApprox[f_,{t0_,x0_},dt_,final_] :=
Module[{steps,G},
steps = Floor[final/dt];
G[{tt_,xx_}] := {tt+dt,xx + f[{tt,xx}] dt};
NestList[ G , {t0,x0} , steps ]]

Notebook Overview

This NoteBook uses Euler's approximation  to approximate  x[t]  when  x  is given by a starting value, x[t0] = x0, and a differential equation

(dx/dt)= f[t,x]

This method is  the hand computation used in solving the S-I-R model of Chapter 2.

The Method

The way Euler's method works is as follows.  The the  increment (or microscope) approximation for the change in a smooth function is

x[t +δt] = x[t] + x'[t] δt + ϵ δt

with  ϵ≈0 small when  δt≈0 is small.Since  x'[t]  is given by the formula  f[t,x[t]],if we know a starting value  x[(Subscript[t, 0])]=(Subscript[x, 0]),we can substitute this in  f[t0,(Subscript[x, 0])] for  x'[t0],

x[(Subscript[t, 0])+δt]=x[(Subscript[t, 0])]+f[(Subscript[t, 0]),(Subscript[x, 0])] δt+ϵ δt
or
x[(Subscript[t, 0])+δt]≈x[(Subscript[t, 0])]+f[(Subscript[t, 0]),(Subscript[x, 0])] δt

This procedure can now be continued to find  x[((Subscript[t, 0])+δt) + δt], using (Subscript[x, 1])= the approximation for  x[(Subscript[t, 0]) +δt]  from the previous line and (Subscript[t, 1]) = (Subscript[t, 0]) + δt.  This gives

x[((Subscript[t, 0])+δt)+δt]=x[(Subscript[t, 0])+δt]+f[(Subscript[t, 0])+dt,x[(Subscript[t, 0])+δt]] δt+ϵ dt
so
x[t1+δt]=x[(Subscript[t, 1])]+f[(Subscript[t, 1]),(Subscript[x, 1])] dt+(Subscript[ϵ, 2]) dt
or
x[(Subscript[t, 0])+2 dt]≈x[(Subscript[t, 1])]+f[(Subscript[t, 1]),(Subscript[x, 1])] δt=(Subscript[x, 2])

Why Does This Converge?

Continuing this way, we can approximate  x[(Subscript[t, 0])], x[(Subscript[t, 0])+δt], x[(Subscript[t, 0])+2δt], x[(Subscript[t, 0])+3δt], x[(Subscript[t, 0])+4δt],...  The errors at each step are of the form η δt,

η δt = (f[(Subscript[t, j]),(Subscript[x, j])] - f[(Subscript[t, j]),x[(Subscript[t, j])]]) δt + ϵ δt = ((f[(Subscript[t, j]),(Subscript[x, j])] - f[(Subscript[t, j]),x[(Subscript[t, j])]]) + ϵ) δt

and we make  1/δt steps,the total error is roughly the maximum of the  η's.Each  η-error is small,because  ϵ  is small and  (Subscript[x, j])≈x[(Subscript[t, j])],so  f[(Subscript[t, j]),(Subscript[x, j])]≈f[(Subscript[t, j]),x[(Subscript[t, j])]],by continuity of f[t,x].

This recursive process is a function,

Input  {t,x} -> Output {t+δt, x[t+δt]}

Let's call this function G:

G[{t,x}] = {t+δt,x[t+δt] } = {t+δt, x[t] + f[t,x[t]] δt}

Recursive computation means we do the "G" operation over and over,

{(Subscript[t, 0]),x[(Subscript[t, 0])]} = {(Subscript[t, 0]),(Subscript[x, 0])}
{(Subscript[t, 0])+δt,x[(Subscript[t, 0])+δt]} = G[{(Subscript[t, 0]),x[(Subscript[t, 0])]}]
{(Subscript[t, 0])+δt+δt, x[(Subscript[t, 0])+δt+δt]} = G[{(Subscript[t, 0])+δt],x[(Subscript[t, 0])+δt]}] = G[ G[{(Subscript[t, 0]),x[(Subscript[t, 0])]}] ]
{(Subscript[t, 0])+3 δt],x[(Subscript[t, 0])+3 δt]} = G[{(Subscript[t, 0])+2 δt,x[(Subscript[t, 0])+2 δt]}] =  G[G[ G[{(Subscript[t, 0]),x[(Subscript[t, 0])]}] ]]
	etc.

Summary

Euler's approximation is a clear recursive process, but converges slowly.  Its biggest advantage is that it only uses "where you start" and "how you change" (in the form of the increment approximation) to find the solution.  Each recursive step changes  t  by  δt and  x by the amount dictated from the differential equation.

Symbolic solutions or Mathematica's NDSolve[.] can be used to get accurate approximate solutions, but these methods don't use the description of change as clearly.

References to the Text

Chapter 4

In Chapter 4 we study the idea of using a differential equation to determine the value of temperature as a function of time.  The EulerApprox program uses the increment (or microscope) approximation to recursively compute a change such as the decrease in temperature (i.e. -dC/dt) for the cool canary problem.  Refer to Exercise 4.2.1 below to recall what this initial value problem looks like in Mathematica.

Chapter 8

In Chapter 8 we begin the study of the natural exponential function

y[t]  = (e^t)

by "defining" it to be the solution of the initial value problem:

y[0] = 1
	dy/dt = y

We study that definition by approximating the solution to the initial value problem.   The discrete system that moves in small steps of size δt is Euler's method and in this specific case results in the approximation

(e^t)≈(((1+(t/(δt))))^(1/δt))

Chapter 21

In Chapter 21, we return to Euler's method as part of the basic theory of differential equations.

Step-by-step Computation

Enter the following computation to see what happens.  In this case, the process function  G  adds  dt  to  t  and  f[x,t] dt  to  x.  We start at  t = t0, so get times t0, t0+dt, t0 + 02 dt, ... and we start at  x = x0, getting  positions  given in the Output.

Clear[t0,x0,G,f];
G[{t_,x_}]:={t+dt,x+f[t,x] dt};
dt=0.5;
finalTime=1.5;
solnList=NestList[G,{t0,x0},Floor[(finalTime/dt)]]

Now let's use a specific differential equation,  dx/dt = Sin[x*t], and specific initial conditions:

Clear[f]
f[t_,x_]:=Sin[x t]
t0=1.`;
x0=1.`;
dt=0.5`;
finalTime=2;
solnList=NestList[G,{t0,x0},Floor[(finalTime/dt)]]
firstPlot=ListPlot[solnList,Joined->True,PlotStyle->{Red}]

Now let's solve with a smaller step size (leaving all else the same):

dt=0.25`;
solnList=NestList[G,{t0,x0},Floor[(finalTime/dt)]];
secondPlot=ListPlot[solnList,Joined->True,PlotStyle->{DarkGreen}];Show[{firstPlot,secondPlot}]

Smaller yet:

dt=0.05`;
solnList=NestList[G,{t0,x0},Floor[(finalTime/dt)]];
thirdPlot=ListPlot[solnList,Joined->True,PlotStyle->{Blue}];Show[firstPlot,secondPlot,thirdPlot]

Complete Sample Computations

The differential equation for "logistic" growth at per capita rate  r  and carrying capacity  c  is

(dx/dt) = r x (1 - x/c)

Here are some solutions at different initial values and step sizes:

Clear[f,G,x,t,c,r];
f[t_,x_]:=r x (1-(x/c));
c=2;
r=1;
G[{t_,x_}]:={t+dt,x+f[t,x] dt};
t0=0;
x0=1;
finalTime=10;
dt=1.`;
steps=Floor[(finalTime/dt)];
Print["Time step = ",dt]
solnList=NestList[G,{t0,x0},steps];
ListPlot[solnList,Joined->True]

A smaller time step:

dt=0.1`;
steps=Floor[(finalTime/dt)];
Print["Time step = ",dt]
solnList=NestList[G,{t0,x0},steps];
ListPlot[solnList,Joined->True]

A tiny time step:

dt=0.01`;
steps=Floor[(finalTime/dt)];
Print["Time step = ",dt]
solnList=NestList[G,{t0,x0},steps];
ListPlot[solnList,Joined->True]

Animation of Convergence

First we can "package" Euler's method as a Mathematica Module of the Initialization.  The form of input is
	eulerApprox[ Rate function, {initial pair}, time step, final time]

Using this Module takes the form:

Clear[f,t,x];
f[{t_,x_}]:=x (1-(x/2));
soln1=eulerApprox[f,{0,1},0.5,6.0]
soln2=eulerApprox[f,{0,3},0.5,6.0]

plot1=ListPlot[soln1,PlotStyle->{Red}]
plot2=ListPlot[soln2,PlotStyle->{Blue}]
Show[plot1,plot2,PlotRange->{0,3}]

Now we animate solutions with these two initial conditions as the time step gets shorter.

Clear[f,t,x];
f[{t_,x_}]:=x (1-(x/2));
Manipulate[
dt=1./(2^n);
soln1=eulerApprox[f,{0,1},dt,6.`];plot1=ListPlot[soln1,PlotStyle->{Red}];soln2=eulerApprox[f,{0,3},dt,6.`];plot2=ListPlot[soln2,PlotStyle->{Blue}];Show[plot1,plot2,Ticks->None,PlotRange->{-1,4},AxesLabel->{"t","x"}]
,{n,1,5}]

Exercises

Exercise 4.2.1 The Dead Canary by Euler's Method

<Enter> the cell below to see graphically what how the temperature of the dead canary changes with respect to time for

Clear[f,G,t,c];
f[c_,t_]:=-k c;
k=0.03`;
G[{t_,c_}]:={t+dt,c+f[c,t] dt};
t0=0;
c0=75;
finalTime=10.`;
dt=0.01`;
steps=Floor[(finalTime/dt)];
solnList=NestList[G,{t0,c0},steps];
ListPlot[solnList,Joined->True,PlotRange->{0,75},AxesLabel->{"t","c"}]
Print["The final time-temp pair =",Take[solnList,-1]]

Exercise 2.3.4 & 4.2.2: Geometry of Euler's Method

To sketch the solution of dy/dx = -y  starting from  y = 1 when x = 0, we start at (x,y) = (0,1) and move at slope - y = -1, say we move distance 0.25,

step1={{0,1},{0.25`,0.75`}}
ListPlot[step1,Ticks->None,AxesLabel->{"x","y"},PlotRange->{{0,2},{0,1}},Joined->True]

Nest, we move from {0.25,0.75} at slope -y = -0.75.  This means as  x  goes from  0.25 to 0.50,  y  goes from  0.75 to  0.75 - 0.75*0.25:

step2={{0,1},{0.25`,0.75`},{0.5`,0.75`-0.75` 0.25`}}
ListPlot[step2,Ticks->None,AxesLabel->{"x","y"},PlotRange->{{0,2},{0,1}},Joined->True]

We need a formula so we can make the machine do these steps.  We always change  x by adding dx, in our case, starting with x = 0 and  dx = 0.25, so x+dx = 0.25, then 0.50, ...

Input = x,	->	Output = x + dx

We change  y  so that  (Change in y)/(Change in x) = (Change in y)/dx = -y.  But the change in  y  is the new value minus the old value, so

(Output y  - Input y)/dx = -y
		
	Output y = Input y + (-y) dx

We express this input -> output function by

G[{x,y}] = { x+dx , y - y dx }
In Mathematica notation:

Clear[G,x,y,dx];
G[{x_,y_}]:={x+dx,y-y dx};
dx=0.25`;
soln1too=NestList[G,{0,1},1]
dots=ListPlot[soln1too,Ticks->None,AxesLabel->{"x","y"},PlotRange->{{0,2},{0,1}},PlotStyle->{PointSize[0.02`],Red},DisplayFunction->Identity];,lines=ListPlot[soln1too,Ticks->None,AxesLabel->{"x","y"},PlotRange->{{0,2},{0,1}},PlotStyle->{Red},Joined->True,DisplayFunction->Identity];,Show[dots,lines,DisplayFunction->$DisplayFunction]

soln2too=NestList[G,{0,1},2]
dots=ListPlot[soln2too,Ticks->None,AxesLabel->{"x","y"},PlotRange->{{0,2},{0,1}},PlotStyle->{PointSize[0.02`],Red},DisplayFunction->Identity];,lines=ListPlot[soln2too,Ticks->None,AxesLabel->{"x","y"},PlotRange->{{0,2},{0,1}},PlotStyle->{Red},Joined->True,DisplayFunction->Identity];,Show[dots,lines,DisplayFunction->$DisplayFunction]

soln3=NestList[G,{0,1},3]
dots=ListPlot[soln3,Ticks->None,AxesLabel->{"x","y"},PlotRange->{{0,2},{0,1}},PlotStyle->{PointSize[0.02`],Red},DisplayFunction->Identity];,lines=ListPlot[soln3,Ticks->None,AxesLabel->{"x","y"},PlotRange->{{0,2},{0,1}},PlotStyle->{Red},Joined->True,DisplayFunction->Identity];,Show[dots,lines,DisplayFunction->$DisplayFunction]

ListPlot[NestList[G,{0,1},8],Ticks->None,AxesLabel->{"x","y"},PlotRange->{{0,2},{0,1}},Joined->True,PlotStyle->{Blue}]

Exercise 21.3.1 & Problem 21.4

Copy and modify the computation above to find solutions to the logistic equation:

dx = x (1 - x/2) dt
for
	t0 = 0	& x0 = 0.1
		& x0 = 0.5
		& x0 = 1.0
		& x0 = 1.5
		& x0 = 2.5
		& x0 = 3.0

Euler's Method in 2-D: The S-I-R model revisited

Euler approximation works when applying the concept to the S-I-R model.  <Enter> the cell below to see how the fraction of susceptible and infectious people change with respect to time with the given initial conditions.

Print["Spread of rubella in fractions of a population."];
Print["Contact number, c = ",c,", Infectious period 1/b =",(1/b)];
Print["Daily contact rate, a =",a];
c=6.8`;
b=(1/(11.`));
a=b c;
sRate[s_,i_]:=-a s i;
iRate[s_,i_]:=(a s-b) i;
G[{t_,s_,i_}]:={t+dt,s+sRate[s,i] dt,i+iRate[s,i] dt};
t0=0.`;
s0=(1/(3.`));
i0=(1/(30.`));
initialValues={t0,s0,i0};
dt=N[(1/2)];
finalTime=N[90];
steps=Floor[(finalTime/dt)];
Print["Increments of dt = ",dt," for ",steps," recursive steps."];
solnList=NestList[G,initialValues,steps];
suscs=Flatten[Take[Transpose[solnList],{2}],1];
times=Flatten[Take[Transpose[solnList],{1}],1];
infes=Flatten[Take[Transpose[solnList],-1],1];
sPairs=Transpose[{times,suscs}];
iPairs=Transpose[{times,infes}];
eulerPlot=ListPlot[sPairs,Joined->True,AxesLabel->{"t","s"}]
ListPlot[iPairs,Joined->True,AxesLabel->{"t","i"}]
Print["The last susceptible value:"];
Spair=Flatten[Take[sPairs,-1]];
eulerSdata=Spair[[2]];
Print[eulerSdata];