Runge-Kutta method: Difference between revisions
Added Easylang
Thundergnat (talk | contribs) (→{{header|C}}: Re add (erroneously?) C example, removed with no explanation) |
(Added Easylang) |
||
(46 intermediate revisions by 25 users not shown) | |||
Line 6:
This equation has an exact solution:
:<math>y(t) = \tfrac{1}{16}(t^2 +4)^2</math>
;Task
Demonstrate the commonly used explicit [[wp:Runge–Kutta_methods#Common_fourth-order_Runge.E2.80.93Kutta_method|fourth-order Runge–Kutta method]] to solve the above differential equation.
* Solve the given differential equation over the range <math>t = 0 \ldots 10</math> with a step value of <math>\delta t=0.1</math> (101 total points, the first being given)
* Print the calculated values of <math>y</math> at whole numbered <math>t</math>'s (<math>0.0, 1.0, \ldots 10.0</math>) along with error as compared to the exact solution.
;Method summary
Starting with a given <math>y_n</math> and <math>t_n</math> calculate:
Line 19 ⟶ 23:
:<math>y_{n+1} = y_n + \tfrac{1}{6} (\delta y_1 + 2\delta y_2 + 2\delta y_3 + \delta y_4)</math>
:<math>t_{n+1} = t_n + \delta t\quad</math>
<br><br>
=={{header|11l}}==
{{trans|Python}}
<syntaxhighlight lang="11l">F rk4(f, x0, y0, x1, n)
V vx = [0.0] * (n + 1)
V vy = [0.0] * (n + 1)
V h = (x1 - x0) / Float(n)
V x = x0
V y = y0
vx[0] = x
vy[0] = y
L(i) 1..n
V k1 = h * f(x, y)
V k2 = h * f(x + 0.5 * h, y + 0.5 * k1)
V k3 = h * f(x + 0.5 * h, y + 0.5 * k2)
V k4 = h * f(x + h, y + k3)
vx[i] = x = x0 + i * h
vy[i] = y = y + (k1 + k2 + k2 + k3 + k3 + k4) / 6
R (vx, vy)
F f(Float x, Float y) -> Float
R x * sqrt(y)
V (vx, vy) = rk4(f, 0.0, 1.0, 10.0, 100)
L(x, y) zip(vx, vy)[(0..).step(10)]
print(‘#2.1 #4.5 #2.8’.format(x, y, y - (4 + x * x) ^ 2 / 16))</syntaxhighlight>
{{out}}
<pre>
0.0 1.00000 0.00000000
1.0 1.56250 -1.45721892e-7
2.0 4.00000 -9.194792e-7
3.0 10.56250 -0.00000291
4.0 24.99999 -0.00000623
5.0 52.56249 -0.00001082
6.0 99.99998 -0.00001659
7.0 175.56248 -0.00002352
8.0 288.99997 -0.00003157
9.0 451.56246 -0.00004072
10.0 675.99995 -0.00005098
</pre>
=={{header|Action!}}==
Calculations on a real Atari 8-bit computer take quite long time. It is recommended to use an emulator capable with increasing speed of Atari CPU.
{{libheader|Action! Tool Kit}}
{{libheader|Action! Real Math}}
<syntaxhighlight lang="action!">INCLUDE "D2:PRINTF.ACT" ;from the Action! Tool Kit
INCLUDE "H6:REALMATH.ACT"
DEFINE PTR="CARD"
REAL one,two,four,six
PROC Init()
IntToReal(1,one)
IntToReal(2,two)
IntToReal(4,four)
IntToReal(6,six)
RETURN
PROC Fun=*(REAL POINTER x,y,res)
DEFINE JSR="$20"
DEFINE RTS="$60"
[JSR $00 $00 ;JSR to address set by SetFun
RTS]
PROC SetFun(PTR p)
PTR addr
addr=Fun+1 ;location of address of JSR
PokeC(addr,p)
RETURN
PROC Rate(REAL POINTER x,y,res)
REAL tmp
Sqrt(y,tmp) ;tmp=sqrt(y)
RealMult(x,tmp,res) ;res=x*sqrt(y)
RETURN
PROC RK4(PTR f REAL POINTER dx,x,y,res)
REAL k1,k2,k3,k4,dx2,k12,k22,tmp1,tmp2,tmp3
SetFun(f)
Fun(x,y,tmp1) ;tmp1=f(x,y)
RealMult(dx,tmp1,k1) ;k1=dx*f(x,y)
RealDiv(dx,two,dx2) ;dx2=dx/2
RealDiv(k1,two,k12) ;k12=k1/2
RealAdd(x,dx2,tmp1) ;tmp1=x+dx/2
RealAdd(y,k12,tmp2) ;tmp2=y+k1/2
Fun(tmp1,tmp2,tmp3) ;tmp3=f(x+dx/2,y+k1/2)
RealMult(dx,tmp3,k2) ;k2=dx*f(x+dx/2,y+k1/2)
RealDiv(k2,two,k22) ;k22=k2/2
RealAdd(y,k22,tmp2) ;tmp2=y+k2/2
Fun(tmp1,tmp2,tmp3) ;tmp3=f(x+dx/2,y+k2/2)
RealMult(dx,tmp3,k3) ;k3=dx*f(x+dx/2,y+k2/2)
RealAdd(x,dx,tmp1) ;tmp1=x+dx
RealAdd(y,k3,tmp2) ;tmp2=y+k3
Fun(tmp1,tmp2,tmp3) ;tmp3=f(x+dx,y+k3)
RealMult(dx,tmp3,k4) ;k4=dx*f(x+dx,y+k3)
RealAdd(k2,k3,tmp1) ;tmp1=k2+k3
RealMult(two,tmp1,tmp2) ;tmp2=2*k2+2*k3
RealAdd(k1,tmp2,tmp1) ;tmp3=k1+2*k2+2*k3
RealAdd(tmp1,k4,tmp2) ;tmp2=k1+2*k2+2*k3+k4
RealDiv(tmp2,six,tmp1) ;tmp1=(k1+2*k2+2*k3+k4)/6
RealAdd(y,tmp1,res) ;res=y+(k1+2*k2+2*k3+k4)/6
RETURN
PROC Calc(REAL POINTER x,res)
REAL tmp1,tmp2
RealMult(x,x,tmp1) ;tmp1=x*x
RealDiv(tmp1,four,tmp2) ;tmp2=x*x/4
RealAdd(tmp2,one,tmp1) ;tmp1=x*x/4+1
Power(tmp1,two,res) ;res=(x*x/4+1)^2
RETURN
PROC RelError(REAL POINTER a,b,res)
REAL tmp
RealDiv(a,b,tmp) ;tmp=a/b
RealSub(tmp,one,res) ;res=a/b-1
RETURN
PROC Main()
REAL x0,x1,x,dx,y,y2,err,tmp1,tmp2
CHAR ARRAY s(20)
INT i,n
Put(125) PutE() ;clear the screen
MathInit()
Init()
PrintF("%-2S %-11S %-8S%E","x","y","rel err")
IntToReal(0,x0)
IntToReal(10,x1)
ValR("0.1",dx)
RealSub(x1,x0,tmp1) ;tmp1=x1-x0
RealDiv(tmp1,dx,tmp2) ;tmp2=(x1-x0)/dx
n=RealToInt(tmp2) ;n=(x1-x0)/dx
i=0
IntToReal(1,y)
DO
IntToReal(i,tmp1) ;tmp1=i
RealMult(dx,tmp1,tmp2) ;tmp2=i*dx
RealAdd(x0,tmp2,x) ;x=x0+i*dx
IF i MOD 10=0 THEN
Calc(x,y2)
RelError(y,y2,err)
StrR(x,s) PrintF("%-2S ",s)
StrR(y,s) PrintF("%-11S ",s)
StrR(err,s) PrintF("%-8S%E",s)
FI
i==+1
IF i>n THEN EXIT FI
RK4(rate,dx,x,y,tmp1) ;tmp1=rk4(rate,dx,x0+dx*(i-1),y)
RealAssign(tmp1,y) ;y=rk4(rate,dx,x0+dx*(i-1),y)
OD
RETURN</syntaxhighlight>
{{out}}
[https://gitlab.com/amarok8bit/action-rosetta-code/-/raw/master/images/Runge-Kutta_method.png Screenshot from Atari 8-bit computer]
<pre>
x y rel err
0 1 0
1 1.56249977 -1.3E-07
2 3.99999882 -2.9E-07
3 10.56249647 -2.9E-07
4 24.99999228 -2.9E-07
5 52.56248607 -2.0E-07
6 99.99997763 -2.1E-07
7 175.562459 -1.8E-07
8 288.999935 -1.9E-07
9 451.562406 0
10 675.999869 -1.4E-07
</pre>
=={{header|Ada}}==
<
with Ada.Numerics.Generic_Elementary_Functions;
procedure RungeKutta is
Line 72 ⟶ 261:
Runge (yprime'Access, t_arr, y_arr, dt);
Print (t_arr, y_arr, 10);
end RungeKutta;</
{{out}}
<pre>y(0.0) = 1.00000000 Error: 0.00000E+00
Line 85 ⟶ 274:
y(9.0) = 451.56245928 Error: 4.07232E-05
y(10.0) = 675.99994902 Error: 5.09833E-05</pre>
=={{header|ALGOL 68}}==
<syntaxhighlight lang="algol68">
BEGIN
PROC rk4 = (PROC (REAL, REAL) REAL f, REAL y, x, dx) REAL :
Line 113 ⟶ 303:
OD
END
</syntaxhighlight>
{{out}}
<pre>
Line 128 ⟶ 318:
9.0000000 451.5625000 451.5624593 -9.0183e-08
10.0000000 676.0000000 675.9999490 -7.5419e-08
</pre>
=={{header|ALGOL W}}==
{{Trans|ALGOL 68}}
As originally defined, the signature of a procedure parameter could not be specified in Algol W (as here), modern compilers may require parameter specifications for the "f" parameter of rk4.
<syntaxhighlight lang="algolw">begin
real procedure rk4 ( real procedure f ; real value y, x, dx ) ;
begin % Fourth-order Runge-Kutta method %
real dy1, dy2, dy3, dy4;
dy1 := dx * f(x, y);
dy2 := dx * f(x + dx / 2.0, y + dy1 / 2.0);
dy3 := dx * f(x + dx / 2.0, y + dy2 / 2.0);
dy4 := dx * f(x + dx, y + dy3);
y + (dy1 + 2.0 * dy2 + 2.0 * dy3 + dy4) / 6.0
end rk4;
real x0, x1, y0, dx;
integer numPoints;
x0 := 0; x1 := 10; y0 := 1.0; % Boundary conditions. %
dx := 0.1; % Step size. %
numPoints := entier ((x1 - x0) / dx + 0.5); % Add 0.5 for rounding errors. %
begin
real procedure dyByDx ( real value x, y ) ; x * sqrt(y); % Differential equation. %
real array y ( 0 :: numPoints); y(0) := y0; % Grid and starting point. %
for i := 1 until numPoints do y(i) := rk4 (dyByDx, y(i-1), x0 + dx * (i - 1), dx);
write( " x true y calc y relative error" );
for i := 0 step 10 until numPoints do begin
real x, trueY;
x := x0 + dx * i;
trueY := (x * x + 4.0) ** 2 / 16.0;
write( r_format := "A", r_w := 12, r_d := 7, s_w := 3, x, trueY, y( i )
, r_format := "S", r_w := 12, y( i ) / trueY - 1
)
end for_i
end
end.</syntaxhighlight>
{{out}}
<pre>
x true y calc y relative error
0.0000000 1.0000000 1.0000000 0.0000e+000
1.0000000 1.5625000 1.5624998 -9.3262e-008
2.0000000 4.0000000 3.9999990 -2.2986e-007
3.0000000 10.5625000 10.5624971 -2.7546e-007
4.0000000 25.0000000 24.9999937 -2.4939e-007
5.0000000 52.5625000 52.5624891 -2.0584e-007
6.0000000 100.0000000 99.9999834 -1.6594e-007
7.0000000 175.5625000 175.5624764 -1.3395e-007
8.0000000 289.0000000 288.9999684 -1.0922e-007
9.0000000 451.5625000 451.5624592 -9.0182e-008
10.0000000 676.0000000 675.9999490 -7.5419e-008
</pre>
=={{header|APL}}==
<syntaxhighlight lang="apl">
∇RK4[⎕]∇
∇
Line 153 ⟶ 392:
[2] ⎕←'T' 'RK4 Y' 'ERROR'⍪TABLE,TABLE[;2]-{((4+⍵*2)*2)÷16}TABLE[;1]
∇
</syntaxhighlight>
{{out}}
<pre>
Line 172 ⟶ 411:
=={{header|AWK}}==
<syntaxhighlight lang="awk">
# syntax: GAWK -f RUNGE-KUTTA_METHOD.AWK
# converted from BBC BASIC
Line 192 ⟶ 431:
exit(0)
}
</syntaxhighlight>
{{out}}
<pre>
Line 210 ⟶ 449:
=={{header|BASIC}}==
==={{header|BASIC256}}===
<syntaxhighlight lang="basic256">y = 1
for i = 0 to 100
t = i / 10
if t = int(t) then
actual = ((t ^ 2 + 4) ^ 2) / 16
print "y("; int(t); ") = "; left(string(y), 13), "Error = "; left(string(actual - y), 13)
end if
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 = y + 0.1 * (k1 + 2 * (k2 + k3) + k4) / 6
next i
end</syntaxhighlight>
==={{header|BBC BASIC}}===
<
FOR i% = 0 TO 100
t = i% / 10
Line 225 ⟶ 483:
k4 = (t + 0.10) * SQR(y + 0.10 * k3)
y += 0.1 * (k1 + 2 * (k2 + k3) + k4) / 6
NEXT i%</
{{out}}
<pre>y(0) = 1 Error = 0
Line 241 ⟶ 499:
==={{header|IS-BASIC}}===
<
110 LET Y=1
120 FOR T=0 TO 10 STEP .1
Line 250 ⟶ 508:
170 LET K4=(T+.1)*SQR(Y+.1*K3)
180 LET Y=Y+.1*(K1+2*(K2+K3)+K4)/6
190 NEXT</
==={{header|QBasic}}===
{{works with|QBasic|1.1}}
{{works with|QuickBasic|4.5}}
<syntaxhighlight lang="qbasic">y! = 1
FOR i = 0 TO 100
t = i / 10
IF t = INT(t) THEN
actual! = ((t ^ 2 + 4) ^ 2) / 16
PRINT USING "y(##) = ###.###### Error = "; t; y;
PRINT actual - y
END IF
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)
y = y + .1 * (k1 + 2 * (k2 + k3) + k4) / 6
NEXT i</syntaxhighlight>
==={{header|True BASIC}}===
{{works with|QBasic}}
<syntaxhighlight lang="qbasic">LET y = 1
FOR i = 0 TO 100
LET t = i / 10
IF t = INT(t) THEN
LET actual = ((t ^ 2 + 4) ^ 2) / 16
PRINT "y("; STR$(t); ") ="; y ; TAB(20); "Error = "; actual - y
END IF
LET k1 = t * SQR(y)
LET k2 = (t + 0.05) * SQR(y + 0.05 * k1)
LET k3 = (t + 0.05) * SQR(y + 0.05 * k2)
LET k4 = (t + 0.10) * SQR(y + 0.10 * k3)
LET Y = Y + 0.1 * (k1 + 2 * (k2 + k3) + k4) / 6
NEXT i
END</syntaxhighlight>
=={{header|C}}==
<
#include <stdlib.h>
#include <math.h>
Line 289 ⟶ 587:
return 0;
}</
{{out}} (errors are relative)
<pre>
Line 307 ⟶ 605:
</pre>
=={{header|C sharp|C#}}==
<
using System;
Line 402 ⟶ 700:
}
}
}</
=={{header|C++}}==
Using Lambdas
<
* compiled with
* g++
*
* g++ -std=c++14 -o rk4 %
*
*/
# include <iostream>
# include <math.h>
auto rk4(double f(double, double))
{
return [f](double t, double y, double dt) -> double {
double dy1 { dt
return ( dy1 + 2*dy2 + 2*dy3 + dy4 ) / 6;
};
}
int main(void)
{
constexpr
double TIME_MAXIMUM { 10.0 },
T_START { 0.0 },
Y_START { 1.0 },
DT { 0.1 },
WHOLE_TOLERANCE { 1e-12 };
for (
auto y { Y_START }, t { T_START };
t <= TIME_MAXIMUM;
y
)
if (ceilf(t)-t < WHOLE_TOLERANCE)
printf("y(%4.1f)\t=%12.6f \t error: %12.6e\n", t, y, std::fabs(y - pow(t*t+4,2)/16));
return 0;
}</syntaxhighlight>
=={{header|Common Lisp}}==
<
(let ((h (float (/ (- x-end x) n) 1d0))
k1 k2 k3 k4)
Line 484 ⟶ 780:
(7.999999999999988d0 288.9999684347983d0 -3.156520000402452d-5)
(8.999999999999984d0 451.56245927683887d0 -4.072315812209126d-5)
(9.99999999999998d0 675.9999490167083d0 -5.0983286655537086d-5))</
=={{header|Crystal}}==
{{trans|Run Basic and Ruby output}}
<syntaxhighlight lang="ruby">y, t = 1, 0
while t <= 10
k1 = t * Math.sqrt(y)
k2 = (t + 0.05) * Math.sqrt(y + 0.05 * k1)
k3 = (t + 0.05) * Math.sqrt(y + 0.05 * k2)
k4 = (t + 0.1) * Math.sqrt(y + 0.1 * k3)
printf("y(%4.1f)\t= %12.6f \t error: %12.6e\n", t, y, (((t**2 + 4)**2 / 16) - y )) if (t.round - t).abs < 1.0e-5
y += 0.1 * (k1 + 2 * (k2 + k3) + k4) / 6
t += 0.1
end</syntaxhighlight>
{{out}}
<pre>
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
</pre>
=={{header|D}}==
{{trans|Ada}}
<
alias FP = real;
Line 525 ⟶ 850:
t_arr[i], y_arr[i],
calc_err(t_arr[i], y_arr[i]));
}</
{{out}}
<pre>y(0.0) = 1.00000000 Error: 0
Line 540 ⟶ 865:
=={{header|Dart}}==
<
num RungeKutta4(Function f, num t, num y, num dt){
Line 566 ⟶ 891:
t += dt;
}
}</
{{out}}
<pre>
Line 580 ⟶ 905:
y(9.00) = 451.56245928 Error = 9.0182772312e-8
y(10.0) = 675.99994902 Error = 7.5419063100e-8
</pre>
=={{header|EasyLang}}==
{{trans|BASIC256}}
<syntaxhighlight>
numfmt 6 0
y = 1
for i = 0 to 100
t = i / 10
if t = floor t
h = t * t + 4
actual = h * h / 16
print "y(" & t & ") = " & y & " Error = " & 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
.
</syntaxhighlight>
=={{header|EDSAC order code}}==
The EDSAC subroutine library had two Runge-Kutta subroutines: G1 for 35-bit values and G2 for 17-bit values. A demo of G1 is given here. Setting up the parameters is rather complicated, but after that it's just a matter of calling G1 once for every step in the Runge-Kutta process.
Since EDSAC real numbers are restricted to -1 <= x < 1, the values in the Rosetta Code task have to be scaled down. For comparison with other languages it's convenient to divide the y values by 1000. With 100 steps, a convenient time interval is 1/128.
G1 can solve equations in several variables, say y_1, ..., y_n. The user must provide an auxiliary subroutine which calculates dy_1/dt, ..., dy_n/dt from y_1, ..., y_n. If the derivatives also depend on t (as in the Rosetta Code task) it's necessary to add a dummy y variable which is identical with t.
<syntaxhighlight lang="edsac">
[Demo of EDSAC library subroutine G1: Runge-Kutta solution of differential equations.
Full description is in Wilkes, Wheeler & Gill, 1951 edn, pages 32-34, 86-87, 132-134.
Before using G1, we need to fix n, m, a, b, c, d, as defined in WWG pages 86-87:
n = number of equations (2 for the Rosetta Code example).
2^m = multiplier for the hy', as large as possible without causing numeric overflow;
with the scaling chosen here, m = 5.
Variables y are stored in n consecutive long locations, the last of which is aD.
Scaled derivatives (2^m)hy' in n consecutive long locations, the last of which is bD.
G1 uses working variables in n consecutive long locations, the last of which is cD.
d = address of user-supplied auxiliary subroutine, which calculates the (2^m)hy'.
For convenience, keep G1 and its storage together. Start at (say) 400 and place:
variables y at 400D, 402D;
scaled derivatives at 404D, 406D;
workspace for G1 at 408D, 410D;
G1 itself at 412.
If the base address is placed in location 51 at load time, all the above
addresses can be accessed via the G parameter:]
T 51 K
P 400 F
[Now set up the 6 preset parameters specified in WWG:]
T 45 K
P 2#G [H parameter: P a D]
P 4 F [N parameter: P 2n F]
P 4 F [M parameter: P (b-a) F, or V (2048-a+b) F if a > b]
P 4 F [& parameter: P (c-b) F, or V (2048-b+c) F if b > c]
P 8 F [L parameter: P 2^(m-2) F]
P 300 F [X parameter: P d F]
[For other addresses in the program we can optionally use some more parameters:]
T 52 K
P 120 F [A parameter: main routine]
P 56 F [B parameter: print subroutine P1 from EDSAC library]
P 350 F [C parameter: constants for Rosetta code example]
P 78 F [V parameter: square root subroutine]
[Library subroutine to read constants; runs at load time and is then overwritten.
R5, for decimal fractions, seems to be unavailable (lost?), so the values are
here read in as 35-bit integers (i.e. times 2^34) by R2.
Values are: 0.001, initial value of y
(2^23)/(10^7) and 25/(2^10) for use in calculations
0.5/(10^9) for rounding to 9 d.p. (print routine P1 doesn't do this)]
GKT20FVDL8FA40DUDTFI40FA40FS39FG@S2FG23FA5@T5@E4@E13Z
T#C
17179869F14411518808F419430400F9#
TZ
[Library subroutine M3; prints header at load time and is then overwritten.]
PFGKIFAFRDLFUFOFE@A6FG@E8FEZPF
*SCALED!FOR!EDSAC@&!!TIME!!!!!!!!!Y!VIA!RK!!!!!Y!DIRECT@&
....PK [end text with some blank tape]
[Runge-Kutta: auxiliary subroutine to calculate (2^m)*h*(dy1/dt) and (2^m)*h*(dy2/dt)
from y1, y2, where y1 is the function y in Rosetta Code (but scaled) and y2 = t.
For the Rosetta code example we're using m = 5, h = 2^(-7)]
E25K TX GK
A3F T20@ [set up return as usual]
H2#G V2#G TD [acc := t^2, temp store in 0D]
H#G VD LD YF TD [y1 times t^2, shift left, round, temp store in 0D]
H2#C VD YF T4D [times (2^23)/(10^7), round, to 4D for square root]
[14] A14@ GV A4D T4#G [call square root, result in 4D, copy to (2^m)hy']
A21@ T6#G [1/4, i.e. (2^m)h with m and h as above, to (2^m)ht']
[20] ZF [overwritten by jump back to caller]
[21] RF [constant 1/4]
[Main routine, with two subroutines in the same address block as the main routine.]
E25K TA GK
[0] #F [figures shift on teleprinter]
[1] MF [decimal point (in figures mode)]
[2] !F @F &F [space, carriage return, line feed,]
[5] K4096F [null char]
[6] P100F [constant: nr of Runge-Kutta steps (in address field)]
[7] PF [negative count of Runge-Kutta steps]
[8] P10F [constant: number of steps between printed values]
[9] PF [negative count of steps between printed values]
[Enter with acc = 0]
[10] O@ [set teleprinter to figures]
S6@ T7@ [init negative count of R-K steps]
S8@ T9@ [init negative count of print steps]
[Before using library subroutine G1, clear its working registers (WWG page 33)]
T8#G T10#G
[Set up initial values of y1 and y2 (where y2 = t)]
A#C T#G [load 0.001 from constants section, store in y1]
T2#G [y2 = t = 0]
[20] A20@ G40@ [call subroutine to print initial values]
[Loop round Runge-Kutta steps]
[22] TF A23@ G12G [clear accumulator, call G1 for Runge-Kutta step]
A9@ A2F U9@ [update negative print count]
G33@ [skip printing if not reached 0]
S8@ T9@ [reset negative print count]
A31@ G40@ [call subroutine to print values]
[33] TF [clear accumulator]
A7@ A2F U7@ [increment negative count of Runge-Kutta steps]
G22@ [loop till count = 0]
O5@ ZF [flush teleprinter buffer; stop]
[Subroutine to print y1 as calculated (1) by Runge-Kutta (2) direct from formula]
[40] A3F T71@ [set up return as usual]
A2#G TD [latest t (= y2) from Runge-Kutta, to 0D for printing]
[44] A44@ G72@ [call subroutine to print t]
O2@ O2@ [followed by 2 spaces]
A#G TD [latest y1 from Runge-Kutta, to 0D for printing]
[50] A50@ G72@ [call subroutine to print y1]
O2@ O2@ [followed by 2 spaces]
A 4#C [load constant 25/(2^10)]
H2#G V2#G TD [add t^2, temp store result in 0D]
HD VD LD YF TD [square, shift 1 left, round, result to 0D]
H2#C VD YF TD [times (2^23)/(10^7), round, to 0D for printing]
[67] A67@ G72@ [call subroutine to print y]
O3@ O4@ [print CR, LF]
[71] ZF [overwritten by jump back to caller]
[Second-level subroutine to print number in 0D to 9 decimal places]
[72] A3F T82@ [set up return as usual]
AD A6#C TD [load number, add decimal rounding, to 0D for printing]
O81@ O1@ [print '0.' since P1 doesn't do so]
A79@ GB [call library subroutine P1 for printing]
[81] P9F [parameter for P1, 9 decimals]
[82] ZF [overwritten by jump back to caller]
[Library subroutine G1 for Runge-Kutta process. 66 locations, even address.]
E25K T12G
GKT4#ZH682DT6#ZPNT12#Z!1405DT14#ZTHT16#ZT2HTZA3FT61@A31@G63@&FT6ZPN
T8ZMMO&H4@A20@E23@T14ZAHT16ZA2HT18ZH12#@S12#@T12#@E28@H4#@T4DUFS38@
A25@T38@S6#@A16#@U46#@A8@U37@A9@U55@A24@T39@ZFR1057#@ZFYFU6DV6DRLYF
UDZFZFADLDADLLS6DN4DYFZFA46#@S14#@G29@A65@S11@ZFA35@U65@GXZF
[Replacement for library routine S2 (square root). 38 locations, even address.
Advantages: More accurate for small values of the argument.
Calculates sqrt(0) without going into an infinite loop.
Disadvantages: Longer and slower than S2 (calculates one bit at a time).]
E25K TV
GKA3FT31@A4DG32@A33@T36#@T4DA33@RDU34#@RDS4DS33@A36#@G22@T36#@A4DS34#@
T4DA36#@A33@G25@TFA36#@S33@A36#@T36#@A34#@RDYFG9@ZFZFK4096FPFPFPFPF
[Library subroutine P1 - print a single positive number. 21 locations.
Prints number in 0D to n places of decimals, where
n is specified by 'P n F' pseudo-order after subroutine call.]
E25K TB
GKA18@U17@S20@T5@H19@PFT5@VDUFOFFFSFL4FTDA5@A2FG6@EFU3FJFM1F
[Define entry point in main routine]
E25K TA GK
E10Z PF [enter at relative address 10 with accumulator = 0]
</syntaxhighlight>
{{out}}
<pre>
SCALED FOR EDSAC
TIME Y VIA RK Y DIRECT
0.000000000 0.001000000 0.001000000
0.078125000 0.001562499 0.001562500
0.156250000 0.003999998 0.004000000
0.234375000 0.010562495 0.010562500
0.312500000 0.024999992 0.025000000
0.390625000 0.052562487 0.052562500
0.468750000 0.099999981 0.100000000
0.546875000 0.175562474 0.175562500
0.625000000 0.288999965 0.289000000
0.703125000 0.451562456 0.451562500
0.781250000 0.675999945 0.676000000
</pre>
=={{header|ERRE}}==
<syntaxhighlight lang="erre">
PROGRAM RUNGE_KUTTA
Line 608 ⟶ 1,122:
Y+=DELTA_T*(K1+2*(K2+K3)+K4)/6
END FOR
END PROGRAM</
{{out}}
<pre>
Line 624 ⟶ 1,138:
</pre>
=={{header|
<syntaxhighlight lang="Excel">
//Worksheet formula to manage looping
=LET(
T₊, SEQUENCE(11, 1, 0, 1),
T, DROP(T₊, -1),
τ, SEQUENCE(1 / δt, 1, 0, δt),
calculated, SCAN(1, T, LAMBDA(y₀, t, REDUCE(y₀, t + τ, RungaKutta4λ(Dλ)))),
calcs, VSTACK(1, calculated),
exact, f(T₊),
HSTACK(T₊, calcs, exact, (exact - calcs) / exact)
)
//Lambda function passed to RungaKutta4λ to evaluate derivatives
Dλ(y,t)
= LAMBDA(y,t, t * SQRT(y))
//Curried Lambda function with derivative function D and y, t as parameters
RungaKutta4λ(Dλ)
= LAMBDA(D,
LAMBDA(yᵣ, tᵣ,
LET(
δy₁, δt * D(yᵣ, tᵣ),
δy₂, δt * D(yᵣ + δy₁ / 2, tᵣ + δt / 2),
δy₃, δt * D(yᵣ + δy₂ / 2, tᵣ + δt / 2),
δy₄, δt * D(yᵣ + δy₃, tᵣ + δt),
yᵣ₊₁, yᵣ + (δy₁ + 2 * δy₂ + 2 * δy₃ + δy₄) / 6,
yᵣ₊₁
)
)
)
//Lambda function returning the exact solution
f(t)
= LAMBDA(t, (1/16) * (t^2 + 4)^2 )
</syntaxhighlight>
{{out}}
<pre>
Time Calculated Exact Rel Error
0.00 1.000000 1.000000 0.00E+00
1.00 1.562500 1.562500 9.33E-08
2.00 3.999999 4.000000 2.30E-07
3.00 10.562497 10.562500 2.75E-07
4.00 24.999994 25.000000 2.49E-07
5.00 52.562489 52.562500 2.06E-07
6.00 99.999983 100.000000 1.66E-07
7.00 175.562476 175.562500 1.34E-07
8.00 288.999968 289.000000 1.09E-07
9.00 451.562459 451.562500 9.02E-08
10.00 675.999949 676.000000 7.54E-08
</pre>
=={{header|F_Sharp|F#}}==
{{works with|F# interactive (fsi.exe)}}
<syntaxhighlight lang="fsharp">
open System
Line 647 ⟶ 1,218:
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) )</
{{out}}
Line 665 ⟶ 1,236:
=={{header|Fortran}}==
<
implicit none
integer, parameter :: dp = kind(1d0)
real(
real(dp) :: 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)**2/16)
do while (t < tstop)
k1 = dt*f(t, y)
k2 = dt*f(t+dt/2, y+k1/2)
k3 = dt*f(t+dt/2, y+k2/2)
k4 = dt*f(t+dt, y+k3)
y = y+(k1+2*(k2+k3)+k4)/6
t = t+dt
if (abs(nint(t)-t) <= 1d-12) then
write (6, '(a,f4.1,a,f12.8,a,es13.6)') 'y(', t, ') = ', y, ' error = ', &
abs(y-(t**2+4)**2/16)
end if
end do
contains
function f(t,y)
real(dp), intent(in) :: t, y
real(dp) :: f
f = t*sqrt(y)
end function f
end program rungekutta</syntaxhighlight>
{{out}}
<pre>
Line 709 ⟶ 1,283:
y(10.0) = 675.99994902 Error = 5.098329E-05
</pre>
=={{header|FreeBASIC}}==
{{trans|BBC BASIC}}
<
' compile with: fbc -s console
' translation of BBC BASIC
Line 741 ⟶ 1,316:
Print : Print "hit any key to end program"
Sleep
End</
{{out}}
<pre>y(0) = 1 Error = 0
Line 756 ⟶ 1,331:
=={{header|FutureBasic}}==
<
def fn dydx( x as double, y as double ) as double = x * sqr(y)
def fn exactY( x as long ) as double = ( x ^2 + 4 ) ^2 / 16
long i
double h, k1, k2, k3, k4, x, y, result
h = 0.1
y = 1
for i = 0 to 100
x = i * h
if x == int(x)
result = fn exactY( x )
print "y("; mid$( str$(x), 2, len$(str$(x) )); ") = "; y, "Error = "; result - y
end if
k1 = h * fn dydx( x, y )
k2 = h * fn dydx( x + h / 2, y + k1 / 2 )
k3 = h * fn dydx( x + h / 2, y + k2 / 2 )
k4 = h * fn dydx( x + h, y + k3 )
y = y + 1 / 6 * ( k1 + 2 * k2 + 2 * k3 + k4 )
next
HandleEvents</syntaxhighlight>
Output:
<pre>
Line 804 ⟶ 1,374:
=={{header|Go}}==
{{works with|Go1}}
<
import (
Line 860 ⟶ 1,430:
func printErr(t, y float64) {
fmt.Printf("y(%.1f) = %f Error: %e\n", t, y, math.Abs(actual(t)-y))
}</
{{out}}
<pre>
Line 877 ⟶ 1,447:
=={{header|Groovy}}==
<syntaxhighlight lang="groovy">
class Runge_Kutta{
static void main(String[] args){
Line 902 ⟶ 1,472:
static def dy(def x){return x*0.1;}
}
</syntaxhighlight>
{{out}}
<pre>
Line 916 ⟶ 1,486:
y(9.0)=451.56245927683966 Error:4.07231603389846E-5
y(10.0)=675.9999490167097 Error:5.098329029351589E-5
</pre>
=={{header|Hare}}==
<syntaxhighlight lang="hare">use fmt;
use math;
export fn main() void = {
rk4_driver(&f, 0.0, 10.0, 1.0, 0.1);
};
fn rk4_driver(func: *fn(_: f64, _: f64) f64, t_init: f64, t_final: f64, y_init: f64, h: f64) void = {
let n = ((t_final - t_init) / h): int;
let tn: f64 = t_init;
let yn: f64 = y_init;
let i: int = 1;
fmt::printfln("{: 2} {: 18} {: 21}", "t", "y(t)", "absolute error")!;
fmt::printfln("{: 2} {: 18} {: 21}", tn, yn, math::absf64(exact(tn) - yn))!;
for (i <= n; i += 1) {
yn = rk4(func, tn, yn, h);
tn = t_init + (i: f64)*h;
if (i % 10 == 0) {
fmt::printfln("{: 2} {: 18} {: 21}\t", tn, yn, math::absf64(exact(tn) - yn))!;
};
};
};
fn rk4(func: *fn(_: f64, _: f64) f64, t: f64, y: f64, h: f64) f64 = {
const k1 = func(t, y);
const k2 = func(t + 0.5*h, y + 0.5*h*k1);
const k3 = func(t + 0.5*h, y + 0.5*h*k2);
const k4 = func(t + h, y + h*k3);
return y + h/6.0 * (k1 + 2.0*k2 + 2.0*k3 + k4);
};
fn f(t: f64, y: f64) f64 = {
return t * math::sqrtf64(y);
};
fn exact(t: f64) f64 = {
return 1.0/16.0 * math::powf64(t*t + 4.0, 2.0);
};</syntaxhighlight>
{{out}}
<pre>
t y(t) absolute error
0 1 0
1 1.562499854278108 1.4572189210859676e-7
2 3.9999990805207997 9.194792003341945e-7
3 10.56249709043755 2.909562450525982e-6
4 24.999993765090633 6.23490936746407e-6
5 52.56248918030258 1.0819697422448371e-5
6 99.99998340540358 1.659459641700778e-5
7 175.56247648227125 2.3517728749311573e-5
8 288.9999684347985 3.156520148195341e-5
9 451.5624592768396 4.072316039582802e-5
10 675.9999490167097 5.098329029351589e-5
</pre>
Line 922 ⟶ 1,550:
Using GHC 7.4.1.
<
:: Floating a
=> a -> a -> a
Line 946 ⟶ 1,574:
(print . (\(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)</
Example executed in GHCi:
<
(0,1.0,0.0)
(1,1.5624998542781088,1.4572189122041834e-7)
Line 960 ⟶ 1,588:
(8,288.99996843479926,3.1565204153594095e-5)
(9,451.562459276841,4.0723166534917254e-5)
(10,675.9999490167125,5.098330132113915e-5)</
(See [[Euler method#Haskell]] for implementation of simple general ODE-solver)
Line 966 ⟶ 1,594:
Or, disaggregated a little, and expressed in terms of a single scanl:
<
rk4 y x dx =
let f x y = x * sqrt y
Line 1,010 ⟶ 1,638:
where
justifyLeft n c s = take n (s ++ replicate n c)
justifyRight n c s = drop (length s) (replicate n c ++ s)</
{{Out}}
<pre>y (0) = 1.0 ±0.0
Line 1,026 ⟶ 1,654:
=={{header|J}}==
'''Solution:'''
<
NB. y is: y(ta) , ta , tb , tstep
NB. u is: function to solve
Line 1,043 ⟶ 1,671:
end.
T ,. Y
)</
'''Example:'''
<
fyp=: (* %:)/ NB. f'(t,y)
report_whole=: (10 * i. >:10)&{ NB. report at whole-numbered t values
Line 1,061 ⟶ 1,689:
8 289 _3.15652e_5
9 451.562 _4.07232e_5
10 676 _5.09833e_5</
'''Alternative solution:'''
The following solution replaces the for loop as well as the calculation of the increments (ks) with an accumulating suffix.
<
'Y0 a b h'=. 4{. y
T=. a + i.@>:&.(%&h) b-a
Line 1,082 ⟶ 1,710:
ks=. (x * [: u y + (* x&,))/\. tableau
({:y) + 6 %~ +/ 1 2 2 1 * ks
)</
Use:
Line 1,090 ⟶ 1,718:
Translation of [[Runge-Kutta_method#Ada|Ada]] via [[Runge-Kutta_method#D|D]]
{{works with|Java|8}}
<
import java.util.function.BiFunction;
Line 1,127 ⟶ 1,755:
calc_err(t_arr[i], y_arr[i]));
}
}</
<pre>y(0,0) = 1,00000000 Error: 0,000000
Line 1,143 ⟶ 1,771:
=={{header|JavaScript}}==
===ES5===
<syntaxhighlight lang="javascript">
function rk4(y, x, dx, f) {
var k1 = dx * f(x, y),
Line 1,180 ⟶ 1,808:
steps += 1;
}
</syntaxhighlight>
{{out}}
<pre>
Line 1,197 ⟶ 1,825:
===ES6===
<
'use strict';
Line 1,359 ⟶ 1,987:
// MAIN ---
return main();
})();</
{{Out}}
<pre>y (0) = 1.0 ± 0e+0
Line 1,377 ⟶ 2,005:
They use "while" and/or "until" as defined in recent versions of jq (after version 1.4).
To use either of the two programs with jq 1.4, simply include the lines in the following block:
<
def _until: if cond then . else (next|_until) end;
_until;
Line 1,383 ⟶ 2,011:
def while(cond; update):
def _while: if cond then ., (update | _while) else empty end;
_while;</
===The Example Differential Equation and its Exact Solution===
<
def yprime: .[0] * (.[1] | sqrt);
Line 1,393 ⟶ 2,021:
. as $t
| (( $t*$t) + 4 )
| . * . / 16;</
===dy/dt===
The first solution presented here uses the terminology and style of the
'''Generic filters:'''
<
def round(n):
(if . < 0 then -1 else 1 end) as $s
Line 1,407 ⟶ 2,035:
# Is the input an integer?
def integerq: ((. - ((.+.01) | floor)) | abs) < 0.01;</
'''dy(f)'''
<
# Input: [t, y]; yp is a filter that accepts [t,y] as input
Line 1,423 ⟶ 2,051:
# Input: [t,y]
def dy(f): runge_kutta(f);</
''' Example''':
<
[0,1]
| while( .[0] <= 10;
Line 1,434 ⟶ 2,062:
"y(\($t|round(1))) = \($y|round(10000)) ± \( ($t|actual) - $y | abs)"
else empty
end</
{{out}}
<
y(0) = 1 ± 0
y(1) = 1.5625 ± 1.4572189210859676e-07
Line 1,451 ⟶ 2,079:
real 0m0.048s
user 0m0.013s
sys 0m0.006s</
===newRK4Step===
The second solution follows the nomenclature and style of the Go solution on this page.
Line 1,461 ⟶ 2,090:
The heart of the program is the filter newRK4Step(yp), which is of type ypStepFunc and performs a single
step of the fourth-order Runge-Kutta method, provided yp is of type ypFunc.
<
def newRK4Step(yp):
.[0] as $t | .[1] as $y | .[2] as $dt
Line 1,503 ⟶ 2,132:
# main(t0; y0; tFinal; dtPrint)
main(0; 1; 10; 1)</
{{out}}
<
y(0) = 1 with error: 0
y(1) = 1.562499854278108 with error: 1.4572189210859676e-07
Line 1,520 ⟶ 2,149:
real 0m0.023s
user 0m0.014s
sys 0m0.006s</
=={{header|Julia}}==
{{works with|Julia|0.6}}
=== Using lambda expressions ===
{{trans|Python}}
<
theoric(t) = (t ^ 2 + 4.0) ^ 2 / 16.0
Line 1,549 ⟶ 2,178:
y += δy(t, y, δt)
t += δt
end</
{{out}}
Line 1,564 ⟶ 2,193:
y(10.0) = 675.999949 error: 5.098329e-05</pre>
=== Alternative version ===
{{trans|Python}}
<
vx = Vector{Float64}(undef, n + 1)
vy = Vector{Float64}(undef, n + 1)
vx[1] = x = x₀
vy[1] = y = y₀
Line 1,586 ⟶ 2,215:
for (x, y) in Iterators.take(zip(vx, vy), 10)
@printf("%4.1f %10.5f %+12.4e\n", x, y, y - theoric(x))
end</
=={{header|Kotlin}}==
<
typealias Y = (Double) -> Double
Line 1,623 ⟶ 2,252:
val yd = fun(t: Double, yt: Double) = t * Math.sqrt(yt)
rungeKutta4(0.0, 10.0, 0.1, y, yd)
}</
{{out}}
Line 1,643 ⟶ 2,272:
=={{header|Liberty BASIC}}==
<syntaxhighlight lang="lb">
'[RC] Runge-Kutta method
'initial conditions
Line 1,675 ⟶ 2,304:
exactY=(x^2 + 4)^2 / 16
end function
</syntaxhighlight>
{{Out}}
<pre>
Line 1,691 ⟶ 2,320:
</pre>
=={{header|
<syntaxhighlight lang="lua">local df = function (t, y)
-- derivative of function by value y at time t
return t*y^0.5
end
local dt = 0.1
local y = 1
print ("t", "realY"..' ', "y", ' '.."error")
print ("---", "-------"..' ', "---------------", ' '.."--------------------")
for i = 0, 100 do
local t = i*dt
if t%1 == 0 then
local realY = (t*t+4)^2/16
print (t, realY..' ', y, ' '..realY-y)
end
local dy1 = df(t, y)
local dy2 = df(t+dt/2, y+dt/2*dy1)
local dy3 = df(t+dt/2, y+dt/2*dy2)
local dy4 = df(t+dt, y+dt*dy3)
y = y + dt*(dy1+2*dy2+2*dy3+dy4)/6
end</syntaxhighlight>
{{Out}}
<pre>t realY y error
--- ------- --------------- --------------------
0.0 1.0 1 0.0
1.0 1.5625 1.5624998542781 1.457218921086e-007
2.0 4.0 3.9999990805208 9.1947919989011e-007
3.0 10.5625 10.562497090438 2.9095624469733e-006
4.0 25.0 24.999993765091 6.2349093639114e-006
5.0 52.5625 52.562489180303 1.0819697415343e-005
6.0 100.0 99.999983405404 1.6594596417008e-005
7.0 175.5625 175.56247648227 2.3517728749312e-005
8.0 289.0 288.9999684348 3.156520142511e-005
9.0 451.5625 451.56245927684 4.0723160338985e-005
10.0 676.0 675.99994901671 5.0983290293516e-005
</pre>
=={{header|Mathematica}}/{{header|Wolfram Language}}==
<syntaxhighlight 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}]
Line 1,711 ⟶ 2,381:
solution = NestList[phi, {0, 1}, 101];
Table[{y[[1]], y[[2]], Abs[y[[2]] - 1/16 (y[[1]]^2 + 4)^2]},
{y, solution[[1 ;; 101 ;; 10]]}]
=={{header|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 [http://www.mathworks.com/matlabcentral/answers/98293-is-there-a-fixed-step-ordinary-differential-equation-ode-solver-in-matlab-8-0-r2012b#answer_107643 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.
<
figure
hold on
Line 1,752 ⟶ 2,421:
y(k+1) = y(k)+(dy1+2*dy2+2*dy3+dy4)/6;
end
end</
{{out}}
<pre>
Line 1,770 ⟶ 2,439:
=={{header|Maxima}}==
<
'diff(y, x) = x * sqrt(y);
ode2(%, y, x);
Line 1,809 ⟶ 2,478:
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]);</
=={{header|МК-61/52}}==
Line 1,822 ⟶ 2,491:
''Input:'' 1/2 (h/2) - Р5, 1 (y<sub>0</sub>) - Р8 and Р7, 0 (t<sub>0</sub>) - Р6.
=={{header|Nim}}==
<
proc fn(t, y: float): float =
Line 1,843 ⟶ 2,511:
echo "y(", cur_t, ") = ", cur_y, ", error = ", solution(cur_t) - cur_y
let dy1 = step * fn(cur_t, cur_y)
let dy2 = step * fn(cur_t + 0.5 * step, cur_y + 0.5 * dy1)
let dy3 = step * fn(cur_t + 0.5 * step, cur_y + 0.5 * dy2)
let dy4 = step * fn(cur_t + step, cur_y + dy3)
import math, strformat
proc fn(t, y: float): float =
result = t * math.sqrt(y)
proc solution(t: float): float =
result = (t^2 + 4)^2 / 16
proc rk(start, stop, step: float) =
let nsteps = int(round((stop - start) / step)) + 1
let delta = (stop - start) / float(nsteps - 1)
var cur_y = 1.0
for i in 0..<nsteps:
let cur_t = start + delta * float(i)
if abs(cur_t - math.round(cur_t)) < 1e-5:
echo &"y({cur_t}) = {cur_y}, error = {solution(cur_t) - cur_y}"
let dy1 = step * fn(cur_t, cur_y)
let dy2 = step * fn(cur_t + 0.5 * step, cur_y + 0.5 * dy1)
Line 1,848 ⟶ 2,538:
let dy4 = step * fn(cur_t + step, cur_y + dy3)
cur_y += (dy1 + 2
rk(start =
cur_y += (dy1 + 2.0 * (dy2 + dy3) + dy4) </syntaxhighlight>
{{out}}
<pre>y(0.0) = 1.0, error = 0.0
y(1.0) = 1.562499854278108, error = 1.457218921085968e-
y(2.0) = 3.9999990805208, error = 9.194792003341945e-
y(3.0) = 10.56249709043755, error = 2.909562448749625e-
y(4.0) = 24.99999376509064, error = 6.234909363911356e-
y(5.0) = 52.
y(6.0) = 99.99998340540358, error = 1.659459641700778e-
y(7.0) = 175.5624764822713, error = 2.351772874931157e-
y(8.0) = 288.9999684347986, error = 3.156520142510999e-
y(9.0) = 451.5624592768397, error = 4.07231603389846e-
y(10.0) = 675.9999490167097, error = 5.098329029351589e-
=={{header|Objeck}}==
<
function : Main(args : String[]) ~ Nil {
x0 := 0.0; x1 := 10.0; dx := .1;
Line 1,900 ⟶ 2,591:
return x * y->SquareRoot();
}
}</
Output:
Line 1,918 ⟶ 2,609:
=={{header|OCaml}}==
<
let exact t = let u = 0.25*.t*.t +. 1.0 in u*.u
Line 1,933 ⟶ 2,624:
if n < 102 then loop h (n+1) (rk4_step (y,t) h)
let _ = loop 0.1 1 (1.0, 0.0)</
{{out}}
<pre>t = 0.000000, y = 1.000000, err = 0
Line 1,948 ⟶ 2,639:
=={{header|Octave}}==
<
#Applying the Runge-Kutta method (This code must be implement on a different file than the main one).
Line 1,978 ⟶ 2,669:
fprintf('%d \t %.5f \t %.5f \t %.4g \n',i,f(i),Yn(1+i*10),f(i)-Yn(1+i*10));
endfor
</syntaxhighlight>
{{out}}
<pre>
Line 1,996 ⟶ 2,687:
=={{header|PARI/GP}}==
{{trans|C}}
<
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
Line 2,011 ⟶ 2,702:
)
};
go()</
{{out}}
<pre>x y rel. err.
Line 2,031 ⟶ 2,722:
This code has been compiled using Free Pascal 2.6.2.
<
uses sysutils;
Line 2,101 ⟶ 2,792:
RungeKutta(@YPrime, tArr, yArr, dt);
Print(tArr, yArr, 10);
end.</
{{out}}
<pre>y( 0.0) = 1.00000000 Error: 0.00000E+000
Line 2,123 ⟶ 2,814:
Notice how we have to use sprintf to deal with floating point rounding. See perlfaq4.
<
my ($yp, $dt) = @_;
sub {
Line 2,144 ⟶ 2,835:
printf "y(%2.0f) = %12f ± %e\n", $t, $y, abs($y - ($t**2 + 4)**2 / 16)
if sprintf("%.4f", $t) =~ /0000$/;
}</
{{out}}
<pre>y( 0) = 1.000000 ± 0.000000e+00
Line 2,197 ⟶ 2,852:
=={{header|Phix}}==
{{trans|ERRE}}
<!--<syntaxhighlight lang="phix">(phixonline)-->
<span style="color: #008080;">with</span> <span style="color: #008080;">javascript_semantics</span>
<span style="color: #008080;">constant</span> <span style="color: #000000;">dt</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">0.1</span>
<span style="color: #004080;">atom</span> <span style="color: #000000;">y</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">1.0</span>
<span style="color: #7060A8;">printf</span><span style="color: #0000FF;">(</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span><span style="color: #008000;">" x true/actual y calculated y relative error\n"</span><span style="color: #0000FF;">)</span>
<span style="color: #7060A8;">printf</span><span style="color: #0000FF;">(</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span><span style="color: #008000;">" --- ------------- ------------- --------------\n"</span><span style="color: #0000FF;">)</span>
<span style="color: #008080;">for</span> <span style="color: #000000;">i</span><span style="color: #0000FF;">=</span><span style="color: #000000;">0</span> <span style="color: #008080;">to</span> <span style="color: #000000;">100</span> <span style="color: #008080;">do</span>
<span style="color: #004080;">atom</span> <span style="color: #000000;">t</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">i</span><span style="color: #0000FF;">*</span><span style="color: #000000;">dt</span>
<span style="color: #008080;">if</span> <span style="color: #004080;">integer</span><span style="color: #0000FF;">(</span><span style="color: #000000;">t</span><span style="color: #0000FF;">)</span> <span style="color: #008080;">then</span>
<span style="color: #004080;">atom</span> <span style="color: #000000;">act</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">power</span><span style="color: #0000FF;">(</span><span style="color: #000000;">t</span><span style="color: #0000FF;">*</span><span style="color: #000000;">t</span><span style="color: #0000FF;">+</span><span style="color: #000000;">4</span><span style="color: #0000FF;">,</span><span style="color: #000000;">2</span><span style="color: #0000FF;">)/</span><span style="color: #000000;">16</span>
<span style="color: #7060A8;">printf</span><span style="color: #0000FF;">(</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span><span style="color: #008000;">"%4.1f %14.9f %14.9f %.9e\n"</span><span style="color: #0000FF;">,{</span><span style="color: #000000;">t</span><span style="color: #0000FF;">,</span><span style="color: #000000;">act</span><span style="color: #0000FF;">,</span><span style="color: #000000;">y</span><span style="color: #0000FF;">,</span><span style="color: #7060A8;">abs</span><span style="color: #0000FF;">(</span><span style="color: #000000;">y</span><span style="color: #0000FF;">-</span><span style="color: #000000;">act</span><span style="color: #0000FF;">)})</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">if</span>
<span style="color: #004080;">atom</span> <span style="color: #000000;">k1</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">t</span><span style="color: #0000FF;">*</span><span style="color: #7060A8;">sqrt</span><span style="color: #0000FF;">(</span><span style="color: #000000;">y</span><span style="color: #0000FF;">),</span>
<span style="color: #000000;">k2</span> <span style="color: #0000FF;">=</span> <span style="color: #0000FF;">(</span><span style="color: #000000;">t</span><span style="color: #0000FF;">+</span><span style="color: #000000;">dt</span><span style="color: #0000FF;">/</span><span style="color: #000000;">2</span><span style="color: #0000FF;">)*</span><span style="color: #7060A8;">sqrt</span><span style="color: #0000FF;">(</span><span style="color: #000000;">y</span><span style="color: #0000FF;">+</span><span style="color: #000000;">dt</span><span style="color: #0000FF;">/</span><span style="color: #000000;">2</span><span style="color: #0000FF;">*</span><span style="color: #000000;">k1</span><span style="color: #0000FF;">),</span>
<span style="color: #000000;">k3</span> <span style="color: #0000FF;">=</span> <span style="color: #0000FF;">(</span><span style="color: #000000;">t</span><span style="color: #0000FF;">+</span><span style="color: #000000;">dt</span><span style="color: #0000FF;">/</span><span style="color: #000000;">2</span><span style="color: #0000FF;">)*</span><span style="color: #7060A8;">sqrt</span><span style="color: #0000FF;">(</span><span style="color: #000000;">y</span><span style="color: #0000FF;">+</span><span style="color: #000000;">dt</span><span style="color: #0000FF;">/</span><span style="color: #000000;">2</span><span style="color: #0000FF;">*</span><span style="color: #000000;">k2</span><span style="color: #0000FF;">),</span>
<span style="color: #000000;">k4</span> <span style="color: #0000FF;">=</span> <span style="color: #0000FF;">(</span><span style="color: #000000;">t</span><span style="color: #0000FF;">+</span><span style="color: #000000;">dt</span><span style="color: #0000FF;">)*</span><span style="color: #7060A8;">sqrt</span><span style="color: #0000FF;">(</span><span style="color: #000000;">y</span><span style="color: #0000FF;">+</span><span style="color: #000000;">dt</span><span style="color: #0000FF;">*</span><span style="color: #000000;">k3</span><span style="color: #0000FF;">)</span>
<span style="color: #000000;">y</span> <span style="color: #0000FF;">+=</span> <span style="color: #000000;">dt</span><span style="color: #0000FF;">*(</span><span style="color: #000000;">k1</span><span style="color: #0000FF;">+</span><span style="color: #000000;">2</span><span style="color: #0000FF;">*(</span><span style="color: #000000;">k2</span><span style="color: #0000FF;">+</span><span style="color: #000000;">k3</span><span style="color: #0000FF;">)+</span><span style="color: #000000;">k4</span><span style="color: #0000FF;">)/</span><span style="color: #000000;">6</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">for</span>
<!--</syntaxhighlight>-->
{{out}}
<pre>
Line 2,231 ⟶ 2,889:
=={{header|PL/I}}==
<syntaxhighlight lang="pl/i">
Runge_Kutta: procedure options (main); /* 10 March 2014 */
declare (y, dy1, dy2, dy3, dy4) float (18);
Line 2,257 ⟶ 2,915:
end Runge_kutta;
</syntaxhighlight>
{{out}}
<pre>
Line 2,277 ⟶ 2,935:
{{works with|PowerShell|4.0}}
<syntaxhighlight lang="powershell">
function Runge-Kutta (${function:F}, ${function:y}, $y0, $t0, $dt, $tEnd) {
function RK ($tn,$yn) {
Line 2,314 ⟶ 2,972:
$tEnd = 10
Runge-Kutta F y $y0 $t0 $dt $tEnd
</syntaxhighlight>
<b>Output:</b>
<pre>
Line 2,331 ⟶ 2,989:
10 675.99994901671 5.09832902935159E-05
</pre>
=={{header|PureBasic}}==
{{trans|BBC Basic}}
<
Define.i i
Define.d y=1.0, k1=0.0, k2=0.0, k3=0.0, k4=0.0, t=0.0
Line 2,351 ⟶ 3,010:
Print("Press return to exit...") : Input()
EndIf
End</
{{out}}
<pre>y( 0) = 1.0000 Error = 0.0000000000
Line 2,367 ⟶ 3,026:
=={{header|Python}}==
<syntaxhighlight lang="python">from math import sqrt
def rk4(f, x0, y0, x1, n):
Line 2,440 ⟶ 3,060:
8.0 288.99997 -3.1565e-05
9.0 451.56246 -4.0723e-05
10.0 675.99995 -5.0983e-05</
=={{header|R}}==
<
vx <- double(n + 1)
vy <- double(n + 1)
Line 2,475 ⟶ 3,095:
[9,] 8 288.999968 -3.156520e-05
[10,] 9 451.562459 -4.072316e-05
[11,] 10 675.999949 -5.098329e-05</
=={{header|Racket}}==
Line 2,482 ⟶ 3,102:
The Runge-Kutta method
<
(define (RK4 F δt)
(λ (t y)
Line 2,491 ⟶ 3,111:
(list (+ t δt)
(+ y (* 1/6 (+ δy1 (* 2 δy2) (* 2 δy3) δy4))))))
</syntaxhighlight>
The method modifier which divides each time-step into ''n'' sub-steps:
<
(define ((step-subdivision n method) F h)
(λ (x . y) (last (ODE-solve F (cons x y)
Line 2,500 ⟶ 3,120:
#:step (/ h n)
#:method method))))
</syntaxhighlight>
Usage:
<
(define (F t y) (* t (sqrt y)))
Line 2,514 ⟶ 3,134:
(match-define (list t y) s)
(printf "t=~a\ty=~a\terror=~a\n" t y (- y (exact-solution t))))
</syntaxhighlight>
{{out}}
<pre>
Line 2,532 ⟶ 3,152:
Graphical representation:
<
> (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)")
</syntaxhighlight>
[[File:runge-kutta.png]]
=={{header|Raku}}==
(formerly Perl 6)
{{Works with|rakudo|2016.03}}
<syntaxhighlight lang="raku" line>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, δy($t, $y, δt))
) {
printf "y(%2d) = %12f ± %e\n", $t, $y, abs($y - ($t**2 + 4)**2 / 16)
if $t %% 1;
}</syntaxhighlight>
{{out}}
<pre>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</pre>
=={{header|REXX}}==
<big><big> y'(t) = t<sup>2</sup> √<span style="text-decoration: overline"> y(t) </span></big></big>
The exact solution: <big><big> y(t) = (t<sup>2</sup>+4)<sup>2</sup> ÷ 16 </big></big>
<syntaxhighlight lang="rexx">/*REXX program uses the Runge─Kutta method to solve the equation: y'(t) = t² √[y(t)] */
numeric digits 40; f= digits() % 4 /*use 40 decimal digs, but only show 10*/
x0= 0; x1= 10; dx= .1 /*define variables: X0 X1 DX */
n=1 + (x1-x0) / dx
y.=1; do m=1 for n-1; p= m - 1; y.m= RK4(dx, x0 + dx*p, y.p)
end /*m*/ /* [↑] use 4th order Runge─Kutta. */
w= digits() % 2
say center('X', f, "═") center('Y', w+2, "═") center("relative error", w+8, '═') /*hdr*/
do i=0 to n-1 by 10; x= (x0 + dx*i) / 1;
say center(x, f) fmt(y.i) left('', 2 + ($>=0) ) fmt($)
end /*i*/ /*└┴┴┴───◄─────── aligns positive #'s. */
exit /*stick a fork in it, we're all done. */
/*──────────────────────────────────────────────────────────────────────────────────────*/
fmt: parse arg z; z= right( format(
jus= has. & \hasE; T=
return translate( right(z, (z>=0) + w + 5*hasE + 2*(jus & (z<0) ) ),
/*──────────────────────────────────────────────────────────────────────────────────────*/
RK4: procedure; parse arg dx,x,y; dxH= dx/2; k1= dx * (x ) * sqrt(y )
k2= dx * (x + dxH) * sqrt(y + k1/2)
k3= dx * (x + dxH) * sqrt(y + k2/2)
k4= dx * (x + dx ) * sqrt(y + k3 )
return y + (k1 + k2*2 + k3*2 + k4) / 6
/*──────────────────────────────────────────────────────────────────────────────────────*/
Line 2,574 ⟶ 3,233:
numeric digits; parse value format(x,2,1,,0) 'E0' with g 'E' _ .; g=g * .5'e'_ % 2
do j=0 while h>9; m.j=h; h=h%2+1; end /*j*/
do k=j+5 to 0 by -1; numeric digits m.k; g=(g+x/g)*.5; end /*k*/; return g</
Programming note: the '''fmt''' function is used to
align the output with attention paid to the different ways some
Line 2,612 ⟶ 3,271:
=={{header|Ring}}==
<
decimals(8)
y = 1.0
Line 2,626 ⟶ 3,285:
y += 0.1 * (k1 + 2 * (k2 + k3) + k4) / 6
next
</syntaxhighlight>
Output:
Line 2,644 ⟶ 3,303:
=={{header|Ruby}}==
<
return ->(t,y,dt){
->(dy1 ){
Line 2,670 ⟶ 3,329:
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</
{{Out}}
<pre>
Line 2,687 ⟶ 3,346:
=={{header|Run BASIC}}==
<
while t <= 10
k1 = t * sqr(y)
Line 2,698 ⟶ 3,357:
t = t + .1
wend
end</
{{out}}
<pre>y( 0) = 1.0000000 Error =0
Line 2,715 ⟶ 3,374:
=={{header|Rust}}==
This is a translation of the javascript solution with some minor differences.
<
let k1 = dx * fx(x, y);
let k2 = dx * fx(x + dx / 2.0, y + k1 / 2.0);
Line 2,748 ⟶ 3,407:
x = ((x * 10.0) + (step * 10.0)) / 10.0;
}
}</
<pre>
y(0): 1.0000000000 0E0
Line 2,764 ⟶ 3,423:
=={{header|Scala}}==
<
val f = (t: Double, y: Double) => t * Math.sqrt(y) // Runge-Kutta solution
val g = (t: Double) => Math.pow(t * t + 4, 2) / 16 // Exact solution
Line 2,790 ⟶ 3,449:
}
}
}</
<pre>
y( 0.0) = 1.00000000 Error: 0.00000e+00
Line 2,806 ⟶ 3,465:
=={{header|Sidef}}==
{{trans|
<
func (t, y, δt) {
var a = (δt * yp(t, y));
Line 2,827 ⟶ 3,486:
y += δy(t, y, δt);
t += δt;
}</
{{out}}
<pre>
Line 2,844 ⟶ 3,503:
=={{header|Standard ML}}==
<
let
val dy1 = dt * y'(tn,yn)
Line 2,884 ⟶ 3,543:
(* Run the suggested test case *)
val () = test 0.0 1.0 0.1 101 10 testy testy'</
{{out}}
<pre>Time: 0.0
Line 2,942 ⟶ 3,601:
=={{header|Stata}}==
<
h = (t1-t0)/(n-1)
a = J(n, 2, 0)
Line 2,982 ⟶ 3,641:
10 | 9 451.5624593 -.0000407232 |
11 | 10 675.999949 -.0000509833 |
+----------------------------------------------+</
=={{header|Swift}}==
{{trans|C}}
<
func rk4(dx: Double, x: Double, y: Double, f: (Double, Double) -> Double) -> Double {
Line 3,023 ⟶ 3,682:
print(String(format: "%2g %11.6g %11.5g", x, y[i], y[i]/y2 - 1))
}</
{{out}}
Line 3,041 ⟶ 3,700:
=={{header|Tcl}}==
<
# Hack to bring argument function into expression
Line 3,073 ⟶ 3,732:
printvals $t $y
}
}</
{{out}}
<pre>y(0.0) = 1.00000000 Error: 0.00000000e+00
Line 3,086 ⟶ 3,745:
y(9.0) = 451.56245928 Error: 4.07231581e-05
y(10.0) = 675.99994902 Error: 5.09832864e-05</pre>
=={{header|V (Vlang)}}==
{{trans|Ring}}
<syntaxhighlight lang="Zig">
import math
fn main() {
mut t, mut k1, mut k2, mut k3, mut k4, mut y := 0.0, 0.0, 0.0, 0.0, 0.0, 1.0
for i in 0..101 {
t = i / 10.0
if t == math.floor(t) {
actual := math.pow((math.pow(t, 2) + 4), 2)/16
println("y(${t:.0}) = ${y:.8f} error = ${(actual - y):.8f}")
}
k1 = t * math.sqrt(y)
k2 = (t + 0.05) * math.sqrt(y + 0.05 * k1)
k3 = (t + 0.05) * math.sqrt(y + 0.05 * k2)
k4 = (t + 0.10) * math.sqrt(y + 0.10 * k3)
y += 0.1 * (k1 + 2 * (k2 + k3) + k4) / 6
}
}
</syntaxhighlight>
{{out}}
<pre>
y(0) = 1.00000000 error = 0.00000000
y(1) = 1.56249985 error = 0.00000015
y(2) = 3.99999908 error = 0.00000092
y(3) = 10.56249709 error = 0.00000291
y(4) = 24.99999377 error = 0.00000623
y(5) = 52.56248918 error = 0.00001082
y(6) = 99.99998341 error = 0.00001659
y(7) = 175.56247648 error = 0.00002352
y(8) = 288.99996843 error = 0.00003157
y(9) = 451.56245928 error = 0.00004072
y(10) = 675.99994902 error = 0.00005098
</pre>
=={{header|Wren}}==
{{trans|Kotlin}}
{{libheader|Wren-fmt}}
<syntaxhighlight lang="wren">import "./fmt" for Fmt
var rungeKutta4 = Fn.new { |t0, tz, dt, y, yd|
var tn = t0
var yn = y.call(tn)
var z = ((tz - t0)/dt).truncate
for (i in 0..z) {
if (i % 10 == 0) {
var exact = y.call(tn)
var error = yn - exact
Fmt.print("$4.1f $10f $10f $9f", tn, yn, exact, error)
}
if (i == z) break
var dy1 = dt * yd.call(tn, yn)
var dy2 = dt * yd.call(tn + 0.5 * dt, yn + 0.5 * dy1)
var dy3 = dt * yd.call(tn + 0.5 * dt, yn + 0.5 * dy2)
var dy4 = dt * yd.call(tn + dt, yn + dy3)
yn = yn + (dy1 + 2.0 * dy2 + 2.0 * dy3 + dy4) / 6.0
tn = tn + dt
}
}
System.print(" T RK4 Exact Error")
System.print("---- --------- ---------- ---------")
var y = Fn.new { |t|
var x = t * t + 4.0
return x * x / 16.0
}
var yd = Fn.new { |t, yt| t * yt.sqrt }
rungeKutta4.call(0, 10, 0.1, y, yd)</syntaxhighlight>
{{out}}
<pre>
T RK4 Exact Error
---- --------- ---------- ---------
0.0 1.000000 1.000000 0.000000
1.0 1.562500 1.562500 -0.000000
2.0 3.999999 4.000000 -0.000001
3.0 10.562497 10.562500 -0.000003
4.0 24.999994 25.000000 -0.000006
5.0 52.562489 52.562500 -0.000011
6.0 99.999983 100.000000 -0.000017
7.0 175.562476 175.562500 -0.000024
8.0 288.999968 289.000000 -0.000032
9.0 451.562459 451.562500 -0.000041
10.0 675.999949 676.000000 -0.000051
</pre>
=={{header|XPL0}}==
<syntaxhighlight lang "XPL0">func real Y_(T, Y);
real T, Y;
return T*sqrt(Y);
def DT = 0.1;
real T, Y, Exact, DY1, DY2, DY3, DY4;
[Text(0, " T RK Exact Error^m^j");
T:= 0.; Y:= 1.;
repeat if Mod(T+.001, 1.) < .01 then
[Format(2, 1);
RlOut(0, T);
Format(5, 7);
RlOut(0, Y);
Exact:= sq(T*T+4.)/16.;
RlOut(0, Exact);
RlOut(0, Y-Exact);
CrLf(0);
];
DY1:= DT * Y_(T, Y);
DY2:= DT * Y_(T+DT/2., Y+DY1/2.);
DY3:= DT * Y_(T+DT/2., Y+DY2/2.);
DY4:= DT * Y_(T+DT, Y+DY3);
Y:= Y + (DY1 + 2.*DY2 + 2.*DY3 + DY4) / 6.;
T:= T + DT;
until T > 10.;
]</syntaxhighlight>
{{out}}
<pre>
T RK Exact Error
0.0 1.0000000 1.0000000 0.0000000
1.0 1.5624999 1.5625000 -0.0000001
2.0 3.9999991 4.0000000 -0.0000009
3.0 10.5624971 10.5625000 -0.0000029
4.0 24.9999938 25.0000000 -0.0000062
5.0 52.5624892 52.5625000 -0.0000108
6.0 99.9999834 100.0000000 -0.0000166
7.0 175.5624765 175.5625000 -0.0000235
8.0 288.9999684 289.0000000 -0.0000316
9.0 451.5624593 451.5625000 -0.0000407
10.0 675.9999490 676.0000000 -0.0000510
</pre>
=={{header|zkl}}==
{{trans|OCaml}}
<
fcn exact(t){ u:=0.25*t*t + 1.0; u*u }
Line 3,104 ⟶ 3,894:
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
}</
{{out}}
<pre>
|