Compare Euler' s Method and the Backward Euler Method

We will compare the solutions given by Euler’s Method and the Backward Euler Method for the following initial value problem:

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

Remember that the exact solution is y(t)=(e^(-10t)). As we saw last time, Euler’s method is unstable for Δt>0.2:

y[0]=1;
f[t_,y_]:=-10y;
Δt=0.22;
numsteps=10;
(* this Do loop computes the Euler's method approximation *)
Do[y[k+1]=f[k*Δt,y[k]]Δt+y[k],{k,0,numsteps-1}];
(* now return a table of points ((Subscript[t, k]),(Subscript[y, k])) *)
eulerSol=Table[{k*Δt,y[k]},{k,0,numsteps}]

eulerPlot=ListPlot[eulerSol,Joined->True,Mesh->All,AxesLabel->{"t","y[t]"}]

Now let's try the Backward Euler Method:

Clear[y]
y[0]=1;
Δt=0.22;
numsteps=10;
(* this Do loop computes the Backward Euler Method approximation *)
Do[y[k+1]=y[k]/(1+10Δt),{k,0,numsteps-1}];
(* now return a table of points ((Subscript[t, k]),(Subscript[y, k])) *)
backwardEulerSol=Table[{k*Δt,y[k]},{k,0,numsteps}]
TableForm[backwardEulerSol]

backwardEulerPlot=ListPlot[backwardEulerSol,Joined->True,Mesh->All,AxesLabel->{"t","y[t]"},PlotStyle->Orange]

Show[eulerPlot,backwardEulerPlot]

Improved Euler’s Method

Consider the initial-value problem:

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

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

First, set up some variables:

f[t_,y_]:=1-y;
y[0]=0;
Δt=0.1;
t[k_]:=k*Δt

Now run a loop that computes the approximation:

Do[
 ytemp=y[k]+f[ t[k],y[k] ]*Δt; (* compute the temporary value of (Subscript[y, k+1]) *)
y[k+1]=y[k]+(f[ t[k],y[k] ]+f[ t[k+1],ytemp ]) Δt/2 (* actual approximation of (Subscript[y, k+1]) *)
, {k,0,9}]

TableForm[Table[{t[k],y[k]},{k,0,10}]]

Plot the approximate solution:

Show[
ListPlot[Table[{t[k],y[k]},{k,0,100}]],
Plot[1-Exp[-t],{t,0,1},PlotStyle->Green]
]

Our approximation at t=1 is:

y[10]

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

yexact=N[1-Exp[-1]]

The error is:

yexact-y[10]

Error Comparison: Euler, Backward Euler, and Improved Euler

For ease of comparison, we will define a Module that computes the Euler’s method approximation:

eulerApprox[f_,y0_,Δt_,numsteps_]:=Module[{y,t},
y[0]=y0;
t[k_]:=k*Δt;
Do[y[k+1]=f[t[k],y[k]]Δt+y[k],{k,0,numsteps-1}];
y[numsteps] (*this is the return value *)
]

We will also define a Module that computes Improved Euler approximation:

improvedEulerApprox[f_,y0_,Δt_,numsteps_]:=Module[{y,t, ytemp},
y[0]=y0;
t[k_]:=k*Δt;
Do[ ytemp=y[k]+f[ t[k],y[k] ]*Δt;
y[k+1]=y[k]+(f[ t[k],y[k] ]+f[ t[k+1],ytemp ]) Δt/2,
{k,0,numsteps-1}];
y[numsteps] (*this is the return value *)
]

We will also define a Module that computes Backward Euler approximation. Since the Backward Euler method is implicit, we will use NSolve to compute each step of the approximation.

backwardEulerApprox[f_,y0_,Δt_,numsteps_]:=Module[{y,t},
y[0]=y0;
t[k_]:=k*Δt;
Do[
val=NSolve[s== y[k]+f[t[k+1],s]Δt,s][[1,1]];
y[k+1]=s/.val
,{k,0,numsteps-1}];
y[numsteps] (*this is the return value *)
]

Now compute the errors obtained when using these methods to solve the following problem for y(1):

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

f[t_,y_]:=1-y
exact=1-Exp[-1];
Abs[exact-eulerApprox[f,0,0.1,10]]
Abs[exact-backwardEulerApprox[f,0,0.1,10]]
Abs[exact-improvedEulerApprox[f,0,0.1,10]]

We can display the errors for different step sizes in a table.

steps={10,100,1000,10000};
TableForm[
Table[
{Abs[exact-eulerApprox[f,0,1./n,n]],
Abs[exact-backwardEulerApprox[f,0,1./n,n]],
Abs[exact-improvedEulerApprox[f,0,1./n,n]]},
{n,steps}],
TableHeadings->{steps,{"Euler","Backward","Improved"}}]

Another example:

Now compute the errors obtained when using these methods to solve the following problem for y(1):

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

f[t_,y_]:=2y-t
exact = (1+Exp[2])3/4; (* exact value of y(1)  *)
Abs[exact-eulerApprox[f,1,0.1,10]]
Abs[exact-backwardEulerApprox[f,1,0.1,10]]
Abs[exact-improvedEulerApprox[f,1,0.1,10]]

See how error decreases with step size for Improved Euler’s method:

steps={10,20,100,200,1000,2000};
TableForm[
Table[
{Abs[exact-eulerApprox[f,1,1./n,n]],
Abs[exact-backwardEulerApprox[f,1,1./n,n]],
Abs[exact-improvedEulerApprox[f,1,1./n,n]]},
{n,steps}],
TableHeadings->{steps,{"Euler","Backward","Improved"}}]

.00072/4