• Outline of Lesson 15
• Initialization
• 15.A Explicit Implementation of Improved Euler's Method
  • 15.A-1 Example 1
  • 15.A-2 Example 2
• 15.B and Improved Euler's Method
  • 15.B-1 Example 1 (revisited)
  • 15.B-2 Example 2 (revisited)
• 15.C Example 3

ORDINARY DIFFERENTIAL EQUATIONS POWERTOOL

Lesson 15 -- Improved Euler's Method

Prof. Douglas B. Meade

Industrial Mathematics Institute

Department of Mathematics

University of South Carolina

Columbia, SC 29208

URL:   http://www.math.sc.edu/~meade/

E-mail: meade@math.sc.edu

Copyright  2001  by Douglas B. MeadeAll rights reserved

-------------------------------------------------------------------

Outline of Lesson 15

15.A Explicit Implementation of Improved Euler's Method

                15.A-1 Example 1

                15.A-2 Example 2

15.B dsolve and Improved Euler's Method

                15.B-1 Example 1 (revisited)

               15.B-2 Example 2 (revisited)

15.C Example 3

>

Initialization

>

restart;

>	with( DEtools ):

>	with( plots ):

Warning, the name changecoords has been redefined

>	interface( rtablesize=25 );

>

15.A Explicit Implementation of Improved Euler's Method

There is some inconsistency in the use of the names Improved Euler's Method and Modified Euler's Method. Typically, the Improved Euler's Method is the method also known as the Trapezoid Method or Heun's Method.

The Improved Euler's method for the solution of a first-order IVP

,     

can be summarized by the formulae

                                                       

                                          

  ( )  

where h is the stepsize.

By comparison, the Modified Euler's Method is typically defined to be

                                 

                  

   

A simple implementation of the Improved Euler's method that accepts the function F, initial time , initial position , stepsize , and number of steps as input would be

>	ImprEuler0 := proc( F, x0, y0, h, N )

>	local i, L, X, X1, X2;

>	X := evalf( < x0 | y0 > );

>

L := X;

>	for i from 1 to N do

>	X1 := X + < h | h*F(X[1],X[2]) >;

>	X2 := X + < h | h*F(X1[1],X1[2]) >;

>	X := (X1+X2)/2;

>	L := L, X;

>

end do;

>	return < L >;

>

end proc:

>

As a test, compute the Improved Euler's Method solution to the IVP

,   

on the interval with = 0.1.

>	ImprEuler0( (x,y)->-2*y, 0, 1, 0.1, 10 );

>

A more sophisticated implementation of the Improved Euler's method would accept as input the ODE, the initial condition, the interval on which the solution should be computed, and the number of steps. In this case, the implementation could appear as

>	ImprEuler := proc( ode, ic, domain, N )

>	local h, i, t, y, F, L, X, X1, X2;

>	t := lhs(domain);

>	y := op(0,lhs(ic));

>	h := ( op(2,rhs(domain))-op(1,rhs(domain)) )/N;

>	F := unapply( subs( y(t)=_y, solve( ode, diff(y(t),t) ) ), (t,_y) );

>	X := evalf( < op(lhs(ic)) | rhs(ic) > );

>

L := X;

>	for i from 1 to N do

>	X1 := X + < h | h*F(X[1],X[2]) >;

>	X2 := X + < h | h*F(X1[1],X1[2]) >;

>	X := (X1 + X2)/2;

>	L := L, X;

>

end do;

>	return < <t|y>, L >

>

end proc:

>

15.A-1 Example 1

On the interval , the approximate solution to the IVP

>	ode1 := diff( y(x), x ) = -2*y(x);

>	ic1 := y(0)=1;

>

by the Improved Euler's Method with 10 subdivisions would be obtained with the command

>	ImprEuler( ode1, ic1, x=0..1, 10 );

>

15.A-2 Example 2

For a second example, use the Improved Euler's Method with = 2, 4, and 8 subdivisions to find an approximate value for where is the solution of the IVP

>	ode2 := diff( y(t), t )*sin(t) + y(t) = 3;

>	ic2 := y(1)=2;

>

The parameters for the numeric solution are

>

a2 := 1;

>

b2 := 2;

>	N2 := [ 2, 4, 8 ];

>

and the list of the three decreasing stepsizes is

>	H2 := map( n->(b2-a2)/n, N2 );

>

Then, the three approximate values for are

>	for k from 1 to 3 do sol2 := ImprEuler( ode2, ic2, t=a2..b2, N2[k] ): Y2[k] := sol2[N2[k]+2,2]; print(` Y2`[k] = Y2[k]); od:

>

Alternatively, these results are summarized in the following table.

>	v0 := < 'h'   | 'N'   | 'Y(2)' >:

>	v1 := < H2[1] | N2[1] | Y2[1]  >:

>	v2 := < H2[2] | N2[2] | Y2[2]  >:

