Runge-Kutta method

Given the example Differential equation:

y'(t)=t\times {\sqrt {y(t)}}

With initial condition:

t_{0}=0

and

y_{0}=y(t_{0})=y(0)=1

This equation has an exact solution:

$y(t)={\tfrac {1}{16}}(t^{2}+4)^{2}$
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 $t=0\ldots 10$ with a step value of $\delta t=0.1$ (101 total points, the first being given)
Print the calculated values of $y$ at whole numbered $t$ 's ( $0.0,1.0,\ldots 10.0$ ) along with error as compared to the exact solution.

Method summary

Starting with a given $y_{n}$ and $t_{n}$ calculate:

\delta y_{1}=\delta t\times y'(t_{n},y_{n})\quad

\delta y_{2}=\delta t\times y'(t_{n}+{\tfrac {1}{2}}\delta t,y_{n}+{\tfrac {1}{2}}\delta y_{1})

\delta y_{3}=\delta t\times y'(t_{n}+{\tfrac {1}{2}}\delta t,y_{n}+{\tfrac {1}{2}}\delta y_{2})

\delta y_{4}=\delta t\times y'(t_{n}+\delta t,y_{n}+\delta y_{3})\quad

then:

y_{n+1}=y_{n}+{\tfrac {1}{6}}(\delta y_{1}+2\delta y_{2}+2\delta y_{3}+\delta y_{4})

t_{n+1}=t_{n}+\delta t\quad

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

D

Translation of: Ada

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

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

void runge(in FP function(in FP, in FP) pure nothrow yp_func,

          ref FPs t, ref FPs y, in FP dt) pure nothrow {
   foreach (n; 0 .. t.length - 1) {
       FP dy1 = dt * yp_func(t[n], y[n]);
       FP dy2 = dt * yp_func(t[n] + dt / 2.0, y[n] + dy1 / 2.0);
       FP dy3 = dt * yp_func(t[n] + dt / 2.0, y[n] + dy2 / 2.0);
       FP 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 {

   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 * sqrt(y), t_arr, y_arr, dt);

   foreach (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

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

Tcl

<lang tcl>package require Tcl 8.5

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