Runge-Kutta method

From Rosetta Code
Task
Runge-Kutta method
You are encouraged to solve this task according to the task description, using any language you may know.

Given the example Differential equation:

With initial condition:

and

This equation has an exact solution:

Task

Demonstrate the commonly used explicit fourth-order Runge–Kutta method to solve the above differential equation.

  • Solve the given differential equation over the range with a step value of (101 total points, the first being given)
  • Print the calculated values of at whole numbered 's () along with error as compared to the exact solution.
Method summary

Starting with a given and calculate:

then:

Ada

<lang Ada>with Ada.Text_IO; use Ada.Text_IO; with Ada.Numerics.Generic_Elementary_Functions; procedure RungeKutta is

  type Floaty is digits 15;
  type Floaty_Array is array (Natural range <>) of Floaty;
  package FIO is new Ada.Text_IO.Float_IO(Floaty); use FIO;
  type Derivative is access function(t, y : Floaty) return Floaty;
  package Math is new Ada.Numerics.Generic_Elementary_Functions (Floaty);
  function calc_err (t, calc : Floaty) return Floaty;
  
  procedure Runge (yp_func : Derivative; t, y : in out Floaty_Array;
                   dt : Floaty) is
     dy1, dy2, dy3, dy4 : Floaty;
  begin
     for n in t'First .. t'Last-1 loop
        dy1 := dt * yp_func(t(n), y(n));
        dy2 := dt * yp_func(t(n) + dt / 2.0, y(n) + dy1 / 2.0);
        dy3 := dt * yp_func(t(n) + dt / 2.0, y(n) + dy2 / 2.0);
        dy4 := dt * yp_func(t(n) + dt, y(n) + dy3);
        t(n+1) := t(n) + dt;
        y(n+1) := y(n) + (dy1 + 2.0 * (dy2 + dy3) + dy4) / 6.0;
     end loop;
  end Runge;
  
  procedure Print (t, y : Floaty_Array; modnum : Positive) is begin
     for i in t'Range loop
        if i mod modnum = 0 then
           Put("y(");   Put (t(i), Exp=>0, Fore=>0, Aft=>1);
           Put(") = "); Put (y(i), Exp=>0, Fore=>0, Aft=>8);
           Put(" Error:"); Put (calc_err(t(i),y(i)), Aft=>5);
           New_Line;
        end if;
     end loop;
  end Print;
  function yprime (t, y : Floaty) return Floaty is begin
     return t * Math.Sqrt (y);
  end yprime;
  function calc_err (t, calc : Floaty) return Floaty is
     actual : constant Floaty := (t**2 + 4.0)**2 / 16.0;
  begin return abs(actual-calc);
  end calc_err;   
  
  dt : constant Floaty := 0.10;
  N : constant Positive := 100;
  t_arr, y_arr : Floaty_Array(0 .. N);

