Runge-Kutta method
Given the example Differential equation:
![Task](http://static.miraheze.org/rosettacodewiki/thumb/b/ba/Rcode-button-task-crushed.png/64px-Rcode-button-task-crushed.png)
You are encouraged to solve this task according to the task description, using any language you may know.
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
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>
- include <stdlib.h>
- 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
<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
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
<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
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>
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
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
<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
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
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
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)")
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
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