mth3060104m9.html

1.4 Numerical Technique: Euler's Method

September 9, 1998

Example One

Consider the initial value problem , , for from to . We'll start with a step size .

> h:=0.5:          # This is our step-size.
n:=ceil(5.0/h): # This calculates how many steps
                # of size h we need to take, to get to 5.0.
t:=0.0:
w:=0.0:         # The initial condition.
data[0]:=[t,w]: # We store the initial point so we can
                # plot it later.
for i from 1 to n do    # This is how to make a block of
                       # code get executed n times.
  mk := (3-w)*(w+1):
  w := w + h*mk:
  t := t + h:
  data[i] := [t,w]: # Store the current point so we can
                     # plot it later.
od:              # This ("do" backwards) marks the end
                # of the "do" block.
t:='t': w:='w': h:='h': # Clean up

datalist1 := [ seq( data[i], i=0..n ) ];
                # This puts all the data into a list
                # suitable for plotting.
plot( datalist1, 0..5, 0..4, style=POINT, color=red);
                # This plots the (approximation to the)
                # solution, with a horizontal (t) range
                # from 0 to 5, and a vertical (w) range
                # of 0 to 4. The POINT style means each
                # data point is marked with a symbol of
                # some kind. If "style=POINT" were omitted,
                # the data points would be joined with
                # straight line segments (try it!).

From what we can figure out qualitatively, this seems odd. Let's look at the direction field.

> with(DEtools):
plot1:=plot( datalist1, 0..5, 0..4,style=point,color=red):
plot2:=dfieldplot(diff(w(t),t)=(3-w)*(w+1), [w(t)],t=0..5,w=0..4,color=blue,arrows=line):
with(plots):
display([plot1,plot2]);

We do better if we decrease the step size. Let's try .

> h:=0.1:          # This is our new step-size.
n:=ceil(5.0/h):
t:=0.0:
w:=0.0:
data[0]:=[t,w]:
for i from 1 to n do
  mk := (3-w)*(w+1):
  w := w + h*mk:
  t := t + h:
  data[i] := [t,w]:
od:
t:='t': w:='w': h:='h':
datalist1 := [ seq( data[i], i=0..n ) ];
plot1:=plot( datalist1, 0..5, 0..4, color=red):
display([plot1,plot2]);

Of course, we can actually solve this initial value problem.

> simplify(dsolve({diff(w(t),t)=(3-w(t))*(w(t)+1),w(0)=0},w(t)));

And plot the solution and the numerical approximation on the same plot.

> plot3:=plot(3*(exp(4*t)-1)/(3+exp(4*t)),t=0..5,color=green):
display([plot3,plot2,plot1]);

Example Two

Let us consider the initial value problem , . Suppose we want to estimate . We can start with a step size of .

> h:=0.5:
n:=ceil(3.0/h):
t:=0.0:
y:=1.0:
data[0]:=[t,y]:
for i from 1 to n do
  mk := sin(y):
  y := y + h*mk:
  t := t + h:
  data[i] := [t,y]:
od:
t:='t': y:='y': h:='h':
datalist1 := [ seq( data[i], i=0..n ) ];
plot1:=plot( datalist1, 0..3, 0..4,style=point, color=red):
display([plot1]);

Now let's try it with a step size half as big.

> h:=0.25:
n:=ceil(3.0/h):
t:=0.0:
y:=1.0:
data[0]:=[t,y]:
for i from 1 to n do
  mk := sin(y):
  y := y + h*mk:
  t := t + h:
  data[i] := [t,y]:
od:
t:='t': y:='y': h:='h':
datalist2 := [ seq( data[i], i=0..n ) ];
plot2:=plot( datalist2, 0..3, 0..4,style=point, color=blue):
display([plot2,plot1]);

> y31:=op(2,op(nops(datalist1),datalist1));
y32:=op(2,op(nops(datalist2),datalist2));
deltay3:=y32-y31;
y3:=y32+deltay3;

Let's see how this compares with the exact solution.

> dsolve({diff(y(t),t)=sin(y(t)),y(0)=1},y(t));

Okay, so now we substitute in into the two "solutions" presented here.

