Example

The following procedure implements the Fourth Order Runge-Kutta method for a system of two ODE's:

>

> RK4sys2:=proc(f, g, x0, alpha, beta, h, N)

> local n,ku1,ku2,ku3,ku4,kw1,kw2,kw3,kw4;

> x(0):=x0:u(0):=alpha:w(0):=beta:

> for n from 0 to N-1 do

> ku1:=h*f(x(n),u(n),w(n)):

> kw1:=h*g(x(n),u(n),w(n)):

> ku2:=h*f(x(n)+h/2,u(n)+ku1/2,w(n)+kw1/2):

> kw2:=h*g(x(n)+h/2,u(n)+ku1/2,w(n)+kw1/2):

> ku3:=h*f(x(n)+h/2,u(n)+ku2/2,w(n)+kw2/2):

> kw3:=h*g(x(n)+h/2,u(n)+ku2/2,w(n)+kw2/2):

> ku4:=h*f(x(n)+h,u(n)+ku3,w(n)+kw3):

> kw4:=h*g(x(n)+h,u(n)+ku3,w(n)+kw3):

> u(n+1) := u(n) + (1/6)*(ku1+2*ku2+2*ku3+ku4):

> w(n+1) := w(n) + (1/6)*(kw1+2*kw2+2*kw3+kw4):

> x(n+1) := x(n) + h:

> od:

> print(x(N),u(N),w(N));

> end:

>

Then consider the initial value problem

, in , ,

with the initial conditions

, .

Using

,

this can be written as

,

,

with

,

.

>

Using the initial values , , we find for the approximation to the solution:

>

> f:=(t,y,nu)->nu;

> g:=(t,y,nu)->-2*nu/t+2*y/t^2+sin(ln(t))/t^2;

> RK4sys2(f,g,1.,1.,0.,0.2,5);

>

The exact value of is

> yexact:=1.464721061;

So the absolute error is

> abserr:=abs(yexact-u(5));

>

This can be improved by using a smaller step size:

>

> RK4sys2(f,g,1.,1.,0.,0.01,100);abserr:=abs(yexact-u(100));

>