begin

  t_arr(0) := 0.0;
  y_arr(0) := 1.0;
  Runge (yprime'Access, t_arr, y_arr, dt);
  Print (t_arr, y_arr, 10);

end RungeKutta;</lang>

Output:
y(0.0) = 1.00000000 Error: 0.00000E+00
y(1.0) = 1.56249985 Error: 1.45722E-07
y(2.0) = 3.99999908 Error: 9.19479E-07
y(3.0) = 10.56249709 Error: 2.90956E-06
y(4.0) = 24.99999377 Error: 6.23491E-06
y(5.0) = 52.56248918 Error: 1.08197E-05
y(6.0) = 99.99998341 Error: 1.65946E-05
y(7.0) = 175.56247648 Error: 2.35177E-05
y(8.0) = 288.99996843 Error: 3.15652E-05
y(9.0) = 451.56245928 Error: 4.07232E-05
y(10.0) = 675.99994902 Error: 5.09833E-05

AWK

<lang AWK>

  1. syntax: GAWK -f RUNGE-KUTTA_METHOD.AWK
  2. converted from BBC BASIC

BEGIN {

   print(" t    y         error")
   y = 1
   for (i=0; i<=100; i++) {
     t = i / 10
     if (t == int(t)) {
       actual = ((t^2+4)^2) / 16
       printf("%2d %12.7f %g\n",t,y,actual-y)
     }
     k1 = t * sqrt(y)
     k2 = (t + 0.05) * sqrt(y + 0.05 * k1)
     k3 = (t + 0.05) * sqrt(y + 0.05 * k2)
     k4 = (t + 0.10) * sqrt(y + 0.10 * k3)
     y += 0.1 * (k1 + 2 * (k2 + k3) + k4) / 6
   }
   exit(0)

} </lang>

output:

 t    y         error
 0    1.0000000 0
 1    1.5624999 1.45722e-007
 2    3.9999991 9.19479e-007
 3   10.5624971 2.90956e-006
 4   24.9999938 6.23491e-006
 5   52.5624892 1.08197e-005
 6   99.9999834 1.65946e-005
 7  175.5624765 2.35177e-005
 8  288.9999684 3.15652e-005
 9  451.5624593 4.07232e-005
10  675.9999490 5.09833e-005

BASIC

BBC BASIC

<lang bbcbasic> y = 1.0

     FOR i% = 0 TO 100
       t = i% / 10
 
       IF t = INT(t) THEN
         actual = ((t^2 + 4)^2) / 16
         PRINT "y("; t ") = "; y TAB(20) "Error = ";  actual - y
       ENDIF
 
       k1 =  t * SQR(y)
       k2 = (t + 0.05) * SQR(y + 0.05 * k1)
       k3 = (t + 0.05) * SQR(y + 0.05 * k2)
       k4 = (t + 0.10) * SQR(y + 0.10 * k3)
       y += 0.1 * (k1 + 2 * (k2 + k3) + k4) / 6
     NEXT i%</lang>

Output:

y(0) = 1            Error = 0
y(1) = 1.56249985   Error = 1.45721892E-7
y(2) = 3.99999908   Error = 9.19479201E-7
y(3) = 10.5624971   Error = 2.90956245E-6
y(4) = 24.9999938   Error = 6.23490936E-6
y(5) = 52.5624892   Error = 1.08196974E-5
y(6) = 99.9999834   Error = 1.65945964E-5
y(7) = 175.562476   Error = 2.35177287E-5
y(8) = 288.999968   Error = 3.15652015E-5
y(9) = 451.562459   Error = 4.07231605E-5
y(10) = 675.999949  Error = 5.09832905E-5

C

<lang c>#include <stdio.h>

  1. include <stdlib.h>
  2. include <math.h>

double rk4(double(*f)(double, double), double dx, double x, double y) { double k1 = dx * f(x, y), k2 = dx * f(x + dx / 2, y + k1 / 2), k3 = dx * f(x + dx / 2, y + k2 / 2), k4 = dx * f(x + dx, y + k3); return y + (k1 + 2 * k2 + 2 * k3 + k4) / 6; }

double rate(double x, double y) { return x * sqrt(y); }

int main(void) { double *y, x, y2; double x0 = 0, x1 = 10, dx = .1; int i, n = 1 + (x1 - x0)/dx; y = malloc(sizeof(double) * n);

for (y[0] = 1, i = 1; i < n; i++) y[i] = rk4(rate, dx, x0 + dx * (i - 1), y[i-1]);

printf("x\ty\trel. err.\n------------\n"); for (i = 0; i < n; i += 10) { x = x0 + dx * i; y2 = pow(x * x / 4 + 1, 2); printf("%g\t%g\t%g\n", x, y[i], y[i]/y2 - 1); }

return 0;

}</lang>output (errors are relative)
x       y       rel. err.
------------
0       1       0
1       1.5625  -9.3262e-08
2       4       -2.2987e-07
3       10.5625 -2.75462e-07
4       25      -2.49396e-07
5       52.5625 -2.05844e-07
6       100     -1.65946e-07
7       175.562 -1.33956e-07
8       289     -1.09222e-07
9       451.562 -9.01828e-08
10      676     -7.54191e-08

D

Translation of: Ada

<lang d>import std.stdio, std.math, std.typecons;

alias FP = real; alias FPs = Typedef!(FP[101]);

void runge(in FP function(in FP, in FP)

          pure nothrow @safe @nogc yp_func,
          ref FPs t, ref FPs y, in FP dt) pure nothrow @safe @nogc {
   foreach (immutable n; 0 .. t.length - 1) {
       immutable FP
           dy1 = dt * yp_func(t[n], y[n]),
           dy2 = dt * yp_func(t[n] + dt / 2.0, y[n] + dy1 / 2.0),
           dy3 = dt * yp_func(t[n] + dt / 2.0, y[n] + dy2 / 2.0),
           dy4 = dt * yp_func(t[n] + dt, y[n] + dy3);
       t[n + 1] = t[n] + dt;
       y[n + 1] = y[n] + (dy1 + 2.0 * (dy2 + dy3) + dy4) / 6.0;
   }

}

FP calc_err(in FP t, in FP calc) pure nothrow @safe @nogc {

   immutable FP actual = (t ^^ 2 + 4.0) ^^ 2 / 16.0;
   return abs(actual - calc);

}

void main() {

   enum FP dt = 0.10;
   FPs t_arr, y_arr;
   t_arr[0] = 0.0;
   y_arr[0] = 1.0;
   runge((t, y) => t * y.sqrt, t_arr, y_arr, dt);
   foreach (immutable i; 0 .. t_arr.length)
       if (i % 10 == 0)
           writefln("y(%.1f) = %.8f Error: %.6g",
                    t_arr[i], y_arr[i],
                    calc_err(t_arr[i], y_arr[i]));

}</lang>

Output:
y(0.0) = 1.00000000 Error: 0
y(1.0) = 1.56249985 Error: 1.45722e-07
y(2.0) = 3.99999908 Error: 9.19479e-07
y(3.0) = 10.56249709 Error: 2.90956e-06
y(4.0) = 24.99999377 Error: 6.23491e-06
y(5.0) = 52.56248918 Error: 1.08197e-05
y(6.0) = 99.99998341 Error: 1.65946e-05
y(7.0) = 175.56247648 Error: 2.35177e-05
y(8.0) = 288.99996843 Error: 3.15652e-05
y(9.0) = 451.56245928 Error: 4.07232e-05
y(10.0) = 675.99994902 Error: 5.09833e-05

Dart

<lang dart>import 'dart:math' as Math;

num RungeKutta4(Function f, num t, num y, num dt){

 num k1 = dt * f(t,y);
 num k2 = dt * f(t+0.5*dt, y + 0.5*k1);
 num k3 = dt * f(t+0.5*dt, y + 0.5*k2);
 num k4 = dt * f(t + dt, y + k3);
 return y + (1/6) * (k1 + 2*k2 + 2*k3 + k4);

}

void main(){

 num t  = 0;
 num dt = 0.1;
 num tf = 10;
 num totalPoints = ((tf-t)/dt).floor()+1;
 num y  = 1;
 Function f  = (num t, num y) => t * Math.sqrt(y);
 Function actual = (num t) => (1/16) * (t*t+4)*(t*t+4);
 for (num i = 0; i <= totalPoints; i++){
   num relativeError = (actual(t) - y)/actual(t);
   if (i%10 == 0){
     print('y(${t.round().toStringAsPrecision(3)}) = ${y.toStringAsPrecision(11)}  Error = ${relativeError.toStringAsPrecision(11)}');
   }
   y  = RungeKutta4(f, t, y, dt);
   t += dt;
 }

}</lang>

Output:
y(0.00) = 1.0000000000  Error = 0.0000000000
y(1.00) = 1.5624998543  Error = 9.3262010950e-8
y(2.00) = 3.9999990805  Error = 2.2986980086e-7
y(3.00) = 10.562497090  Error = 2.7546153479e-7
y(4.00) = 24.999993765  Error = 2.4939637555e-7
y(5.00) = 52.562489180  Error = 2.0584442034e-7
y(6.00) = 99.999983405  Error = 1.6594596090e-7
y(7.00) = 175.56247648  Error = 1.3395644308e-7
y(8.00) = 288.99996843  Error = 1.0922214534e-7
y(9.00) = 451.56245928  Error = 9.0182772312e-8
y(10.0) = 675.99994902  Error = 7.5419063100e-8


Fortran

<lang fortran>program rungekutta implicit none real(kind=kind(1.0D0)) :: t,dt,tstart,tstop real(kind=kind(1.0D0)) :: y,k1,k2,k3,k4 tstart =0.0D0 ; tstop =10.0D0 ; dt = 0.1D0 y = 1.0D0 t = tstart write(6,'(A,f4.1,A,f12.8,A,es13.6)') 'y(',t,') = ',y,' Error = '&

          &,abs(y-(t**2+4.0d0)**2/16.0d0)

do; if ( t .ge. tstop ) exit

  k1 = f (t           , y                 )
  k2 = f (t+0.5D0 * dt, y +0.5D0 * dt * k1)
  k3 = f (t+0.5D0 * dt, y +0.5D0 * dt * k2)
  k4 = f (t+        dt, y +        dt * k3)
  y = y + dt *( k1 + 2.0D0 *( k2 + k3 ) + k4 )/6.0D0
  t = t + dt
  if(abs(real(nint(t))-t) .le. 1.0D-12) then
     write(6,'(A,f4.1,A,f12.8,A,es13.6)') 'y(',t,') = ',y,' Error = '&
          &,abs(y-(t**2+4.0d0)**2/16.0d0)
  end if

end do contains

 function f (t,y)
   implicit none
   real(kind=kind(1.0D0)),intent(in) :: y,t
   real(kind=kind(1.0D0)) :: f
   f = t*sqrt(y)
 end function f

end program rungekutta </lang>

Output:
y( 0.0) =   1.00000000 Error =  0.000000E+00
y( 1.0) =   1.56249985 Error =  1.457219E-07
y( 2.0) =   3.99999908 Error =  9.194792E-07
y( 3.0) =  10.56249709 Error =  2.909562E-06
y( 4.0) =  24.99999377 Error =  6.234909E-06
y( 5.0) =  52.56248918 Error =  1.081970E-05
y( 6.0) =  99.99998341 Error =  1.659460E-05
y( 7.0) = 175.56247648 Error =  2.351773E-05
y( 8.0) = 288.99996843 Error =  3.156520E-05
y( 9.0) = 451.56245928 Error =  4.072316E-05
y(10.0) = 675.99994902 Error =  5.098329E-05

F#

<lang F Sharp> open System

let y'(t,y) = t * sqrt(y)

let RungeKutta4 t0 y0 t_max dt =

   let dy1(t,y) = dt * y'(t,y)
   let dy2(t,y) = dt * y'(t+dt/2.0, y+dy1(t,y)/2.0)
   let dy3(t,y) = dt * y'(t+dt/2.0, y+dy2(t,y)/2.0)
   let dy4(t,y) = dt * y'(t+dt, y+dy3(t,y))
   (t0,y0) |> Seq.unfold (fun (t,y) ->
       if ( t <= t_max) then Some((t,y), (Math.Round(t+dt, 6), y + ( dy1(t,y) + 2.0*dy2(t,y) + 2.0*dy3(t,y) + dy4(t,y))/6.0)) 
       else None
       )

let y_exact t = (pown (pown t 2 + 4.0) 2)/16.0

RungeKutta4 0.0 1.0 10.0 0.1

   |> Seq.filter (fun (t,y) -> t % 1.0 = 0.0 )
   |> Seq.iter (fun (t,y) -> Console.WriteLine("y({0})={1}\t(relative error:{2})", t, y, (y / y_exact(t))-1.0) )</lang>
Output:
y(0)=1			(relative error:0)
y(1)=1.56249985427811	(relative error:-9.32620110027926E-08)
y(2)=3.9999990805208	(relative error:-2.29869800194571E-07)
y(3)=10.5624970904376	(relative error:-2.75461533583155E-07)
y(4)=24.9999937650906	(relative error:-2.49396374552013E-07)
y(5)=52.5624891803026	(relative error:-2.05844421730106E-07)
y(6)=99.9999834054036	(relative error:-1.65945964192282E-07)
y(7)=175.562476482271	(relative error:-1.33956447156969E-07)
y(8)=288.999968434799	(relative error:-1.09222150213029E-07)
y(9)=451.56245927684	(relative error:-9.01827772459285E-08)
y(10)=675.99994901671	(relative error:-7.54190684348899E-08)

Go

Works with: Go1

<lang go>package main

import (

   "fmt"
   "math"

)

type ypFunc func(t, y float64) float64 type ypStepFunc func(t, y, dt float64) float64

// newRKStep takes a function representing a differential equation // and returns a function that performs a single step of the forth-order // Runge-Kutta method. func newRK4Step(yp ypFunc) ypStepFunc {

   return func(t, y, dt float64) float64 {
       dy1 := dt * yp(t, y)
       dy2 := dt * yp(t+dt/2, y+dy1/2)
       dy3 := dt * yp(t+dt/2, y+dy2/2)
       dy4 := dt * yp(t+dt, y+dy3)
       return y + (dy1+2*(dy2+dy3)+dy4)/6
   }

}