> evalf(subs(t=3.0,arctan(2*exp(t+ln(-1/2*(2+2*sqrt(1+tan(1)^2))/tan(1)))/(1+exp(2*t+2*ln(-1/2*(2+2*sqrt(1+tan(1)^2))/tan(1)))),(-exp(2*t+2*ln(-1/2*(2+2*sqrt(1+tan(1)^2))/tan(1)))+1)/(1+exp(2*t+2*ln(-1/2*(2+2*sqrt(1+tan(1)^2))/tan(1)))))));

> y3exact:=evalf(subs(t=3,arctan(2*exp(t+ln(-1/2*(2-2*sqrt(1+tan(1)^2))/tan(1)))/(1+exp(2*t+2*ln(-1/2*(2-2*sqrt(1+tan(1)^2))/tan(1)))),(-exp(2*t+2*ln(-1/2*(2-2*sqrt(1+tan(1)^2))/tan(1)))+1)/(1+exp(2*t+2*ln(-1/2*(2-2*sqrt(1+tan(1)^2))/tan(1)))))));

> y3exact-y3;

Notice that actual overshoots the exact value.

Example Three

Here's an example of how the numerical method can break down.
Consider , .

> plot2:=dfieldplot(diff(y(t),t)=sin(y)*exp(t), [y(t)],t=0..5,y=0..6,color=blue,arrows=line):
display(plot2);

Notice that there is an equilibrium solution at . Let's try Euler's method with a step size of 0.1.

> h:=0.1:
n:=ceil(6.0/h):
t:=0.0:
y:=5.0:
data[0]:=[t,y]:
for i from 1 to n do
  mk := sin(y)*exp(t):
  y := y + h*mk:
  t := t + h:
  data[i] := [t,y]:
od:
t:='t': y:='y': h:='h':
datalist1 := [ seq( data[i], i=0..n ) ];
plot1:=plot( datalist1, 0..6, 0..6,color=red):
display([plot1]);

Note the misbehavior starting near . See what happens when the step size decreases.

Example Four

Here's another example of how Euler's method can "fail".
Consider , . Suppose we try to find the solution on the interval from to .

> h:=0.1:
n:=ceil(1.0/h):
t:=0.0:
y:=1.0:
data[0]:=[t,y]:
for i from 1 to n do
  mk := y^2:
  y := y + h*mk:
  t := t + h:
  data[i] := [t,y]:
od:
t:='t': y:='y': h:='h':
datalist1 := [ seq( data[i], i=0..n ) ];
plot1:=plot( datalist1, 0..1, 0..10,color=red):
display([plot1]);

Now let's try it again with a smaller step size: .

> h:=0.05:
n:=ceil(1.0/h):
t:=0.0:
y:=1.0:
data[0]:=[t,y]:
for i from 1 to n do
  mk := y^2:
  y := y + h*mk:
  t := t + h:
  data[i] := [t,y]:
od:
t:='t': y:='y': h:='h':
datalist2 := [ seq( data[i], i=0..n ) ];
plot2:=plot( datalist2, 0..1, 0..10,color=blue):
display([plot1,plot2]);

The first gave as the value of , the second . According to our rule of thumb, the true solution ought to be somewhere between and

> 9.552668001+(9.552668001-6.128898403);

So let's halve the step size once again.

> h:=0.025:
n:=ceil(1.0/h):
t:=0.0:
y:=1.0:
data[0]:=[t,y]:
for i from 1 to n do
  mk := y^2:
  y := y + h*mk:
  t := t + h:
  data[i] := [t,y]:
od:
t:='t': y:='y': h:='h':
datalist3 := [ seq( data[i], i=0..n ) ];
plot3:=plot( datalist3, 0..1, 0..10,color=green):
display([plot1,plot2,plot3]);

What? ? That's quite a bit out of the range we expect.
What is going on? Well, let's look at the analytic solution.

> dsolve({diff(y(t),t)=y(t)^2,y(0)=1},y(t));

The solution blows up at . We can't expect to get a good numerical approximation to this solution near using Euler's method.

> plot4:=plot(-1/(t-1),t=0..1,color=black):
display([plot1,plot2,plot3,plot4]);