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>
- syntax: GAWK -f RUNGE-KUTTA_METHOD.AWK
- 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>
- 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 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
<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
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
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
Pascal
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
- 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
<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