Runge-Kutta method

Runge-Kutta method is a draft programming task. It is not yet considered ready to be promoted as a complete task, for reasons that should be found in its talk page.

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

Mathematica

<lang Mathematica>DSolve[{y'[t] == t * Sqrt[y[t]], y[0] == 1, y'[0] == 0}, y[t], t]//Simplify -> {{y[t] -> 1/16 (4+t^2)^2}} (* Symbolic expression found : no delta from expected solution *) y[t] -> 1/16*(4+t^2)^2 /. t -> Range[10] ->y[{1,2,3,4,5,6,7,8,9,10}]->{25/16, 4, 169/16, 25, 841/16, 100, 2809/16, 289, 7225/16, 676}</lang>

Python

<lang Python>

2. RHS of the ODE to be solved

def yp(t,y): from math import sqrt return t*sqrt(y)

2. The RK4 routine
4. INPUT ARGUMENTS
5. yp function object providing the RHS of the 1st order ODE
6. tmin, tmax the domain of the solution
7. y0, t0 initial values
8. dt time step over which the solution is sought
10. OUTPUT
11. ti the list of points where the equation has been solved
12. yi the values of the solution found

def RK4(yp,tmin,tmax,y0,t0,dt): yi=[y0] ti=[t0] nmax=int( (tmax-tmin)/dt +1) for n in range(1,nmax,1): tn=ti[n-1] yn=yi[n-1] dy1=dt*yp( tn, yn) dy2=dt*yp( tn+dt/2.0, yn+dy1/2.0) dy3=dt*yp( tn+dt/2.0, yn+dy2/2.0) dy4=dt*yp( tn+dt, yn+dy3) yi.append( yn+(1.0/6.0)*(dy1+2.0*dy2+2.0*dy3+dy4) ) ti.append(tn+dt) return [ti,yi]

[tmin, tmax, dt] = [0.0, 10.0, 0.1] [y0, t0] = [1.0, 0.0] [t,y]=RK4(yp,tmin,tmax,y0,t0,dt) for i in range(0,len(t),10): print ("y(%2.1f)\t= %4.6f \t error:%4.6g")%(t[i], y[i], abs( y[i]- ((t[i]**2 + 4.0)**2)/16.0 ) )

</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

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