ORDINARY DIFFERENTIAL EQUATIONS POWERTOOL

Unit 16 -- Runge-Kutta Method

Outline of Unit 16

Initialization

> restart;

> with( DEtools ):

> with( plots ):

> with( linalg ):

Warning, the name changecoords has been redefined

Warning, the name adjoint has been redefined

Warning, the protected names norm and trace have been redefined and unprotected

>

The Euler and Improved Euler Methods are members of the family of Runge-Kutta methods for solving an IVP

,

The method generally referred to as the second-order Runge-Kutta Method ( RK2 ) is defined by the formula

( )

where h is the stepsize.

A simple implementation of the second-order Runge-Kutta Method that accepts the function F, initial time , initial position , stepsize , and number of steps as input would be

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

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

> X := evalf( [ x0, y0 ] );

> L := X;

> for i from 1 to N do

> f := F(op(X));

> X1 := X + [ h, h*f ];

> Xt := X + 2/3*[ h, h*f ];

> X2 := X + [ h, h*F(op(Xt)) ];

> X := (X1+3*X2)/4;

> L := L, X;

> end do;

> return matrix( N+1, 2, [ L ] );

> end proc;

>

As a test, compute the Runge-Kutta Method solution to

,

on the interval [ 0, 1 ] with =0.1.

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

>

A more sophisticated implementation of the Runge-Kutta 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

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

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

> 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

> f := F(op(X));

> X1 := X + [ h, h*f ];

> Xt := X + 2/3*[ h, h*f ];

> X2 := X + [ h, h*F(op(X1)) ];

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

> L := L, X;

> end do;

> return matrix(N+2,2,[[t,y],L])

> end proc:

>

The approximate solution to the IVP

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

> ic1 := y(0)=1;

>

on the interval [0,1] by RK2 with 10 subdivisions would be obtained with the command

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

>

For a second example, use RK2 with = 2, 4, and 8 subdivisions to find an approximate value for where

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

> ic2 := y(1)=2;

>

> a2 := 1;

> b2 := 2;

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

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

>

> sol2 := RungeKutta2( ode2, ic2, t=a2..b2, N2[1] ):

> Y2[1] := sol2[N2[1]+2,2];

>

> sol2 := RungeKutta2( ode2, ic2, t=a2..b2, N2[2] ):

> Y2[2] := sol2[N2[2]+2,2];

>

> sol2 := RungeKutta2( ode2, ic2, t=a2..b2, N2[3] ):

> Y2[3] := sol2[N2[3]+2,2];

>

These results can be summarized in a table

> v0 := vector( [ 'h', 'N', 'Y(2)' ] ):

> v1 := vector( [ H2[1], N2[1], Y2[1] ] ):

> v2 := vector( [ H2[2], N2[2], Y2[2] ] ):

> v3 := vector( [ H2[3], N2[3], Y2[3] ] ):

> stackmatrix( v0, v1, v2, v3 );

>

Implementations of the Runge-Kutta Methods for a first-order system are not significantly different or more difficult, but will not be considered at this time.

>

To request the use of the second-order Runge-Kutta method in Maple's numerical computations, use method=classical[rk2] . The third- and fourth-Order Runge-Kutta Methods are utilized when method=classical[rk3] or method=classical[rk4] is specified. Also, recall that Maple's default numerical method is the Fehlberg fourth-fifth order Runge-Kutta method ( method=rkf45 ).

To illustrate the use of dsolve to obtain approximations using RK2, revisit the two examples considered in the previous section.

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.

> x0 := 0:

> h1 := 0.1:

> N1 := 10:

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

>

Then, the table of approximate solution values computed using RK2 is

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

> method=classical[rk2], stepsize=h1, value=x_list );

>

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

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

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

> odeplot( sol, [x,y(x)], 0..1 );

>

Construction of the table of results when the second-order Runge-Kutta method is used is accomplished with three separate calls to dsolve :

> sol2 := dsolve( { ode2, ic2 }, y(t), type=numeric,

> method=classical[rk2], stepsize=H2[1] ):

> Y2[1] := eval( y(t), sol2(2) );

>

> sol2 := dsolve( { ode2, ic2 }, y(t), type=numeric,

> method=classical[rk2], stepsize=H2[2] ):

> Y2[2] := eval( y(t), sol2(2) );

>

> sol2 := dsolve( { ode2, ic2 }, y(t), type=numeric,

> method=classical[rk2], stepsize=H2[3] ):

> Y2[3] := eval( y(t), sol2(2) );

>

The final table is

> v0 := vector( [ 'h', 'N', 'Y(2)' ] ):

> v1 := vector( [ H2[1], N2[1], Y2[1] ] ):

> v2 := vector( [ H2[2], N2[2], Y2[2] ] ):

> v3 := vector( [ H2[3], N2[3], Y2[3] ] ):

> stackmatrix( 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 the second-order Runge-Kutta 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 [0,8].

The Runge-Kutta Method with = 4 subdivisions yields

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

> method=classical[rk2], stepsize=2 ):

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

> color=BLUE, numpoints=4 ):

> display( rk2_plot );

>

To determine if this approximation is reasonable, superimpose this solution on the slope field.

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

> display( [ slope_field, rk2_plot ] );

>

On the interval [0,2] the Runge-Kutta solution does not look too bad. However, on [2,4] the Runge-Kutta solution does not follow the slope field and is a much poorer approximation to the true solution. This solution is very similar to the one obtained with the Improved Euler Method. To obtain a reasonable approximation on the entire interval using the Runge-Kutta Method, a smaller stepsize is required.

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

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

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

> color=CYAN, numpoints=10 ):

> display( [ slope_field, rk2_plot2 ] );

>

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.

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

> arrows=NONE, linecolor=GREEN ):

> display( [ slope_field, sol_plot, rk2_plot, rk2_plot2 ] );

>