Euler's Method Error Demonstration

Consider the initial-value problem:

(dy/dt)=1-y,  y(0)=0

We will investigate the error in approximating y(1) using Euler’s method.

First, compute the exact value of y(1):

By separation of variables, we find the exact solution is y(t)=1-(e^(-t)). Thus, the exact value of y(1) is:

y1=1-Exp[-1]

N[y1]

Now we will write a Module to compute the error in approximating y(1) when using Euler’s method with q steps:

EulerError[q_]:=Module[{y},
y[0]=0;
Δt=1/q;
f[t_,y_]:=1-y;

(* the following loop does the work of Euler's method *)
Do[ y[n+1]=y[n]+f[n*Δt,y[n]]*Δt, {n,0,q-1} ];

(* compute the error *)
error = Abs[y1-y[q]];
N[error] (* return the error, as a decimal *)
]

Try it for 20 steps:

EulerError[20]

Try it for 100 steps:

EulerError[100]

Make a table of the errors for various numbers of steps:

steps=Table[100s,{s,1,10}]

Table[EulerError[n],{n,steps}]

TableForm[Table[{n,EulerError[n]},{n,steps}],TableHeadings->{None, {"steps","error"}}]

Make a plot of the error:

lp=ListPlot[ Table[{n,EulerError[n]},{n,steps}] ]

Note that the error is proportional to Δt=(1/n), where n is the number of steps:

Show[{ Plot[0.1845/n,{n,1,1000},PlotStyle->Orange],lp}]

Numerical Instability in Euler’s Method

For some differential equations, Euler’s method requires a very small step size to avoid numerical instability. Consider the initial value problem:

(dy/dt)=-10y,  y(0)=1

Note that the exact solution is y(t)=(e^(-10t)).

Observe what happens when we use Euler’s method with step size Δt>0.2.

Clear[y]
f[t_,y_]:=-10y
y[0]=1
Δt=0.22;
steps=Floor[3/Δt]    (* approximate y(3) *)
Do[ y[n+1] = f[Δt*n, y[n]]*Δt + y[n], {n,0,steps}]
eulerTable = Table[{Δt*n,y[n]},{n,0,steps+1}]

TableForm[eulerTable,TableHeadings-> {None,{"t","y[t]"}}]

We will plot the approximate and exact solutions:

Show[
ListPlot[eulerTable,Joined->True,Mesh->All,AxesLabel->{"t","y[t]"}],
Plot[Exp[-10t],{t,0,3},PlotRange->All,PlotStyle->Green]
]

Try the Backward Euler’s Method

The Backward Euler’s Method uses f((Subscript[t, n+1]),(Subscript[y, n+1])) to approximate the slope of the solution at ((Subscript[t, n]),(Subscript[y, n])). For the previous problem, this produces the recursive formula (Subscript[y, n+1])=((Subscript[y, n])/(1+10Δt)). Let’s see how this works in practice:

Clear[y]

(* fill in backward Euler's method code here *)

backwardEulerTable = Table[{Δt*n, y[n]},{n,0,steps+1}];
TableForm[backwardEulerTable,TableHeadings-> {None,{"t","y[t]"}}]

Show[
ListPlot[backwardEulerTable,Joined->True,Mesh->All,AxesLabel->{"t","y[t]"},PlotRange->{0,1}],
Plot[Exp[-10t],{t,0,3},PlotRange->All,PlotStyle->Green]
]