// example differential equation func yprime(t, y float64) float64 {

   return t * math.Sqrt(y)

}

// exact solution of example func actual(t float64) float64 {

   t = t*t + 4
   return t * t / 16

}

func main() {

   t0, tFinal := 0, 10 // task specifies times as integers,
   dtPrint := 1        // and to print at whole numbers.
   y0 := 1.            // initial y.
   dtStep := .1        // step value.
   t, y := float64(t0), y0
   ypStep := newRK4Step(yprime)
   for t1 := t0 + dtPrint; t1 <= tFinal; t1 += dtPrint {
       printErr(t, y) // print intermediate result
       for steps := int(float64(dtPrint)/dtStep + .5); steps > 1; steps-- {
           y = ypStep(t, y, dtStep)
           t += dtStep
       }
       y = ypStep(t, y, float64(t1)-t) // adjust step to integer time
       t = float64(t1)
   }
   printErr(t, y) // print final result

}

func printErr(t, y float64) {

   fmt.Printf("y(%.1f) = %f Error: %e\n", t, y, math.Abs(actual(t)-y))

}</lang>

Output:
y(0.0) = 1.000000 Error: 0.000000e+00
y(1.0) = 1.562500 Error: 1.457219e-07
y(2.0) = 3.999999 Error: 9.194792e-07
y(3.0) = 10.562497 Error: 2.909562e-06
y(4.0) = 24.999994 Error: 6.234909e-06
y(5.0) = 52.562489 Error: 1.081970e-05
y(6.0) = 99.999983 Error: 1.659460e-05
y(7.0) = 175.562476 Error: 2.351773e-05
y(8.0) = 288.999968 Error: 3.156520e-05
y(9.0) = 451.562459 Error: 4.072316e-05
y(10.0) = 675.999949 Error: 5.098329e-05

Haskell

Using GHC 7.4.1.

<lang haskell>import Data.List

dv :: Floating a => a -> a -> a dv = (. sqrt). (*)

fy t = 1/16 * (4+t^2)^2

rk4 :: (Enum a, Fractional a)=> (a -> a -> a) -> a -> a -> a -> [(a,a)] rk4 fd y0 a h = zip ts $ scanl (flip fc) y0 ts where

 ts = [a,h ..]
 fc t y = sum. (y:). zipWith (*) [1/6,1/3,1/3,1/6]
   $ scanl (\k f -> h * fd (t+f*h) (y+f*k)) (h * fd t y) [1/2,1/2,1]

task = mapM_ print

 $ map (\(x,y)-> (truncate x,y,fy x - y)) 
 $ filter (\(x,_) -> 0== mod (truncate $ 10*x) 10) 
 $ take 101 $ rk4 dv 1.0 0 0.1</lang>

Example executed in GHCi: <lang haskell>*Main> task (0,1.0,0.0) (1,1.5624998542781088,1.4572189122041834e-7) (2,3.9999990805208006,9.194792029987298e-7) (3,10.562497090437557,2.909562461184123e-6) (4,24.999993765090654,6.234909399438493e-6) (5,52.56248918030265,1.0819697635611192e-5) (6,99.99998340540378,1.6594596999652822e-5) (7,175.56247648227165,2.3517730085131916e-5) (8,288.99996843479926,3.1565204153594095e-5) (9,451.562459276841,4.0723166534917254e-5) (10,675.9999490167125,5.098330132113915e-5)</lang>

J

Solution: <lang j>NB.*rk4 a Solve function using Runge-Kutta method NB. y is: y(ta) , ta , tb , tstep NB. u is: function to solve NB. eg: fyp rk4 1 0 10 0.1 rk4=: adverb define

'Y0 a b h'=. 4{. y
T=. a + i.@>:&.(%&h) b - a
Y=. Yt=. Y0
for_t. }: T do.
  ty=. t,Yt
  k1=. h * u ty
  k2=. h * u ty + -: h,k1 
  k3=. h * u ty + -: h,k2 
  k4=. h * u ty + h,k3
  Y=. Y, Yt=. Yt + (%6) * 1 2 2 1 +/@:* k1, k2, k3, k4  
end.