>	v3 := < H2[3] | N2[3] | Y2[3]  >:

>	< v0, v1, v2, v3 >;

>

Implementations of the Improved Euler's Method for a first-order system are not significantly different or more difficult, but will not be considered at this time.

>

15.B dsolve and Improved Euler's Method

To request the use of the Improved Euler's Method in Maple's numerical computations, use method=classical[heunform] . The Modified Euler Method, or Improved Polygon Method, is implemented with method=classical[impoly] .

To illustrate, revisit the two examples considered in the previous section.

15.B-1 Example 1 (revisited)

When an explicit table of values is needed, it is necessary to provide a list of values of the independent variable at which the approximate solution should be reported.  Such a list would be

>

x0 := 0:

>	h1 := 0.1:

>

N1 := 10:

>	x_list := Array( 1..N1+1, i -> x0 + (i-1)*h1 );

>

Then, the table of approximate solution values computed using the Improved Euler's Method is

>	dsolve( { ode1, ic1 }, y(x), type=numeric,

>	method=classical[heunform], stepsize=h1,

>	output=x_list );

>

If the results are to be plotted, then the dsolve and odeplot commands can be used as follows to create Figure 15.1.

>	sol := dsolve( { ode1, ic1 }, y(x), type=numeric,

>	method=classical[heunform], stepsize=h1 ):

>	odeplot( sol, [x,y(x)], 0..1, title="Figure 15.1");

>

15.B-2 Example 2 (revisited)

The table of results can be constructed using three calls to dsolve as follows.  The for-loop reduces the amount of typing needed with three separate executions of the statements.

>	for k from 1 to 3 do  sol2 := dsolve( { ode2, ic2 }, y(t), type=numeric,                  method=classical[heunform], stepsize=H2[k] ):  Y2[k] := eval( y(t), sol2(2) );  print(` Y2`[k] = Y2[k]); od:

>

The resulting table is

>	v0 := < 'h'   | 'N'   | 'Y(2)' >:

>	v1 := < H2[1] | N2[1] | Y2[1]  >:

>	v2 := < H2[2] | N2[2] | Y2[2]  >:

>	v3 := < H2[3] | N2[3] | Y2[3]  >:

>	< v0, v1, v2, v3 >;

>

15.C Example 3

The final example illustrates the use of the full range of Maple tools to obtain, visualize, and analyze an approximate solution to an IVP obtained by the Improved Euler's Method.

Consider the problem of obtaining a solution to the IVP

>	ode3 := diff( y(t) , t ) = sin( y(t) );

>	ic3 := y(0) = 1;

>

on the interval .

The Improved Euler's Method with = 4 subdivisions yields Figure 15.2.

>	sol3ie := dsolve( { ode3, ic3 }, y(t), type=numeric, method=classical[heunform], stepsize=2 ):

>	impreuler_plot := odeplot( sol3ie, [t,y(t)], 0..8, style=line,

>	color=BLUE, numpoints=5 ):

>	display( impreuler_plot, title="Figure 15.2" );

>

To determine if this approximation is reasonable, superimpose this solution on the slope field as done in Figure 15.3.

>	slope_field := DEplot( ode3, y(t), t=0..8, y=0..4, arrows=SMALL ):

>	display( [ slope_field, impreuler_plot ], title="Figure 15.3" );

>

On the interval the Improved Euler solution does not look too bad. However, on the Improved Euler solution does not follow the slope field and is a much poorer approximation to the true solution. Unlike the Euler Method solution in Lesson 14, which crossed the equilibrium solution, the Improved Euler's Method solution does not approach the equilibrium solution fast enough. To obtain a reasonable approximation on the entire interval using the Improved Euler's Method, a smaller stepsize is required.  The result of computing with stepsize = 0.8 appears in Figure 15.4

>	sol3ie2 := dsolve( { ode3, ic3 }, y(t), type=numeric,

>	method=classical[heunform], stepsize=0.8 ):

>	impreuler_plot2 := odeplot( sol3ie2, [t,y(t)], 0..8, style=line,

>	color=CYAN, thickness=3, numpoints=10 ):

>	display( [ slope_field, impreuler_plot2 ], title="Figure 15.4" );

>

To complete the evaluation of this approximate solution, the DEplot command is used to include the (Maple-generated approximate) solution curve for this initial condition.  Figure 15.5 shows the Improved Euler solution with stepsize 2 as the blue curve, the Improved Euler solution with stepsize 0.8 as the cyan curve, and the default (adaptive) numeric solution as the brown curve.

>	sol_plot := DEplot( ode3, y(t), t=0..8, [ [ic3] ],

>	arrows=NONE, linecolor=BROWN ):

>	display( [ slope_field, sol_plot, impreuler_plot, impreuler_plot2 ] );

>

[Back to ODE Powertool Table of Contents]