ORDINARY DIFFERENTIAL EQUATIONS POWERTOOLLesson 14 -- Euler's MethodProf. Douglas B. MeadeIndustrial Mathematics InstituteDepartment of MathematicsUniversity of South CarolinaColumbia, SC 29208 URL: http://www.math.sc.edu/~meade/E-mail: meade@math.sc.eduCopyright \251 2001 by Douglas B. MeadeAll rights reserved-------------------------------------------------------------------

14.A Explicit Implementation of Euler's Method 14.A-1 Example 1 14.A-2 Example 214.B dsolve and Euler's Method 14.B-1 Example 1 (revisited) 14.B-2 Example 2 (revisited)14.C Example 3

restart;with( DEtools ):with( plots ):interface( rtablesize=25 );

Euler's method for the solution of a first-order IVP NiMvKiYlI2R5RyIiIiUjZHhHISIiLSUiRkc2JCUieEctJSJ5RzYjRiw=, NiMvLSUieUc2IyYlInhHNiMiIiEmRiVGKQ== can be summarized by the formulae NiMvJiUieEc2IywmJSJuRyIiIkYpRiksJiZGJTYjRihGKSUiaEdGKQ== NiMvJiUieUc2IywmJSJuRyIiIkYpRiksJiZGJTYjJSJoR0YpKiZGLUYpLSUiRkc2JCYlInhHNiNGKCZGJUY0RilGKQ== where h is the stepsize.A simple implementation of Euler's method that accepts the function F, initial time NiMmJSJ4RzYjIiIh, initial position NiMmJSJ5RzYjIiIh, stepsize NiMlImhH, and number of steps NiMlIk5H as input would beEuler0 := proc( F, x0, y0, h, N ) local i, L, X; X := evalf( < x0 | y0 > ); L := X; for i from 1 to N do X := X + < h | h*F(X[1],X[2]) >; L := L, X; end do; return < L >;end proc:As a test, compute the Euler's Method solution to NiMvKiYlI2R5RyIiIiUjZHhHISIiLCQqJiIiI0YmJSJ5R0YmRig=, NiMvLSUieUc2IyIiISIiIg==on the interval NiM3JCIiISIiIg== with NiMlImhH = 0.1.Euler0( (x,y)->-2*y, 0, 1, 0.1, 10 );A more sophisticated implementation of 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 asEuler := proc( ode, ic, domain, N ) local h, i, t, y, F, L, X; 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 X := X + < h | h*F(X[1],X[2]) >; L := L, X; end do; return < <t|y>, L >end proc:

The approximate solution to the IVPode1 := diff( y(x), x ) = -2*y(x);ic1 := y(0)=1;on the interval NiM3JCIiISIiIg== by Euler's Method with 10 subdivisions would be obtained with the commandEuler( ode1, ic1, x=0..1, 10 );

For a second example, use Euler's Method with NiMlIk5H = 2, 4, and 8 subdivisions to find an approximate value for NiMtJSJ5RzYjIiIj where the IVP isode2 := diff( y(t), t )*sin(t) + y(t) = 3;ic2 := y(1)=2;and the parameters for the numeric solution area2 := 1;b2 := 2;N2 := [ 2, 4, 8 ];The three stepsizes to be used appear in the listH2 := map( n->(b2-a2)/n, N2 );and the results of the calculations are obtained viafor k from 1 to 3 do sol2 := Euler( 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 Euler's Method for a first-order system are not significantly different or more difficult, but will not be considered at this time.

While it is not difficult to implement Euler's Method in Maple, there is no real reason to do so. The dsolve command can be used to obtain approximate solutions to IVPs for first-order ODEs, including systems.To illustrate, revisit Examples 1 and 2.

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.In Example 1, the numeric parameters werex0 := 0;h1 := 0.1;N1 := 10;A list of values of the independent variable NiMlInhH at which the corresponding values of the dependent variable NiMlInlH are to be obtained isx_list := Array( 1..N1+1, i -> x0 + (i-1)*h1 );Then, the table of approximate solution values computed using Euler's Method isdsolve( { ode1, ic1 }, y(x), type=numeric, method=classical, stepsize=h1, output=x_list );If the results are to be plotted, then the dsolve and odeplot commands can be used as follows, resulting in Figure 14.1.sol := dsolve( { ode1, ic1 }, y(x), type=numeric, method=classical, stepsize=h1 ):odeplot( sol, [x,y(x)], 0..1, title="Figure 14.1" );

The three approximations to NiMtJSJ5RzYjIiIj are obtained viafor k from 1 to 3 do sol2 := dsolve( { ode2, ic2 }, y(t), type=numeric, method=classical, stepsize=H2[k] ): Y2[k] := eval( y(t), sol2(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 >;

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 Euler's Method.Consider the problem of obtaining a solution to the IVPode3 := diff( y(t) , t ) = sin( y(t) ) ;ic3 := y(0) = 1;Euler0((t,y)->sin(y),0,1,2,4);sol3e(2); sol3e(4); sol3e(6);on the interval NiM3JCIiISIiKQ==.Euler's Method with NiMlIk5H = 4 subdivisions yields Figure 14.2, provided odeplot is instructed to plot five points since the initial point must be counted.sol3e := dsolve( { ode3, ic3 }, y(t), type=numeric, method=classical, stepsize=2 ):euler_plot := odeplot( sol3e, [t,y(t)], 0..8, style=line, color=BLUE, numpoints=5 ):display( euler_plot, title="Figure 14.2" );To determine if this approximation is reasonable, superimpose this solution on the slope field, as shown in Figure 14.3.slope_field := DEplot( ode3, y(t), t=0..8, y=0..4, arrows=SMALL ):display( [ slope_field, euler_plot ], title="Figure 14.3" );On the interval NiM3JCIiISIiIw== the Euler solution does not look too bad. However, on NiM3JCIiIyIiJQ== the Euler solution does not follow the slope field and is a much poorer approximation to the true solution. Note also that the Euler solution repeatedly crosses the equilibrium solution NiMvLSUieUc2IyUidEclI1BpRw==. To obtain a reasonable approximation on the entire interval using Euler's Method, a smaller stepsize is required. Figure 14.4 shows the result of computing with a stepsize of 0.5.sol3e2 := dsolve( { ode3, ic3 }, y(t), type=numeric, method=classical, stepsize=0.5 ):euler_plot2 := odeplot( sol3e2, [t,y(t)], 0..8, style=line, color=CYAN, thickness=3, numpoints=17 ):display( [ slope_field, euler_plot2 ], title="Figure 14.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 14.5 shows the Euler solution with stepsize 2 as the blue curve, the Euler solution with stepsize 0.5 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, euler_plot, euler_plot2 ], title="Figure 14.5" );

[Back to ODE Powertool Table of Contents]