T ,. Y )</lang> Example: <lang j> fy=: (%16) * [: *: 4 + *: NB. f(t,y)

  fyp=: (* %:)/                         NB. f'(t,y)
  report_whole=: (10 * i. >:10)&{       NB. report at whole-numbered t values
  report_err=: (, {: - [: fy {.)"1      NB. report errors
  report_err report_whole fyp rk4 1 0 10 0.1
0       1           0
1  1.5625 _1.45722e_7
2       4 _9.19479e_7
3 10.5625 _2.90956e_6
4      25 _6.23491e_6
5 52.5625 _1.08197e_5
6     100 _1.65946e_5
7 175.562 _2.35177e_5
8     289 _3.15652e_5
9 451.562 _4.07232e_5

10 676 _5.09833e_5</lang>

Alternative solution:

The following solution replaces the for loop as well as the calculation of the increments (ks) with an accumulating suffix. <lang j>rk4=: adverb define

'Y0 a b h'=. 4{. y
T=. a + i.@>:&.(%&h) b-a
(,. [: h&(u nextY)@,/\. Y0 ,~ }.)&.|. T

)

NB. nextY a Calculate Yn+1 of a function using Runge-Kutta method NB. y is: 2-item numeric list of time t and y(t) NB. u is: function to use NB. x is: step size NB. eg: 0.001 fyp nextY 0 1 nextY=: adverb define

tableau=. 1 0.5 0.5, x * u y
ks=. (x * [: u y + (* x&,))/\. tableau
({:y) + 6 %~ +/ 1 2 2 1 * ks

)</lang>

Use:

report_err report_whole fyp rk4 1 0 10 0.1

JavaScript

<lang JavaScript> function rk4(y, x, dx, f) {

   var k1 = dx * f(x, y),
       k2 = dx * f(x + dx / 2.0,     +y + k1 / 2.0),
       k3 = dx * f(x + dx / 2.0,     +y + k2 / 2.0),
       k4 = dx * f(x + dx,         +y + k3);
   return y + (k1 + 2.0 * k2 + 2.0 * k3 + k4) / 6.0;

}

function f(t, y) {

   return t * Math.sqrt(y);

}

function actual(t) {

   return (1/16) * (t*t+4)*(t*t+4);

}

var y = 1.0,

   x = 0.0,
   step = 0.1,
   steps = 0,
   maxSteps = 101,
   sampleEveryN = 10;

while (steps < maxSteps) {

   if (steps%sampleEveryN === 0) {
       console.log(x + "\t" + y + "\t" + (actual(x) - y).toExponential());
   }
   y = rk4(y, x, step, f);
   // using integer math for the step addition
   // to prevent floating point errors as 0.2 + 0.1 != 0.3
   x = ((x * 10) + (step * 10)) / 10;
   steps += 1;

} </lang>

Output:
0	1			0e+0
1	1.562499854278108	1.4572189210859676e-7
2	3.999999080520799	9.194792007782837e-7
3	10.562497090437551	2.9095624487496252e-6
4	24.999993765090636	6.234909363911356e-6
5	52.562489180302585	1.0819697415342944e-5
6	99.99998340540358	1.659459641700778e-5
7	175.56247648227125	2.3517728749311573e-5
8	288.9999684347986	3.156520142510999e-5
9	451.56245927683966	4.07231603389846e-5
10	675.9999490167097	5.098329029351589e-5

Julia

<lang Julia> function rk4(f)

       return   (t,y,dt)-> 
              ( (dy1   )-> 
              ( (dy2   )-> 
              ( (dy3   )-> 
              ( (dy4   )->( dy1 + 2*dy2 + 2*dy3 + dy4 ) / 6 
              )( dt * f( t +dt  , y + dy3   ) )
              )( dt * f( t +dt/2, y + dy2/2 ) )
              )( dt * f( t +dt/2, y + dy1/2 ) )
              )( dt * f( t      , y         ) )

end

theory(t) = (t^2 + 4.0)^2 / 16.0

tmax = 10.0 ttol = 1.e-5

t0 = 0.0 y0 = 1.0 dt = 0.1

dy = rk4( (t,y) -> t*sqrt(y) )

t = t0 y = y0

while t <= tmax

       if abs(round(t) - t) < ttol 
         @printf( STDOUT,"y(%4.1f)\t= %12.6f \t error: %12.6e\n",t,y,abs(y-theory(t)) )
       end
       y = y + dy(t,y,dt)
       t = t + dt

end </lang>

Output:
y( 0.0)	=     1.000000 	 error: 0.000000e+00
y( 1.0)	=     1.562500 	 error: 1.457219e-07
y( 2.0)	=     3.999999 	 error: 9.194792e-07
y( 3.0)	=    10.562497 	 error: 2.909562e-06
y( 4.0)	=    24.999994 	 error: 6.234909e-06
y( 5.0)	=    52.562489 	 error: 1.081970e-05
y( 6.0)	=    99.999983 	 error: 1.659460e-05
y( 7.0)	=   175.562476 	 error: 2.351773e-05
y( 8.0)	=   288.999968 	 error: 3.156520e-05
y( 9.0)	=   451.562459 	 error: 4.072316e-05
y(10.0)	=   675.999949 	 error: 5.098329e-05

Mathematica

<lang Mathematica>(* Symbolic solution *) DSolve[{y'[t] == t*Sqrt[y[t]], y[0] == 1}, y, t] Table[{t, 1/16 (4 + t^2)^2}, {t, 0, 10}]

(* Numerical solution I (not RK4) *) Table[{t, y[t], Abs[y[t] - 1/16*(4 + t^2)^2]}, {t, 0, 10}] /.

First@NDSolve[{y'[t] == t*Sqrt[y[t]], y[0] == 1}, y, {t, 0, 10}]

(* Numerical solution II (RK4) *) f[{t_, y_}] := {1, t Sqrt[y]} h = 0.1; phi[y_] := Module[{k1, k2, k3, k4},

 k1 = h*f[y];
 k2 = h*f[y + 1/2 k1];
 k3 = h*f[y + 1/2 k2];
 k4 = h*f[y + k3];
 y + k1/6 + k2/3 + k3/3 + k4/6]

solution = NestList[phi, {0, 1}, 101]; Table[{y1, y2, Abs[y2 - 1/16 (y1^2 + 4)^2]},

 {y,  solution1 ;; 101 ;; 10}] 

</lang>

MATLAB

The normally-used built-in solver is the ode45 function, which uses a non-fixed-step solver with 4th/5th order Runge-Kutta methods. The MathWorks Support Team released a package of fixed-step RK method ODE solvers on MATLABCentral. The ode4 function contained within uses a 4th-order Runge-Kutta method. Here is code that tests both ode4 and my own function, shows that they are the same, and compares them to the exact solution. <lang MATLAB>function testRK4Programs

   figure
   hold on
   t = 0:0.1:10;
   y = 0.0625.*(t.^2+4).^2;
   plot(t, y, '-k')
   [tode4, yode4] = testODE4(t);
   plot(tode4, yode4, '--b')
   [trk4, yrk4] = testRK4(t);
   plot(trk4, yrk4, ':r')
   legend('Exact', 'ODE4', 'RK4')
   hold off
   fprintf('Time\tExactVal\tODE4Val\tODE4Error\tRK4Val\tRK4Error\n')
   for k = 1:10:length(t)
       fprintf('%.f\t\t%7.3f\t\t%7.3f\t%7.3g\t%7.3f\t%7.3g\n', t(k), y(k), ...
           yode4(k), abs(y(k)-yode4(k)), yrk4(k), abs(y(k)-yrk4(k)))
   end

end

function [t, y] = testODE4(t)

   y0 = 1;
   y = ode4(@(tVal,yVal)tVal*sqrt(yVal), t, y0);

end

function [t, y] = testRK4(t)

   dydt = @(tVal,yVal)tVal*sqrt(yVal);
   y = zeros(size(t));
   y(1) = 1;
   for k = 1:length(t)-1
       dt = t(k+1)-t(k);
       dy1 = dt*dydt(t(k), y(k));
       dy2 = dt*dydt(t(k)+0.5*dt, y(k)+0.5*dy1);
       dy3 = dt*dydt(t(k)+0.5*dt, y(k)+0.5*dy2);
       dy4 = dt*dydt(t(k)+dt, y(k)+dy3);
       y(k+1) = y(k)+(dy1+2*dy2+2*dy3+dy4)/6;
   end

end</lang>

Output:
Time	ExactVal	ODE4Val		ODE4Error	RK4Val		RK4Error
0	  1.000		  1.000		      0		  1.000		      0
1	  1.563		  1.562		1.46e-007	  1.562		1.46e-007
2	  4.000		  4.000		9.19e-007	  4.000		9.19e-007
3	 10.563		 10.562		2.91e-006	 10.562		2.91e-006
4	 25.000		 25.000		6.23e-006	 25.000		6.23e-006
5	 52.563		 52.562		1.08e-005	 52.562		1.08e-005
6	100.000		100.000		1.66e-005	100.000		1.66e-005
7	175.563		175.562		2.35e-005	175.562		2.35e-005
8	289.000		289.000		3.16e-005	289.000		3.16e-005
9	451.563		451.562		4.07e-005	451.562		4.07e-005
10	676.000		676.000		5.1e-005	676.000		5.1e-005

Maxima

<lang maxima>/* Here is how to solve a differential equation */ 'diff(y, x) = x * sqrt(y); ode2(%, y, x); ic1(%, x = 0, y = 1); factor(solve(%, y)); /* [y = (x^2 + 4)^2 / 16] */

/* The Runge-Kutta solver is builtin */

load(dynamics)$ sol: rk(t * sqrt(y), y, 1, [t, 0, 10, 1.0])$ plot2d([discrete, sol])$

/* An implementation of RK4 for one equation */

rk4(f, x0, y0, x1, n) := block([h, x, y, vx, vy, k1, k2, k3, k4],

  h: bfloat((x1 - x0) / (n - 1)),
  x: x0,
  y: y0,
  vx: makelist(0, n + 1),
  vy: makelist(0, n + 1),
  vx[1]: x0,
  vy[1]: y0,
  for i from 1 thru n do (
     k1: bfloat(h * f(x, y)),
     k2: bfloat(h * f(x + h / 2, y + k1 / 2)),
     k3: bfloat(h * f(x + h / 2, y + k2 / 2)),
     k4: bfloat(h * f(x + h, y + k3)),
     vy[i + 1]: y: y + (k1 + 2 * k2 + 2 * k3 + k4) / 6,
     vx[i + 1]: x: x + h
  ),
  [vx, vy]

)$

[x, y]: rk4(lambda([x, y], x * sqrt(y)), 0, 1, 10, 101)$

plot2d([discrete, x, y])$

s: map(lambda([x], (x^2 + 4)^2 / 16), x)$

for i from 1 step 10 thru 101 do print(x[i], " ", y[i], " ", y[i] - s[i]);</lang>

МК-61/52

ПП	38	П1	ПП	30	П2	ПП	35	П3	2
*	ПП	30	ИП2	ИП3	+	2	*	+	ИП1
+	3	/	ИП7	+	П7	П8	С/П	БП	00
ИП6	ИП5	+	П6	<->	ИП7	+	П8

ИП8	КвКор	ИП6	*

ИП5	*	В/О

Input: 1/2 (h/2) - Р5, 1 (y0) - Р8 and Р7, 0 (t0) - Р6.

OCaml

<lang ocaml>let y' t y = t *. sqrt y let exact t = let u = 0.25*.t*.t +. 1.0 in u*.u

let rk4_step (y,t) h =

 let k1 = h *. y' t y in
 let k2 = h *. y' (t +. 0.5*.h) (y +. 0.5*.k1) in
 let k3 = h *. y' (t +. 0.5*.h) (y +. 0.5*.k2) in
 let k4 = h *. y' (t +. h) (y +. k3) in
 (y +. (k1+.k4)/.6.0 +. (k2+.k3)/.3.0, t +. h)

let rec loop h n (y,t) =

 if n mod 10 = 1 then
   Printf.printf "t = %f,\ty = %f,\terr = %g\n" t y (abs_float (y -. exact t));
 if n < 102 then loop h (n+1) (rk4_step (y,t) h)

let _ = loop 0.1 1 (1.0, 0.0)</lang>

Output:
t = 0.000000,	y = 1.000000,	err = 0
t = 1.000000,	y = 1.562500,	err = 1.45722e-07
t = 2.000000,	y = 3.999999,	err = 9.19479e-07
t = 3.000000,	y = 10.562497,	err = 2.90956e-06
t = 4.000000,	y = 24.999994,	err = 6.23491e-06
t = 5.000000,	y = 52.562489,	err = 1.08197e-05
t = 6.000000,	y = 99.999983,	err = 1.65946e-05
t = 7.000000,	y = 175.562476,	err = 2.35177e-05
t = 8.000000,	y = 288.999968,	err = 3.15652e-05
t = 9.000000,	y = 451.562459,	err = 4.07232e-05
t = 10.000000,	y = 675.999949,	err = 5.09833e-05

Octave

This example is incomplete. Please ensure that it meets all task requirements and remove this message.

Implementing Runge-Kutta in octave is a bit of a hassle. For now we'll showcase how to solve an ordinary differential equation using the lsode function.

<lang octave>function ydot = f(y, t)

   ydot = t * sqrt( y );

endfunction

t = [0:10]'; y = lsode("f", 1, t);

[ t, y, y - 1/16 * (t.**2 + 4).**2 ]</lang>

Output:
ans =

     0.00000     1.00000     0.00000
     1.00000     1.56250     0.00000
     2.00000     4.00000     0.00000
     3.00000    10.56250     0.00000
     4.00000    25.00000     0.00000
     5.00000    52.56250     0.00000
     6.00000   100.00000     0.00000
     7.00000   175.56250     0.00000
     8.00000   289.00001     0.00001
     9.00000   451.56251     0.00001
    10.00000   676.00001     0.00001

PARI/GP

Translation of: C

<lang parigp>rk4(f,dx,x,y)={

 my(k1=dx*f(x,y), k2=dx*f(x+dx/2,y+k1/2), k3=dx*f(x+dx/2,y+k2/2), k4=dx*f(x+dx,y+k3));
 y + (k1 + 2*k2 + 2*k3 + k4) / 6

}; rate(x,y)=x*sqrt(y); go()={

 my(x0=0,x1=10,dx=.1,n=1+(x1-x0)\dx,y=vector(n));
 y[1]=1;
 for(i=2,n,y[i]=rk4(rate, dx, x0 + dx * (i - 1), y[i-1]));
 print("x\ty\trel. err.\n------------");
 forstep(i=1,n,10,
   my(x=x0+dx*i,y2=(x^2/4+1)^2);
   print(x "\t" y[i] "\t" y[i]/y2 - 1)
 )

}; go()</lang>

Output:
x       y       rel. err.
------------
0.100000000     1       -0.00498131231
1.10000000      1.68999982      -0.00383519474
2.10000000      4.40999894      -0.00237694942
3.10000000      11.5599968      -0.00146924588
4.10000000      27.0399933      -0.000961094862
5.10000000      56.2499884      -0.000666538719
6.10000000      106.089982      -0.000485427212
7.10000000      184.959975      -0.000367681962
8.10000000      302.759966      -0.000287408941
9.10000000      470.889955      -0.000230470905

Pascal

Translation of: Ada

This code has been compiled using Free Pascal 2.6.2.

<lang pascal>program RungeKuttaExample;

uses sysutils;

type

   TDerivative = function (t, y : Real) : Real;
   

procedure RungeKutta(yDer : TDerivative;

                    var t, y : array of Real;
                    dt   : Real);

var

   dy1, dy2, dy3, dy4 : Real;
   idx                : Cardinal;

begin

   for idx := Low(t) to High(t) - 1 do
   begin
       dy1 := dt * yDer(t[idx],            y[idx]);
       dy2 := dt * yDer(t[idx] + dt / 2.0, y[idx] + dy1 / 2.0);
       dy3 := dt * yDer(t[idx] + dt / 2.0, y[idx] + dy2 / 2.0);
       dy4 := dt * yDer(t[idx] + dt,       y[idx] + dy3);
       
       t[idx + 1] := t[idx] + dt;
       y[idx + 1] := y[idx] + (dy1 + 2.0 * (dy2 + dy3) + dy4) / 6.0;
   end;

end;

function CalcError(t, y : Real) : Real; var

   trueVal : Real;

begin

   trueVal := sqr(sqr(t) + 4.0) / 16.0;
   CalcError := abs(trueVal - y);

end;

procedure Print(t, y : array of Real;

               modnum : Integer);

var

   idx : Cardinal;

begin

   for idx := Low(t) to High(t) do
   begin
       if idx mod modnum = 0 then
       begin
           WriteLn(Format('y(%4.1f) = %12.8f  Error: %12.6e', 
               [t[idx], y[idx], CalcError(t[idx], y[idx])]));
       end;
   end;

end;

function YPrime(t, y : Real) : Real; begin

   YPrime := t * sqrt(y);

end;

const

   dt = 0.10;
   N = 100;
   

var

   tArr, yArr : array [0..N] of Real;
   

begin

   tArr[0] := 0.0;
   yArr[0] := 1.0;
   
   RungeKutta(@YPrime, tArr, yArr, dt);
   Print(tArr, yArr, 10);

end.</lang>

Output:
y( 0.0) =   1.00000000  Error: 0.00000E+000
y( 1.0) =   1.56249985  Error: 1.45722E-007
y( 2.0) =   3.99999908  Error: 9.19479E-007
y( 3.0) =  10.56249709  Error: 2.90956E-006
y( 4.0) =  24.99999377  Error: 6.23491E-006
y( 5.0) =  52.56248918  Error: 1.08197E-005
y( 6.0) =  99.99998341  Error: 1.65946E-005
y( 7.0) = 175.56247648  Error: 2.35177E-005
y( 8.0) = 288.99996843  Error: 3.15652E-005
y( 9.0) = 451.56245928  Error: 4.07232E-005
y(10.0) = 675.99994902  Error: 5.09833E-005

Perl

There are many ways of doing this. Here we define the runge_kutta function as a function of and , returning a closure which itself takes as argument and returns the next .

Notice how we have to use sprintf to deal with floating point rounding. See perlfaq4. <lang perl>sub runge_kutta {

   my ($yp, $dt) = @_;
   sub {

my ($t, $y) = @_; my @dy = $dt * $yp->( $t , $y ); push @dy, $dt * $yp->( $t + $dt/2, $y + $dy[0]/2 ); push @dy, $dt * $yp->( $t + $dt/2, $y + $dy[1]/2 ); push @dy, $dt * $yp->( $t + $dt , $y + $dy[2] ); return $t + $dt, $y + ($dy[0] + 2*$dy[1] + 2*$dy[2] + $dy[3]) / 6;

   }

}

my $RK = runge_kutta sub { $_[0] * sqrt $_[1] }, .1;

for(

   my ($t, $y) = (0, 1);
   sprintf("%.0f", $t) <= 10;
   ($t, $y) = $RK->($t, $y)

) {

   printf "y(%2.0f) = %12f ± %e\n", $t, $y, abs($y - ($t**2 + 4)**2 / 16)
   if sprintf("%.4f", $t) =~ /0000$/;

}</lang>

Output:
y( 0) =     1.000000 ± 0.000000e+00
y( 1) =     1.562500 ± 1.457219e-07
y( 2) =     3.999999 ± 9.194792e-07
y( 3) =    10.562497 ± 2.909562e-06
y( 4) =    24.999994 ± 6.234909e-06
y( 5) =    52.562489 ± 1.081970e-05
y( 6) =    99.999983 ± 1.659460e-05
y( 7) =   175.562476 ± 2.351773e-05
y( 8) =   288.999968 ± 3.156520e-05
y( 9) =   451.562459 ± 4.072316e-05
y(10) =   675.999949 ± 5.098329e-05

Perl 6

<lang perl6>sub runge-kutta(&yp) {

   return -> \t, \y, \δt {
       my $a = δt * yp( t, y );
       my $b = δt * yp( t + δt/2, y + $a/2 );
       my $c = δt * yp( t + δt/2, y + $b/2 );
       my $d = δt * yp( t + δt, y + $c );
       ($a + 2*($b + $c) + $d) / 6;
   }

}

constant δt = .1; my &δy = runge-kutta { $^t * sqrt($^y) };

loop (

   my ($t, $y) = (0, 1);
   $t <= 10;
   ($t, $y) = ($t + δt, $y + δy($t, $y, δt))

) {

   printf "y(%2d) = %12f ± %e\n", $t, $y, abs($y - ($t**2 + 4)**2 / 16)
   if $t.narrow ~~ Int;

}</lang>

Output:
y( 0) =     1.000000 ± 0.000000e+00
y( 1) =     1.562500 ± 1.457219e-07
y( 2) =     3.999999 ± 9.194792e-07
y( 3) =    10.562497 ± 2.909562e-06
y( 4) =    24.999994 ± 6.234909e-06
y( 5) =    52.562489 ± 1.081970e-05
y( 6) =    99.999983 ± 1.659460e-05
y( 7) =   175.562476 ± 2.351773e-05
y( 8) =   288.999968 ± 3.156520e-05
y( 9) =   451.562459 ± 4.072316e-05
y(10) =   675.999949 ± 5.098329e-05

PL/I

<lang PL/I> Runge_Kutta: procedure options (main); /* 10 March 2014 */

  declare (y, dy1, dy2, dy3, dy4) float (18);
  declare t fixed decimal (10,1);
  declare dt float (18) static initial (0.1);
  y = 1;
  do t = 0 to 10 by 0.1;
     dy1 = dt * ydash(t, y);
     dy2 = dt * ydash(t + dt/2, y + dy1/2);
     dy3 = dt * ydash(t + dt/2, y + dy2/2);
     dy4 = dt * ydash(t + dt,   y + dy3);
     if mod(t, 1.0) = 0 then
        put skip edit('y(', trim(t), ')=', y, ', error = ', abs(y - (t**2 + 4)**2 / 16 ))
                     (3 a, column(9), f(16,10), a, f(13,10));      
     y = y + (dy1 + 2*dy2 + 2*dy3 + dy4)/6;
  end;


ydash: procedure (t, y) returns (float(18));

  declare (t, y) float (18) nonassignable;
  return ( t*sqrt(y) );

end ydash;

end Runge_kutta; </lang> Output:-

y(0.0)=     1.0000000000, error =  0.0000000000
y(1.0)=     1.5624998543, error =  0.0000001457
y(2.0)=     3.9999990805, error =  0.0000009195
y(3.0)=    10.5624970904, error =  0.0000029096
y(4.0)=    24.9999937651, error =  0.0000062349
y(5.0)=    52.5624891803, error =  0.0000108197
y(6.0)=    99.9999834054, error =  0.0000165946
y(7.0)=   175.5624764823, error =  0.0000235177
y(8.0)=   288.9999684348, error =  0.0000315652
y(9.0)=   451.5624592768, error =  0.0000407232
y(10.0)=  675.9999490167, error =  0.0000509833

Python

<lang Python>def RK4(f):

   return lambda t, y, dt: (
           lambda dy1: (
           lambda dy2: (
           lambda dy3: (
           lambda dy4: (dy1 + 2*dy2 + 2*dy3 + dy4)/6
           )( dt * f( t + dt  , y + dy3   ) )

)( dt * f( t + dt/2, y + dy2/2 ) ) )( dt * f( t + dt/2, y + dy1/2 ) ) )( dt * f( t , y ) )

def theory(t): return (t**2 + 4)**2 /16

from math import sqrt dy = RK4(lambda t, y: t*sqrt(y))

t, y, dt = 0., 1., .1 while t <= 10:

   if abs(round(t) - t) < 1e-5:

print("y(%2.1f)\t= %4.6f \t error: %4.6g" % ( t, y, abs(y - theory(t))))

   t, y = t + dt, y + dy( t, y, dt )

</lang>

Output:
y(0.0)	= 1.000000 	 error:    0
y(1.0)	= 1.562500 	 error: 1.45722e-07
y(2.0)	= 3.999999 	 error: 9.19479e-07
y(3.0)	= 10.562497 	 error: 2.90956e-06
y(4.0)	= 24.999994 	 error: 6.23491e-06
y(5.0)	= 52.562489 	 error: 1.08197e-05
y(6.0)	= 99.999983 	 error: 1.65946e-05
y(7.0)	= 175.562476 	 error: 2.35177e-05
y(8.0)	= 288.999968 	 error: 3.15652e-05
y(9.0)	= 451.562459 	 error: 4.07232e-05
y(10.0)	= 675.999949 	 error: 5.09833e-05

Ruby

<lang ruby> def calc_rk4(f)

 return ->(t,y,dt){
        ->(dy1   ){
        ->(dy2   ){
        ->(dy3   ){
        ->(dy4   ){ ( dy1 + 2*dy2 + 2*dy3 + dy4 ) / 6 }.call(
          dt * f.call( t + dt  , y + dy3   ))}.call(
          dt * f.call( t + dt/2, y + dy2/2 ))}.call(
          dt * f.call( t + dt/2, y + dy1/2 ))}.call(
          dt * f.call( t       , y         ))}

end

TIME_MAXIMUM, WHOLE_TOLERANCE = 10.0, 1.0e-5 T_START, Y_START, DT = 0.0, 1.0, 0.10

def my_diff_eqn(t,y) ; t * Math.sqrt(y)  ; end def my_solution(t ) ; (t**2 + 4)**2 / 16  ; end def find_error(t,y) ; (y - my_solution(t)).abs  ; end def is_whole?(t ) ; (t.round - t).abs < WHOLE_TOLERANCE ; end

dy = calc_rk4( ->(t,y){my_diff_eqn(t,y)} )

t, y = T_START, Y_START while t <= TIME_MAXIMUM

 printf("y(%4.1f)\t= %12.6f \t error: %12.6e\n",t,y,find_error(t,y)) if is_whole?(t)
 t, y = t + DT, y + dy.call(t,y,DT)

end </lang>

Output:
y( 0.0)	=     1.000000 	 error: 0.000000e+00
y( 1.0)	=     1.562500 	 error: 1.457219e-07
y( 2.0)	=     3.999999 	 error: 9.194792e-07
y( 3.0)	=    10.562497 	 error: 2.909562e-06
y( 4.0)	=    24.999994 	 error: 6.234909e-06
y( 5.0)	=    52.562489 	 error: 1.081970e-05
y( 6.0)	=    99.999983 	 error: 1.659460e-05
y( 7.0)	=   175.562476 	 error: 2.351773e-05
y( 8.0)	=   288.999968 	 error: 3.156520e-05
y( 9.0)	=   451.562459 	 error: 4.072316e-05
y(10.0)	=   675.999949 	 error: 5.098329e-05

Racket

See Euler method#Racket for implementation of simple general ODE-solver.

The Runge-Kutta method <lang racket> (define (RK4 F δt)

 (λ (t y) 
   (define δy1 (* δt (F t y)))
   (define δy2 (* δt (F (+ t (* 1/2 δt)) (+ y (* 1/2 δy1)))))
   (define δy3 (* δt (F (+ t (* 1/2 δt)) (+ y (* 1/2 δy2)))))
   (define δy4 (* δt (F (+ t δt) (+ y δy1))))
   (list (+ t δt) 
         (+ y (* 1/6 (+ δy1 (* 2 δy2) (* 2 δy3) δy4))))))

</lang>

The method modifier which divides each time-step into n sub-steps: <lang racket> (define ((step-subdivision n method) F h)

 (λ (x . y) (last (ODE-solve F (cons x y) 
                             #:x-max (+ x h) 
                             #:step (/ h n)
                             #:method method))))

</lang>

Usage: <lang racket> (define (F t y) (* t (sqrt y)))

(define (exact-solution t) (* 1/16 (sqr (+ 4 (sqr t)))))

(define numeric-solution

   (ODE-solve F '(0 1) #:x-max 10 #:step 1 #:method (step-subdivision 10 RK4)))

(for ([s numeric-solution])

 (match-define (list t y) s)
 (printf "t=~a\ty=~a\terror=~a\n" t y (- y (exact-solution t))))

</lang>

Output:
t=0	y=1	                error=0
t=1	y=1.562499854278108	error=-1.4572189210859676e-07
t=2	y=3.999999080520799	error=-9.194792007782837e-07
t=3	y=10.562497090437551	error=-2.9095624487496252e-06
t=4	y=24.999993765090636	error=-6.234909363911356e-06
t=5	y=52.562489180302585	error=-1.0819697415342944e-05
t=6	y=99.99998340540358	error=-1.659459641700778e-05
t=7	y=175.56247648227125	error=-2.3517728749311573e-05
t=8	y=288.9999684347986	error=-3.156520142510999e-05
t=9	y=451.56245927683966	error=-4.07231603389846e-05
t=10	y=675.9999490167097	error=-5.098329029351589e-05

Graphical representation:

<lang racket> > (require plot) > (plot (list (function exact-solution 0 10 #:label "Exact solution")

             (points numeric-solution #:label "Runge-Kutta method"))
  #:x-label "t" #:y-label "y(t)")

</lang>

REXX

<lang rexx>/*REXX program uses the Runge-Kutta method to solve the differential */ /* ____ */ /*equation: y'(t)=t²√y(t) which has the exact solution: y(t)=(t²+4)²/16*/

numeric digits 40; d=digits()%2 /*use forty digits, show ½ that. */ x0=0; x1=10; dx=.1; n=1 + (x1-x0) / dx; y.=1

      do m=1  for n-1;                mm=m-1
      y.m=Runge_Kutta(dx, x0+dx*mm, y.mm)
      end   /*m*/

say center(x,13,'─') center(y,d,'─') ' ' center('relative error',d,'─')

      do i=0  to n-1  by 10;         x=(x0+dx*i)/1;       y2=(x*x/4+1)**2
      relE=format(y.i/y2-1,,13)/1;   if relE=0  then relE=' 0'
      say  center(x,13)   right(format(y.i,,12),d)    '  '   left(relE,d)
      end   /*i*/

exit /*stick a fork in it, we're done.*/ /*──────────────────────────────────RATE subroutine─────────────────────*/ rate: return arg(1)*sqrt(arg(2)) /*──────────────────────────────────Runge_Kutta subroutine──────────────*/ Runge_Kutta: procedure; parse arg dx,x,y

                                         k1 = dx * rate(x      , y      )
                                         k2 = dx * rate(x+dx/2 , y+k1/2 )
                                         k3 = dx * rate(x+dx/2 , y+k2/2 )
                                         k4 = dx * rate(x+dx   , y+k3   )

return y + (k1 + 2*k2 + 2*k3 + k4) / 6 /*──────────────────────────────────SQRT subroutine─────────────────────*/ sqrt: procedure; parse arg x; if x=0 then return 0; d=digits()

      numeric digits 11;        g=.sqrtG()
              do j=0 while p>9; m.j=p; p=p%2+1; end;  do k=j+5 to 0 by -1
      if m.k>11  then numeric digits m.k
      g=.5*(g+x/g); end;        numeric digits d;     return g/1

.sqrtG: numeric form; m.=11; p=d+d%4+2

      parse value format(x,2,1,,0) 'E0' with g 'E' _ .; return g*.5'E'_%2</lang>

output

──────X────── ─────────Y──────────   ───relative error───
      0             1.000000000000     0
      1             1.562499854278    -9.3262010934636E-8
      2             3.999999080521    -2.2986980018794E-7
      3            10.562497090438    -2.7546153355698E-7
      4            24.999993765091    -2.4939637458826E-7
      5            52.562489180303    -2.058444217357E-7
      6            99.999983405404    -1.6594596403363E-7
      7           175.562476482271    -1.3395644712797E-7
      8           288.999968434799    -1.092221504E-7
      9           451.562459276840    -9.018277747572E-8
     10           675.999949016709    -7.5419068845528E-8

Run BASIC

<lang Runbasic>y = 1 while t <= 10

  k1	=  t        * sqr(y)
  k2	= (t + .05) * sqr(y + .05 * k1)
  k3	= (t + .05) * sqr(y + .05 * k2)
  k4	= (t + .1)  * sqr(y + .1  * k3)

if right$(using("##.#",t),1) = "0" then print "y(";using("##",t);") ="; using("####.#######", y);chr$(9);"Error ="; (((t^2 + 4)^2) /16) -y

   y = y + .1 *(k1 + 2 * (k2 + k3) + k4) / 6
  t = t + .1

wend end</lang>

Output:
y( 0) =   1.0000000	Error =0
y( 1) =   1.5624999	Error =1.45721892e-7
y( 2) =   3.9999991	Error =9.19479203e-7
y( 3) =  10.5624971	Error =2.90956246e-6
y( 4) =  24.9999938	Error =6.23490939e-6
y( 5) =  52.5624892	Error =1.08196973e-5
y( 6) =  99.9999834	Error =1.65945961e-5
y( 7) = 175.5624765	Error =2.3517728e-5
y( 8) = 288.9999684	Error =3.15652e-5
y( 9) = 451.5624593	Error =4.07231581e-5
y(10) = 675.9999490	Error =5.09832864e-5

Standard ML

<lang sml>fun step y' (tn,yn) dt =

   let
       val dy1 = dt * y'(tn,yn)
       val dy2 = dt * y'(tn + 0.5 * dt, yn + 0.5 * dy1)
       val dy3 = dt * y'(tn + 0.5 * dt, yn + 0.5 * dy2)
       val dy4 = dt * y'(tn + dt, yn + dy3)
   in
       (tn + dt, yn + (1.0 / 6.0) * (dy1 + 2.0*dy2 + 2.0*dy3 + dy4))
   end

(* Suggested test case *) fun testy' (t,y) =

   t * Math.sqrt y

fun testy t =

   (1.0 / 16.0) * Math.pow(Math.pow(t,2.0) + 4.0, 2.0)

(* Test-runner that iterates the step function and prints the results. *) fun test t0 y0 dt steps print_freq y y' =

   let
       fun loop i (tn,yn) =
           if i = steps then ()
           else
               let
                   val (t1,y1) = step y' (tn,yn) dt
                   val y1' = y tn
                   val () = if i mod print_freq = 0 then
                                (print ("Time: " ^ Real.toString tn ^ "\n");
                                 print ("Exact: " ^ Real.toString y1' ^ "\n");
                                 print ("Approx: " ^ Real.toString yn ^ "\n");
                                 print ("Error: " ^ Real.toString (y1' - yn) ^ "\n\n"))
                            else ()
                in
                    loop (i+1) (t1,y1)
               end
   in
       loop 0 (t0,y0)
   end

(* Run the suggested test case *) val () = test 0.0 1.0 0.1 101 10 testy testy'</lang> Output

Time: 0.0
Exact: 1.0
Approx: 1.0
Error: ~1.11022302463E~16

Time: 1.0
Exact: 1.5625
Approx: 1.56249985428
Error: 1.45722452549E~07

Time: 2.0
Exact: 4.0
Approx: 3.99999908052
Error: 9.19479203443E~07

Time: 3.0
Exact: 10.5625
Approx: 10.5624970904
Error: 2.90956245586E~06

Time: 4.0
Exact: 25.0
Approx: 24.9999937651
Error: 6.23490938878E~06

Time: 5.0
Exact: 52.5625
Approx: 52.5624891803
Error: 1.08196973727E~05

Time: 6.0
Exact: 100.0
Approx: 99.9999834054
Error: 1.65945961186E~05

Time: 7.0
Exact: 175.5625
Approx: 175.562476482
Error: 2.35177280956E~05

Time: 8.0
Exact: 289.0
Approx: 288.999968435
Error: 3.15651997767E~05

Time: 9.0
Exact: 451.5625
Approx: 451.562459277
Error: 4.07231581221E~05

Time: 10.0
Exact: 676.0
Approx: 675.999949017
Error: 5.09832866555E~05

Tcl

<lang tcl>package require Tcl 8.5

  1. Hack to bring argument function into expression

proc tcl::mathfunc::dy {t y} {upvar 1 dyFn dyFn; $dyFn $t $y}

proc rk4step {dyFn y* t* dt} {

   upvar 1 ${y*} y ${t*} t
   set dy1 [expr {$dt * dy($t,       $y)}]
   set dy2 [expr {$dt * dy($t+$dt/2, $y+$dy1/2)}]
   set dy3 [expr {$dt * dy($t+$dt/2, $y+$dy2/2)}]
   set dy4 [expr {$dt * dy($t+$dt,   $y+$dy3)}]
   set y [expr {$y + ($dy1 + 2*$dy2 + 2*$dy3 + $dy4)/6.0}]
   set t [expr {$t + $dt}]

}

proc y {t} {expr {($t**2 + 4)**2 / 16}} proc δy {t y} {expr {$t * sqrt($y)}}

proc printvals {t y} {

   set err [expr {abs($y - [y $t])}]
   puts [format "y(%.1f) = %.8f\tError: %.8e" $t $y $err]

}

set t 0.0 set y 1.0 set dt 0.1 printvals $t $y for {set i 1} {$i <= 101} {incr i} {

   rk4step  δy  y t  $dt
   if {$i%10 == 0} {

printvals $t $y

   }

}</lang>

Output:
y(0.0) = 1.00000000	Error: 0.00000000e+00
y(1.0) = 1.56249985	Error: 1.45721892e-07
y(2.0) = 3.99999908	Error: 9.19479203e-07
y(3.0) = 10.56249709	Error: 2.90956245e-06
y(4.0) = 24.99999377	Error: 6.23490939e-06
y(5.0) = 52.56248918	Error: 1.08196973e-05
y(6.0) = 99.99998341	Error: 1.65945961e-05
y(7.0) = 175.56247648	Error: 2.35177280e-05
y(8.0) = 288.99996843	Error: 3.15652000e-05
y(9.0) = 451.56245928	Error: 4.07231581e-05
y(10.0) = 675.99994902	Error: 5.09832864e-05

zkl

Translation of: OCaml

<lang zkl>fcn yp(t,y) { t * y.sqrt() } fcn exact(t){ u:=0.25*t*t + 1.0; u*u }

fcn rk4_step([(y,t)],h){

  k1:=h * yp(t,y);
  k2:=h * yp(t + 0.5*h, y + 0.5*k1);
  k3:=h * yp(t + 0.5*h, y + 0.5*k2);
  k4:=h * yp(t + h, y + k3);
  T(y + (k1+k4)/6.0 + (k2+k3)/3.0, t + h);

}

fcn loop(h,n,[(y,t)]){

  if(n % 10 == 1)
     print("t = %f,\ty = %f,\terr = %g\n".fmt(t,y,(y - exact(t)).abs()));
  if(n < 102) return(loop(h,(n+1),rk4_step(T(y,t),h))) //tail recursion

}</lang>

Output:
loop(0.1,1,T(1.0, 0.0))
t = 0.000000,	y = 1.000000,	err = 0
t = 1.000000,	y = 1.562500,	err = 1.45722e-07
t = 2.000000,	y = 3.999999,	err = 9.19479e-07
t = 3.000000,	y = 10.562497,	err = 2.90956e-06
t = 4.000000,	y = 24.999994,	err = 6.23491e-06
t = 5.000000,	y = 52.562489,	err = 1.08197e-05
t = 6.000000,	y = 99.999983,	err = 1.65946e-05
t = 7.000000,	y = 175.562476,	err = 2.35177e-05
t = 8.000000,	y = 288.999968,	err = 3.15652e-05
t = 9.000000,	y = 451.562459,	err = 4.07232e-05
t = 10.000000,	y = 675.999949,	err = 5.09833e-05