Runge-Kutta method: Difference between revisions
Content deleted Content added
Thundergnat (talk | contribs) m syntax highlighting fixup automation |
|||
Line 28: | Line 28: | ||
{{trans|Python}} |
{{trans|Python}} |
||
< |
<syntaxhighlight lang="11l">F rk4(f, x0, y0, x1, n) |
||
V vx = [0.0] * (n + 1) |
V vx = [0.0] * (n + 1) |
||
V vy = [0.0] * (n + 1) |
V vy = [0.0] * (n + 1) |
||
Line 50: | Line 50: | ||
V (vx, vy) = rk4(f, 0.0, 1.0, 10.0, 100) |
V (vx, vy) = rk4(f, 0.0, 1.0, 10.0, 100) |
||
L(x, y) zip(vx, vy)[(0..).step(10)] |
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))</ |
print(‘#2.1 #4.5 #2.8’.format(x, y, y - (4 + x * x) ^ 2 / 16))</syntaxhighlight> |
||
{{out}} |
{{out}} |
||
Line 71: | Line 71: | ||
{{libheader|Action! Tool Kit}} |
{{libheader|Action! Tool Kit}} |
||
{{libheader|Action! Real Math}} |
{{libheader|Action! Real Math}} |
||
< |
<syntaxhighlight lang="action!">INCLUDE "D2:PRINTF.ACT" ;from the Action! Tool Kit |
||
INCLUDE "H6:REALMATH.ACT" |
INCLUDE "H6:REALMATH.ACT" |
||
Line 191: | Line 191: | ||
RealAssign(tmp1,y) ;y=rk4(rate,dx,x0+dx*(i-1),y) |
RealAssign(tmp1,y) ;y=rk4(rate,dx,x0+dx*(i-1),y) |
||
OD |
OD |
||
RETURN</ |
RETURN</syntaxhighlight> |
||
{{out}} |
{{out}} |
||
[https://gitlab.com/amarok8bit/action-rosetta-code/-/raw/master/images/Runge-Kutta_method.png Screenshot from Atari 8-bit computer] |
[https://gitlab.com/amarok8bit/action-rosetta-code/-/raw/master/images/Runge-Kutta_method.png Screenshot from Atari 8-bit computer] |
||
Line 210: | Line 210: | ||
=={{header|Ada}}== |
=={{header|Ada}}== |
||
< |
<syntaxhighlight lang="ada">with Ada.Text_IO; use Ada.Text_IO; |
||
with Ada.Numerics.Generic_Elementary_Functions; |
with Ada.Numerics.Generic_Elementary_Functions; |
||
procedure RungeKutta is |
procedure RungeKutta is |
||
Line 261: | Line 261: | ||
Runge (yprime'Access, t_arr, y_arr, dt); |
Runge (yprime'Access, t_arr, y_arr, dt); |
||
Print (t_arr, y_arr, 10); |
Print (t_arr, y_arr, 10); |
||
end RungeKutta;</ |
end RungeKutta;</syntaxhighlight> |
||
{{out}} |
{{out}} |
||
<pre>y(0.0) = 1.00000000 Error: 0.00000E+00 |
<pre>y(0.0) = 1.00000000 Error: 0.00000E+00 |
||
Line 276: | Line 276: | ||
=={{header|ALGOL 68}}== |
=={{header|ALGOL 68}}== |
||
<syntaxhighlight lang="algol68"> |
|||
<lang ALGOL68> |
|||
BEGIN |
BEGIN |
||
PROC rk4 = (PROC (REAL, REAL) REAL f, REAL y, x, dx) REAL : |
PROC rk4 = (PROC (REAL, REAL) REAL f, REAL y, x, dx) REAL : |
||
Line 303: | Line 303: | ||
OD |
OD |
||
END |
END |
||
</syntaxhighlight> |
|||
</lang> |
|||
{{out}} |
{{out}} |
||
<pre> |
<pre> |
||
Line 323: | Line 323: | ||
{{Trans|ALGOL 68}} |
{{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. |
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 ) ; |
real procedure rk4 ( real procedure f ; real value y, x, dx ) ; |
||
begin % Fourth-order Runge-Kutta method % |
begin % Fourth-order Runge-Kutta method % |
||
Line 352: | Line 352: | ||
end for_i |
end for_i |
||
end |
end |
||
end.</ |
end.</syntaxhighlight> |
||
{{out}} |
{{out}} |
||
<pre> |
<pre> |
||
Line 370: | Line 370: | ||
=={{header|APL}}== |
=={{header|APL}}== |
||
<syntaxhighlight lang="apl"> |
|||
<lang APL> |
|||
∇RK4[⎕]∇ |
∇RK4[⎕]∇ |
||
∇ |
∇ |
||
Line 392: | Line 392: | ||
[2] ⎕←'T' 'RK4 Y' 'ERROR'⍪TABLE,TABLE[;2]-{((4+⍵*2)*2)÷16}TABLE[;1] |
[2] ⎕←'T' 'RK4 Y' 'ERROR'⍪TABLE,TABLE[;2]-{((4+⍵*2)*2)÷16}TABLE[;1] |
||
∇ |
∇ |
||
</syntaxhighlight> |
|||
</lang> |
|||
{{out}} |
{{out}} |
||
<pre> |
<pre> |
||
Line 411: | Line 411: | ||
=={{header|AWK}}== |
=={{header|AWK}}== |
||
<syntaxhighlight lang="awk"> |
|||
<lang AWK> |
|||
# syntax: GAWK -f RUNGE-KUTTA_METHOD.AWK |
# syntax: GAWK -f RUNGE-KUTTA_METHOD.AWK |
||
# converted from BBC BASIC |
# converted from BBC BASIC |
||
Line 431: | Line 431: | ||
exit(0) |
exit(0) |
||
} |
} |
||
</syntaxhighlight> |
|||
</lang> |
|||
{{out}} |
{{out}} |
||
<pre> |
<pre> |
||
Line 450: | Line 450: | ||
=={{header|BASIC}}== |
=={{header|BASIC}}== |
||
==={{header|BASIC256}}=== |
==={{header|BASIC256}}=== |
||
< |
<syntaxhighlight lang="basic256">y = 1 |
||
for i = 0 to 100 |
for i = 0 to 100 |
||
t = i / 10 |
t = i / 10 |
||
Line 465: | Line 465: | ||
y = y + 0.1 * (k1 + 2 * (k2 + k3) + k4) / 6 |
y = y + 0.1 * (k1 + 2 * (k2 + k3) + k4) / 6 |
||
next i |
next i |
||
end</ |
end</syntaxhighlight> |
||
==={{header|BBC BASIC}}=== |
==={{header|BBC BASIC}}=== |
||
< |
<syntaxhighlight lang="bbcbasic"> y = 1.0 |
||
FOR i% = 0 TO 100 |
FOR i% = 0 TO 100 |
||
t = i% / 10 |
t = i% / 10 |
||
Line 483: | Line 483: | ||
k4 = (t + 0.10) * SQR(y + 0.10 * k3) |
k4 = (t + 0.10) * SQR(y + 0.10 * k3) |
||
y += 0.1 * (k1 + 2 * (k2 + k3) + k4) / 6 |
y += 0.1 * (k1 + 2 * (k2 + k3) + k4) / 6 |
||
NEXT i%</ |
NEXT i%</syntaxhighlight> |
||
{{out}} |
{{out}} |
||
<pre>y(0) = 1 Error = 0 |
<pre>y(0) = 1 Error = 0 |
||
Line 499: | Line 499: | ||
==={{header|IS-BASIC}}=== |
==={{header|IS-BASIC}}=== |
||
< |
<syntaxhighlight lang="is-basic">100 PROGRAM "Runge.bas" |
||
110 LET Y=1 |
110 LET Y=1 |
||
120 FOR T=0 TO 10 STEP .1 |
120 FOR T=0 TO 10 STEP .1 |
||
Line 508: | Line 508: | ||
170 LET K4=(T+.1)*SQR(Y+.1*K3) |
170 LET K4=(T+.1)*SQR(Y+.1*K3) |
||
180 LET Y=Y+.1*(K1+2*(K2+K3)+K4)/6 |
180 LET Y=Y+.1*(K1+2*(K2+K3)+K4)/6 |
||
190 NEXT</ |
190 NEXT</syntaxhighlight> |
||
==={{header|QBasic}}=== |
==={{header|QBasic}}=== |
||
{{works with|QBasic|1.1}} |
{{works with|QBasic|1.1}} |
||
{{works with|QuickBasic|4.5}} |
{{works with|QuickBasic|4.5}} |
||
< |
<syntaxhighlight lang="qbasic">y! = 1 |
||
FOR i = 0 TO 100 |
FOR i = 0 TO 100 |
||
t = i / 10 |
t = i / 10 |
||
Line 528: | Line 528: | ||
k4! = (t + .1) * SQR(y + .1 * k3) |
k4! = (t + .1) * SQR(y + .1 * k3) |
||
y = y + .1 * (k1 + 2 * (k2 + k3) + k4) / 6 |
y = y + .1 * (k1 + 2 * (k2 + k3) + k4) / 6 |
||
NEXT i</ |
NEXT i</syntaxhighlight> |
||
==={{header|True BASIC}}=== |
==={{header|True BASIC}}=== |
||
{{works with|QBasic}} |
{{works with|QBasic}} |
||
< |
<syntaxhighlight lang="qbasic">LET y = 1 |
||
FOR i = 0 TO 100 |
FOR i = 0 TO 100 |
||
LET t = i / 10 |
LET t = i / 10 |
||
Line 547: | Line 547: | ||
LET Y = Y + 0.1 * (k1 + 2 * (k2 + k3) + k4) / 6 |
LET Y = Y + 0.1 * (k1 + 2 * (k2 + k3) + k4) / 6 |
||
NEXT i |
NEXT i |
||
END</ |
END</syntaxhighlight> |
||
=={{header|C}}== |
=={{header|C}}== |
||
< |
<syntaxhighlight lang="c">#include <stdio.h> |
||
#include <stdlib.h> |
#include <stdlib.h> |
||
#include <math.h> |
#include <math.h> |
||
Line 587: | Line 587: | ||
return 0; |
return 0; |
||
}</ |
}</syntaxhighlight> |
||
{{out}} (errors are relative) |
{{out}} (errors are relative) |
||
<pre> |
<pre> |
||
Line 606: | Line 606: | ||
=={{header|C sharp|C#}}== |
=={{header|C sharp|C#}}== |
||
< |
<syntaxhighlight lang="csharp"> |
||
using System; |
using System; |
||
Line 700: | Line 700: | ||
} |
} |
||
} |
} |
||
}</ |
}</syntaxhighlight> |
||
=={{header|C++}}== |
=={{header|C++}}== |
||
Using Lambdas |
Using Lambdas |
||
< |
<syntaxhighlight lang="cpp">/* |
||
* compiled with: |
* compiled with: |
||
* g++ (Debian 8.3.0-6) 8.3.0 |
* g++ (Debian 8.3.0-6) 8.3.0 |
||
Line 745: | Line 745: | ||
return 0; |
return 0; |
||
}</ |
}</syntaxhighlight> |
||
=={{header|Common Lisp}}== |
=={{header|Common Lisp}}== |
||
< |
<syntaxhighlight lang="lisp">(defun runge-kutta (f x y x-end n) |
||
(let ((h (float (/ (- x-end x) n) 1d0)) |
(let ((h (float (/ (- x-end x) n) 1d0)) |
||
k1 k2 k3 k4) |
k1 k2 k3 k4) |
||
Line 780: | Line 780: | ||
(7.999999999999988d0 288.9999684347983d0 -3.156520000402452d-5) |
(7.999999999999988d0 288.9999684347983d0 -3.156520000402452d-5) |
||
(8.999999999999984d0 451.56245927683887d0 -4.072315812209126d-5) |
(8.999999999999984d0 451.56245927683887d0 -4.072315812209126d-5) |
||
(9.99999999999998d0 675.9999490167083d0 -5.0983286655537086d-5))</ |
(9.99999999999998d0 675.9999490167083d0 -5.0983286655537086d-5))</syntaxhighlight> |
||
=={{header|Crystal}}== |
=={{header|Crystal}}== |
||
{{trans|Run Basic and Ruby output}} |
{{trans|Run Basic and Ruby output}} |
||
< |
<syntaxhighlight lang="ruby">y, t = 1, 0 |
||
while t <= 10 |
while t <= 10 |
||
k1 = t * Math.sqrt(y) |
k1 = t * Math.sqrt(y) |
||
Line 794: | Line 794: | ||
y += 0.1 * (k1 + 2 * (k2 + k3) + k4) / 6 |
y += 0.1 * (k1 + 2 * (k2 + k3) + k4) / 6 |
||
t += 0.1 |
t += 0.1 |
||
end</ |
end</syntaxhighlight> |
||
{{out}} |
{{out}} |
||
Line 813: | Line 813: | ||
=={{header|D}}== |
=={{header|D}}== |
||
{{trans|Ada}} |
{{trans|Ada}} |
||
< |
<syntaxhighlight lang="d">import std.stdio, std.math, std.typecons; |
||
alias FP = real; |
alias FP = real; |
||
Line 850: | Line 850: | ||
t_arr[i], y_arr[i], |
t_arr[i], y_arr[i], |
||
calc_err(t_arr[i], y_arr[i])); |
calc_err(t_arr[i], y_arr[i])); |
||
}</ |
}</syntaxhighlight> |
||
{{out}} |
{{out}} |
||
<pre>y(0.0) = 1.00000000 Error: 0 |
<pre>y(0.0) = 1.00000000 Error: 0 |
||
Line 865: | Line 865: | ||
=={{header|Dart}}== |
=={{header|Dart}}== |
||
< |
<syntaxhighlight lang="dart">import 'dart:math' as Math; |
||
num RungeKutta4(Function f, num t, num y, num dt){ |
num RungeKutta4(Function f, num t, num y, num dt){ |
||
Line 891: | Line 891: | ||
t += dt; |
t += dt; |
||
} |
} |
||
}</ |
}</syntaxhighlight> |
||
{{out}} |
{{out}} |
||
<pre> |
<pre> |
||
Line 913: | Line 913: | ||
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. |
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. |
[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. |
Full description is in Wilkes, Wheeler & Gill, 1951 edn, pages 32-34, 86-87, 132-134. |
||
Line 1,058: | Line 1,058: | ||
E25K TA GK |
E25K TA GK |
||
E10Z PF [enter at relative address 10 with accumulator = 0] |
E10Z PF [enter at relative address 10 with accumulator = 0] |
||
</syntaxhighlight> |
|||
</lang> |
|||
{{out}} |
{{out}} |
||
<pre> |
<pre> |
||
Line 1,077: | Line 1,077: | ||
=={{header|ERRE}}== |
=={{header|ERRE}}== |
||
<syntaxhighlight lang="erre"> |
|||
<lang ERRE> |
|||
PROGRAM RUNGE_KUTTA |
PROGRAM RUNGE_KUTTA |
||
Line 1,102: | Line 1,102: | ||
Y+=DELTA_T*(K1+2*(K2+K3)+K4)/6 |
Y+=DELTA_T*(K1+2*(K2+K3)+K4)/6 |
||
END FOR |
END FOR |
||
END PROGRAM</ |
END PROGRAM</syntaxhighlight> |
||
{{out}} |
{{out}} |
||
<pre> |
<pre> |
||
Line 1,120: | Line 1,120: | ||
=={{header|F_Sharp|F#}}== |
=={{header|F_Sharp|F#}}== |
||
{{works with|F# interactive (fsi.exe)}} |
{{works with|F# interactive (fsi.exe)}} |
||
< |
<syntaxhighlight lang="fsharp"> |
||
open System |
open System |
||
Line 1,141: | Line 1,141: | ||
RungeKutta4 0.0 1.0 10.0 0.1 |
RungeKutta4 0.0 1.0 10.0 0.1 |
||
|> Seq.filter (fun (t,y) -> t % 1.0 = 0.0 ) |
|> 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) )</ |
|> Seq.iter (fun (t,y) -> Console.WriteLine("y({0})={1}\t(relative error:{2})", t, y, (y / y_exact(t))-1.0) )</syntaxhighlight> |
||
{{out}} |
{{out}} |
||
Line 1,159: | Line 1,159: | ||
=={{header|Fortran}}== |
=={{header|Fortran}}== |
||
< |
<syntaxhighlight lang="fortran">program rungekutta |
||
implicit none |
implicit none |
||
integer, parameter :: dp = kind(1d0) |
integer, parameter :: dp = kind(1d0) |
||
Line 1,191: | Line 1,191: | ||
f = t*sqrt(y) |
f = t*sqrt(y) |
||
end function f |
end function f |
||
end program rungekutta</ |
end program rungekutta</syntaxhighlight> |
||
{{out}} |
{{out}} |
||
<pre> |
<pre> |
||
Line 1,209: | Line 1,209: | ||
=={{header|FreeBASIC}}== |
=={{header|FreeBASIC}}== |
||
{{trans|BBC BASIC}} |
{{trans|BBC BASIC}} |
||
< |
<syntaxhighlight lang="freebasic">' version 03-10-2015 |
||
' compile with: fbc -s console |
' compile with: fbc -s console |
||
' translation of BBC BASIC |
' translation of BBC BASIC |
||
Line 1,239: | Line 1,239: | ||
Print : Print "hit any key to end program" |
Print : Print "hit any key to end program" |
||
Sleep |
Sleep |
||
End</ |
End</syntaxhighlight> |
||
{{out}} |
{{out}} |
||
<pre>y(0) = 1 Error = 0 |
<pre>y(0) = 1 Error = 0 |
||
Line 1,254: | Line 1,254: | ||
=={{header|FutureBasic}}== |
=={{header|FutureBasic}}== |
||
< |
<syntaxhighlight lang="futurebasic">window 1 |
||
def fn dydx( x as double, y as double ) as double = x * sqr(y) |
def fn dydx( x as double, y as double ) as double = x * sqr(y) |
||
Line 1,279: | Line 1,279: | ||
next |
next |
||
HandleEvents</ |
HandleEvents</syntaxhighlight> |
||
Output: |
Output: |
||
<pre> |
<pre> |
||
Line 1,297: | Line 1,297: | ||
=={{header|Go}}== |
=={{header|Go}}== |
||
{{works with|Go1}} |
{{works with|Go1}} |
||
< |
<syntaxhighlight lang="go">package main |
||
import ( |
import ( |
||
Line 1,353: | Line 1,353: | ||
func printErr(t, y float64) { |
func printErr(t, y float64) { |
||
fmt.Printf("y(%.1f) = %f Error: %e\n", t, y, math.Abs(actual(t)-y)) |
fmt.Printf("y(%.1f) = %f Error: %e\n", t, y, math.Abs(actual(t)-y)) |
||
}</ |
}</syntaxhighlight> |
||
{{out}} |
{{out}} |
||
<pre> |
<pre> |
||
Line 1,370: | Line 1,370: | ||
=={{header|Groovy}}== |
=={{header|Groovy}}== |
||
<syntaxhighlight lang="groovy"> |
|||
<lang Groovy> |
|||
class Runge_Kutta{ |
class Runge_Kutta{ |
||
static void main(String[] args){ |
static void main(String[] args){ |
||
Line 1,395: | Line 1,395: | ||
static def dy(def x){return x*0.1;} |
static def dy(def x){return x*0.1;} |
||
} |
} |
||
</syntaxhighlight> |
|||
</lang> |
|||
{{out}} |
{{out}} |
||
<pre> |
<pre> |
||
Line 1,412: | Line 1,412: | ||
=={{header|Hare}}== |
=={{header|Hare}}== |
||
< |
<syntaxhighlight lang="hare">use fmt; |
||
use math; |
use math; |
||
Line 1,452: | Line 1,452: | ||
fn exact(t: f64) f64 = { |
fn exact(t: f64) f64 = { |
||
return 1.0/16.0 * math::powf64(t*t + 4.0, 2.0); |
return 1.0/16.0 * math::powf64(t*t + 4.0, 2.0); |
||
};</ |
};</syntaxhighlight> |
||
{{out}} |
{{out}} |
||
<pre> |
<pre> |
||
Line 1,473: | Line 1,473: | ||
Using GHC 7.4.1. |
Using GHC 7.4.1. |
||
< |
<syntaxhighlight lang="haskell">dv |
||
:: Floating a |
:: Floating a |
||
=> a -> a -> a |
=> a -> a -> a |
||
Line 1,497: | Line 1,497: | ||
(print . (\(x, y) -> (truncate x, y, fy x - y))) |
(print . (\(x, y) -> (truncate x, y, fy x - y))) |
||
(filter (\(x, _) -> 0 == mod (truncate $ 10 * x) 10) $ |
(filter (\(x, _) -> 0 == mod (truncate $ 10 * x) 10) $ |
||
take 101 $ rk4 dv 1.0 0 0.1)</ |
take 101 $ rk4 dv 1.0 0 0.1)</syntaxhighlight> |
||
Example executed in GHCi: |
Example executed in GHCi: |
||
< |
<syntaxhighlight lang="haskell">*Main> task |
||
(0,1.0,0.0) |
(0,1.0,0.0) |
||
(1,1.5624998542781088,1.4572189122041834e-7) |
(1,1.5624998542781088,1.4572189122041834e-7) |
||
Line 1,511: | Line 1,511: | ||
(8,288.99996843479926,3.1565204153594095e-5) |
(8,288.99996843479926,3.1565204153594095e-5) |
||
(9,451.562459276841,4.0723166534917254e-5) |
(9,451.562459276841,4.0723166534917254e-5) |
||
(10,675.9999490167125,5.098330132113915e-5)</ |
(10,675.9999490167125,5.098330132113915e-5)</syntaxhighlight> |
||
(See [[Euler method#Haskell]] for implementation of simple general ODE-solver) |
(See [[Euler method#Haskell]] for implementation of simple general ODE-solver) |
||
Line 1,517: | Line 1,517: | ||
Or, disaggregated a little, and expressed in terms of a single scanl: |
Or, disaggregated a little, and expressed in terms of a single scanl: |
||
< |
<syntaxhighlight lang="haskell">rk4 :: Double -> Double -> Double -> Double |
||
rk4 y x dx = |
rk4 y x dx = |
||
let f x y = x * sqrt y |
let f x y = x * sqrt y |
||
Line 1,561: | Line 1,561: | ||
where |
where |
||
justifyLeft n c s = take n (s ++ replicate n c) |
justifyLeft n c s = take n (s ++ replicate n c) |
||
justifyRight n c s = drop (length s) (replicate n c ++ s)</ |
justifyRight n c s = drop (length s) (replicate n c ++ s)</syntaxhighlight> |
||
{{Out}} |
{{Out}} |
||
<pre>y (0) = 1.0 ±0.0 |
<pre>y (0) = 1.0 ±0.0 |
||
Line 1,577: | Line 1,577: | ||
=={{header|J}}== |
=={{header|J}}== |
||
'''Solution:''' |
'''Solution:''' |
||
< |
<syntaxhighlight lang="j">NB.*rk4 a Solve function using Runge-Kutta method |
||
NB. y is: y(ta) , ta , tb , tstep |
NB. y is: y(ta) , ta , tb , tstep |
||
NB. u is: function to solve |
NB. u is: function to solve |
||
Line 1,594: | Line 1,594: | ||
end. |
end. |
||
T ,. Y |
T ,. Y |
||
)</ |
)</syntaxhighlight> |
||
'''Example:''' |
'''Example:''' |
||
< |
<syntaxhighlight lang="j"> fy=: (%16) * [: *: 4 + *: NB. f(t,y) |
||
fyp=: (* %:)/ NB. f'(t,y) |
fyp=: (* %:)/ NB. f'(t,y) |
||
report_whole=: (10 * i. >:10)&{ NB. report at whole-numbered t values |
report_whole=: (10 * i. >:10)&{ NB. report at whole-numbered t values |
||
Line 1,612: | Line 1,612: | ||
8 289 _3.15652e_5 |
8 289 _3.15652e_5 |
||
9 451.562 _4.07232e_5 |
9 451.562 _4.07232e_5 |
||
10 676 _5.09833e_5</ |
10 676 _5.09833e_5</syntaxhighlight> |
||
'''Alternative solution:''' |
'''Alternative solution:''' |
||
The following solution replaces the for loop as well as the calculation of the increments (ks) with an accumulating suffix. |
The following solution replaces the for loop as well as the calculation of the increments (ks) with an accumulating suffix. |
||
< |
<syntaxhighlight lang="j">rk4=: adverb define |
||
'Y0 a b h'=. 4{. y |
'Y0 a b h'=. 4{. y |
||
T=. a + i.@>:&.(%&h) b-a |
T=. a + i.@>:&.(%&h) b-a |
||
Line 1,633: | Line 1,633: | ||
ks=. (x * [: u y + (* x&,))/\. tableau |
ks=. (x * [: u y + (* x&,))/\. tableau |
||
({:y) + 6 %~ +/ 1 2 2 1 * ks |
({:y) + 6 %~ +/ 1 2 2 1 * ks |
||
)</ |
)</syntaxhighlight> |
||
Use: |
Use: |
||
Line 1,641: | Line 1,641: | ||
Translation of [[Runge-Kutta_method#Ada|Ada]] via [[Runge-Kutta_method#D|D]] |
Translation of [[Runge-Kutta_method#Ada|Ada]] via [[Runge-Kutta_method#D|D]] |
||
{{works with|Java|8}} |
{{works with|Java|8}} |
||
< |
<syntaxhighlight lang="java">import static java.lang.Math.*; |
||
import java.util.function.BiFunction; |
import java.util.function.BiFunction; |
||
Line 1,678: | Line 1,678: | ||
calc_err(t_arr[i], y_arr[i])); |
calc_err(t_arr[i], y_arr[i])); |
||
} |
} |
||
}</ |
}</syntaxhighlight> |
||
<pre>y(0,0) = 1,00000000 Error: 0,000000 |
<pre>y(0,0) = 1,00000000 Error: 0,000000 |
||
Line 1,694: | Line 1,694: | ||
=={{header|JavaScript}}== |
=={{header|JavaScript}}== |
||
===ES5=== |
===ES5=== |
||
<syntaxhighlight lang="javascript"> |
|||
<lang JavaScript> |
|||
function rk4(y, x, dx, f) { |
function rk4(y, x, dx, f) { |
||
var k1 = dx * f(x, y), |
var k1 = dx * f(x, y), |
||
Line 1,731: | Line 1,731: | ||
steps += 1; |
steps += 1; |
||
} |
} |
||
</syntaxhighlight> |
|||
</lang> |
|||
{{out}} |
{{out}} |
||
<pre> |
<pre> |
||
Line 1,748: | Line 1,748: | ||
===ES6=== |
===ES6=== |
||
< |
<syntaxhighlight lang="javascript">(() => { |
||
'use strict'; |
'use strict'; |
||
Line 1,910: | Line 1,910: | ||
// MAIN --- |
// MAIN --- |
||
return main(); |
return main(); |
||
})();</ |
})();</syntaxhighlight> |
||
{{Out}} |
{{Out}} |
||
<pre>y (0) = 1.0 ± 0e+0 |
<pre>y (0) = 1.0 ± 0e+0 |
||
Line 1,928: | Line 1,928: | ||
They use "while" and/or "until" as defined in recent versions of jq (after version 1.4). |
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: |
To use either of the two programs with jq 1.4, simply include the lines in the following block: |
||
< |
<syntaxhighlight lang="jq">def until(cond; next): |
||
def _until: if cond then . else (next|_until) end; |
def _until: if cond then . else (next|_until) end; |
||
_until; |
_until; |
||
Line 1,934: | Line 1,934: | ||
def while(cond; update): |
def while(cond; update): |
||
def _while: if cond then ., (update | _while) else empty end; |
def _while: if cond then ., (update | _while) else empty end; |
||
_while;</ |
_while;</syntaxhighlight> |
||
===The Example Differential Equation and its Exact Solution=== |
===The Example Differential Equation and its Exact Solution=== |
||
< |
<syntaxhighlight lang="jq"># yprime maps [t,y] to a number, i.e. t * sqrt(y) |
||
def yprime: .[0] * (.[1] | sqrt); |
def yprime: .[0] * (.[1] | sqrt); |
||
Line 1,944: | Line 1,944: | ||
. as $t |
. as $t |
||
| (( $t*$t) + 4 ) |
| (( $t*$t) + 4 ) |
||
| . * . / 16;</ |
| . * . / 16;</syntaxhighlight> |
||
===dy/dt=== |
===dy/dt=== |
||
Line 1,950: | Line 1,950: | ||
'''Generic filters:''' |
'''Generic filters:''' |
||
< |
<syntaxhighlight lang="jq"># n is the number of decimal places of precision |
||
def round(n): |
def round(n): |
||
(if . < 0 then -1 else 1 end) as $s |
(if . < 0 then -1 else 1 end) as $s |
||
Line 1,958: | Line 1,958: | ||
# Is the input an integer? |
# Is the input an integer? |
||
def integerq: ((. - ((.+.01) | floor)) | abs) < 0.01;</ |
def integerq: ((. - ((.+.01) | floor)) | abs) < 0.01;</syntaxhighlight> |
||
'''dy(f)''' |
'''dy(f)''' |
||
< |
<syntaxhighlight lang="jq">def dt: 0.1; |
||
# Input: [t, y]; yp is a filter that accepts [t,y] as input |
# Input: [t, y]; yp is a filter that accepts [t,y] as input |
||
Line 1,974: | Line 1,974: | ||
# Input: [t,y] |
# Input: [t,y] |
||
def dy(f): runge_kutta(f);</ |
def dy(f): runge_kutta(f);</syntaxhighlight> |
||
''' Example''': |
''' Example''': |
||
< |
<syntaxhighlight lang="jq"># state: [t,y] |
||
[0,1] |
[0,1] |
||
| while( .[0] <= 10; |
| while( .[0] <= 10; |
||
Line 1,985: | Line 1,985: | ||
"y(\($t|round(1))) = \($y|round(10000)) ± \( ($t|actual) - $y | abs)" |
"y(\($t|round(1))) = \($y|round(10000)) ± \( ($t|actual) - $y | abs)" |
||
else empty |
else empty |
||
end</ |
end</syntaxhighlight> |
||
{{out}} |
{{out}} |
||
< |
<syntaxhighlight lang="sh">$ time jq -r -n -f rk4.pl.jq |
||
y(0) = 1 ± 0 |
y(0) = 1 ± 0 |
||
y(1) = 1.5625 ± 1.4572189210859676e-07 |
y(1) = 1.5625 ± 1.4572189210859676e-07 |
||
Line 2,002: | Line 2,002: | ||
real 0m0.048s |
real 0m0.048s |
||
user 0m0.013s |
user 0m0.013s |
||
sys 0m0.006s</ |
sys 0m0.006s</syntaxhighlight> |
||
===newRK4Step=== |
===newRK4Step=== |
||
Line 2,013: | Line 2,013: | ||
The heart of the program is the filter newRK4Step(yp), which is of type ypStepFunc and performs a single |
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. |
step of the fourth-order Runge-Kutta method, provided yp is of type ypFunc. |
||
< |
<syntaxhighlight lang="jq"># Input: [t, y, dt] |
||
def newRK4Step(yp): |
def newRK4Step(yp): |
||
.[0] as $t | .[1] as $y | .[2] as $dt |
.[0] as $t | .[1] as $y | .[2] as $dt |
||
Line 2,055: | Line 2,055: | ||
# main(t0; y0; tFinal; dtPrint) |
# main(t0; y0; tFinal; dtPrint) |
||
main(0; 1; 10; 1)</ |
main(0; 1; 10; 1)</syntaxhighlight> |
||
{{out}} |
{{out}} |
||
< |
<syntaxhighlight lang="sh">$ time jq -n -r -f runge-kutta.jq |
||
y(0) = 1 with error: 0 |
y(0) = 1 with error: 0 |
||
y(1) = 1.562499854278108 with error: 1.4572189210859676e-07 |
y(1) = 1.562499854278108 with error: 1.4572189210859676e-07 |
||
Line 2,072: | Line 2,072: | ||
real 0m0.023s |
real 0m0.023s |
||
user 0m0.014s |
user 0m0.014s |
||
sys 0m0.006s</ |
sys 0m0.006s</syntaxhighlight> |
||
=={{header|Julia}}== |
=={{header|Julia}}== |
||
Line 2,079: | Line 2,079: | ||
=== Using lambda expressions === |
=== Using lambda expressions === |
||
{{trans|Python}} |
{{trans|Python}} |
||
< |
<syntaxhighlight lang="julia">f(x, y) = x * sqrt(y) |
||
theoric(t) = (t ^ 2 + 4.0) ^ 2 / 16.0 |
theoric(t) = (t ^ 2 + 4.0) ^ 2 / 16.0 |
||
Line 2,101: | Line 2,101: | ||
y += δy(t, y, δt) |
y += δy(t, y, δt) |
||
t += δt |
t += δt |
||
end</ |
end</syntaxhighlight> |
||
{{out}} |
{{out}} |
||
Line 2,118: | Line 2,118: | ||
=== Alternative version === |
=== Alternative version === |
||
{{trans|Python}} |
{{trans|Python}} |
||
< |
<syntaxhighlight lang="julia">function rk4(f::Function, x₀::Float64, y₀::Float64, x₁::Float64, n) |
||
vx = Vector{Float64}(undef, n + 1) |
vx = Vector{Float64}(undef, n + 1) |
||
vy = Vector{Float64}(undef, n + 1) |
vy = Vector{Float64}(undef, n + 1) |
||
Line 2,138: | Line 2,138: | ||
for (x, y) in Iterators.take(zip(vx, vy), 10) |
for (x, y) in Iterators.take(zip(vx, vy), 10) |
||
@printf("%4.1f %10.5f %+12.4e\n", x, y, y - theoric(x)) |
@printf("%4.1f %10.5f %+12.4e\n", x, y, y - theoric(x)) |
||
end</ |
end</syntaxhighlight> |
||
=={{header|Kotlin}}== |
=={{header|Kotlin}}== |
||
< |
<syntaxhighlight lang="scala">// version 1.1.2 |
||
typealias Y = (Double) -> Double |
typealias Y = (Double) -> Double |
||
Line 2,175: | Line 2,175: | ||
val yd = fun(t: Double, yt: Double) = t * Math.sqrt(yt) |
val yd = fun(t: Double, yt: Double) = t * Math.sqrt(yt) |
||
rungeKutta4(0.0, 10.0, 0.1, y, yd) |
rungeKutta4(0.0, 10.0, 0.1, y, yd) |
||
}</ |
}</syntaxhighlight> |
||
{{out}} |
{{out}} |
||
Line 2,195: | Line 2,195: | ||
=={{header|Liberty BASIC}}== |
=={{header|Liberty BASIC}}== |
||
<syntaxhighlight lang="lb"> |
|||
<lang lb> |
|||
'[RC] Runge-Kutta method |
'[RC] Runge-Kutta method |
||
'initial conditions |
'initial conditions |
||
Line 2,227: | Line 2,227: | ||
exactY=(x^2 + 4)^2 / 16 |
exactY=(x^2 + 4)^2 / 16 |
||
end function |
end function |
||
</syntaxhighlight> |
|||
</lang> |
|||
{{Out}} |
{{Out}} |
||
<pre> |
<pre> |
||
Line 2,244: | Line 2,244: | ||
=={{header|Lua}}== |
=={{header|Lua}}== |
||
< |
<syntaxhighlight lang="lua">local df = function (t, y) |
||
-- derivative of function by value y at time t |
-- derivative of function by value y at time t |
||
return t*y^0.5 |
return t*y^0.5 |
||
Line 2,266: | Line 2,266: | ||
local dy4 = df(t+dt, y+dt*dy3) |
local dy4 = df(t+dt, y+dt*dy3) |
||
y = y + dt*(dy1+2*dy2+2*dy3+dy4)/6 |
y = y + dt*(dy1+2*dy2+2*dy3+dy4)/6 |
||
end</ |
end</syntaxhighlight> |
||
{{Out}} |
{{Out}} |
||
<pre>t realY y error |
<pre>t realY y error |
||
Line 2,285: | Line 2,285: | ||
=={{header|Mathematica}}/{{header|Wolfram Language}}== |
=={{header|Mathematica}}/{{header|Wolfram Language}}== |
||
< |
<syntaxhighlight lang="mathematica">(* Symbolic solution *) |
||
DSolve[{y'[t] == t*Sqrt[y[t]], y[0] == 1}, y, t] |
DSolve[{y'[t] == t*Sqrt[y[t]], y[0] == 1}, y, t] |
||
Table[{t, 1/16 (4 + t^2)^2}, {t, 0, 10}] |
Table[{t, 1/16 (4 + t^2)^2}, {t, 0, 10}] |
||
Line 2,304: | Line 2,304: | ||
solution = NestList[phi, {0, 1}, 101]; |
solution = NestList[phi, {0, 1}, 101]; |
||
Table[{y[[1]], y[[2]], Abs[y[[2]] - 1/16 (y[[1]]^2 + 4)^2]}, |
Table[{y[[1]], y[[2]], Abs[y[[2]] - 1/16 (y[[1]]^2 + 4)^2]}, |
||
{y, solution[[1 ;; 101 ;; 10]]}]</ |
{y, solution[[1 ;; 101 ;; 10]]}]</syntaxhighlight> |
||
=={{header|MATLAB}}== |
=={{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. |
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. |
||
< |
<syntaxhighlight lang="matlab">function testRK4Programs |
||
figure |
figure |
||
hold on |
hold on |
||
Line 2,344: | Line 2,344: | ||
y(k+1) = y(k)+(dy1+2*dy2+2*dy3+dy4)/6; |
y(k+1) = y(k)+(dy1+2*dy2+2*dy3+dy4)/6; |
||
end |
end |
||
end</ |
end</syntaxhighlight> |
||
{{out}} |
{{out}} |
||
<pre> |
<pre> |
||
Line 2,362: | Line 2,362: | ||
=={{header|Maxima}}== |
=={{header|Maxima}}== |
||
< |
<syntaxhighlight lang="maxima">/* Here is how to solve a differential equation */ |
||
'diff(y, x) = x * sqrt(y); |
'diff(y, x) = x * sqrt(y); |
||
ode2(%, y, x); |
ode2(%, y, x); |
||
Line 2,401: | Line 2,401: | ||
s: map(lambda([x], (x^2 + 4)^2 / 16), x)$ |
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]);</ |
for i from 1 step 10 thru 101 do print(x[i], " ", y[i], " ", y[i] - s[i]);</syntaxhighlight> |
||
=={{header|МК-61/52}}== |
=={{header|МК-61/52}}== |
||
Line 2,416: | Line 2,416: | ||
=={{header|Nim}}== |
=={{header|Nim}}== |
||
< |
<syntaxhighlight lang="nim">import math |
||
proc fn(t, y: float): float = |
proc fn(t, y: float): float = |
||
Line 2,464: | Line 2,464: | ||
rk(start = 0, stop = 10, step = 0.1) |
rk(start = 0, stop = 10, step = 0.1) |
||
cur_y += (dy1 + 2.0 * (dy2 + dy3) + dy4) </ |
cur_y += (dy1 + 2.0 * (dy2 + dy3) + dy4) </syntaxhighlight> |
||
{{out}} |
{{out}} |
||
<pre>y(0.0) = 1.0, error = 0.0 |
<pre>y(0.0) = 1.0, error = 0.0 |
||
Line 2,479: | Line 2,479: | ||
=={{header|Objeck}}== |
=={{header|Objeck}}== |
||
< |
<syntaxhighlight lang="objeck">class RungeKuttaMethod { |
||
function : Main(args : String[]) ~ Nil { |
function : Main(args : String[]) ~ Nil { |
||
x0 := 0.0; x1 := 10.0; dx := .1; |
x0 := 0.0; x1 := 10.0; dx := .1; |
||
Line 2,514: | Line 2,514: | ||
return x * y->SquareRoot(); |
return x * y->SquareRoot(); |
||
} |
} |
||
}</ |
}</syntaxhighlight> |
||
Output: |
Output: |
||
Line 2,532: | Line 2,532: | ||
=={{header|OCaml}}== |
=={{header|OCaml}}== |
||
< |
<syntaxhighlight lang="ocaml">let y' t y = t *. sqrt y |
||
let exact t = let u = 0.25*.t*.t +. 1.0 in u*.u |
let exact t = let u = 0.25*.t*.t +. 1.0 in u*.u |
||
Line 2,547: | Line 2,547: | ||
if n < 102 then loop h (n+1) (rk4_step (y,t) h) |
if n < 102 then loop h (n+1) (rk4_step (y,t) h) |
||
let _ = loop 0.1 1 (1.0, 0.0)</ |
let _ = loop 0.1 1 (1.0, 0.0)</syntaxhighlight> |
||
{{out}} |
{{out}} |
||
<pre>t = 0.000000, y = 1.000000, err = 0 |
<pre>t = 0.000000, y = 1.000000, err = 0 |
||
Line 2,562: | Line 2,562: | ||
=={{header|Octave}}== |
=={{header|Octave}}== |
||
< |
<syntaxhighlight lang="octave"> |
||
#Applying the Runge-Kutta method (This code must be implement on a different file than the main one). |
#Applying the Runge-Kutta method (This code must be implement on a different file than the main one). |
||
Line 2,592: | Line 2,592: | ||
fprintf('%d \t %.5f \t %.5f \t %.4g \n',i,f(i),Yn(1+i*10),f(i)-Yn(1+i*10)); |
fprintf('%d \t %.5f \t %.5f \t %.4g \n',i,f(i),Yn(1+i*10),f(i)-Yn(1+i*10)); |
||
endfor |
endfor |
||
</syntaxhighlight> |
|||
</lang> |
|||
{{out}} |
{{out}} |
||
<pre> |
<pre> |
||
Line 2,610: | Line 2,610: | ||
=={{header|PARI/GP}}== |
=={{header|PARI/GP}}== |
||
{{trans|C}} |
{{trans|C}} |
||
< |
<syntaxhighlight 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)); |
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 |
y + (k1 + 2*k2 + 2*k3 + k4) / 6 |
||
Line 2,625: | Line 2,625: | ||
) |
) |
||
}; |
}; |
||
go()</ |
go()</syntaxhighlight> |
||
{{out}} |
{{out}} |
||
<pre>x y rel. err. |
<pre>x y rel. err. |
||
Line 2,645: | Line 2,645: | ||
This code has been compiled using Free Pascal 2.6.2. |
This code has been compiled using Free Pascal 2.6.2. |
||
< |
<syntaxhighlight lang="pascal">program RungeKuttaExample; |
||
uses sysutils; |
uses sysutils; |
||
Line 2,715: | Line 2,715: | ||
RungeKutta(@YPrime, tArr, yArr, dt); |
RungeKutta(@YPrime, tArr, yArr, dt); |
||
Print(tArr, yArr, 10); |
Print(tArr, yArr, 10); |
||
end.</ |
end.</syntaxhighlight> |
||
{{out}} |
{{out}} |
||
<pre>y( 0.0) = 1.00000000 Error: 0.00000E+000 |
<pre>y( 0.0) = 1.00000000 Error: 0.00000E+000 |
||
Line 2,737: | Line 2,737: | ||
Notice how we have to use sprintf to deal with floating point rounding. See perlfaq4. |
Notice how we have to use sprintf to deal with floating point rounding. See perlfaq4. |
||
< |
<syntaxhighlight lang="perl">sub runge_kutta { |
||
my ($yp, $dt) = @_; |
my ($yp, $dt) = @_; |
||
sub { |
sub { |
||
Line 2,758: | Line 2,758: | ||
printf "y(%2.0f) = %12f ± %e\n", $t, $y, abs($y - ($t**2 + 4)**2 / 16) |
printf "y(%2.0f) = %12f ± %e\n", $t, $y, abs($y - ($t**2 + 4)**2 / 16) |
||
if sprintf("%.4f", $t) =~ /0000$/; |
if sprintf("%.4f", $t) =~ /0000$/; |
||
}</ |
}</syntaxhighlight> |
||
{{out}} |
{{out}} |
||
Line 2,775: | Line 2,775: | ||
=={{header|Phix}}== |
=={{header|Phix}}== |
||
{{trans|ERRE}} |
{{trans|ERRE}} |
||
<!--< |
<!--<syntaxhighlight lang="phix">(phixonline)--> |
||
<span style="color: #008080;">with</span> <span style="color: #008080;">javascript_semantics</span> |
<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: #008080;">constant</span> <span style="color: #000000;">dt</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">0.1</span> |
||
Line 2,793: | Line 2,793: | ||
<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: #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> |
<span style="color: #008080;">end</span> <span style="color: #008080;">for</span> |
||
<!--</ |
<!--</syntaxhighlight>--> |
||
{{out}} |
{{out}} |
||
<pre> |
<pre> |
||
Line 2,812: | Line 2,812: | ||
=={{header|PL/I}}== |
=={{header|PL/I}}== |
||
<syntaxhighlight lang="pl/i"> |
|||
<lang PL/I> |
|||
Runge_Kutta: procedure options (main); /* 10 March 2014 */ |
Runge_Kutta: procedure options (main); /* 10 March 2014 */ |
||
declare (y, dy1, dy2, dy3, dy4) float (18); |
declare (y, dy1, dy2, dy3, dy4) float (18); |
||
Line 2,838: | Line 2,838: | ||
end Runge_kutta; |
end Runge_kutta; |
||
</syntaxhighlight> |
|||
</lang> |
|||
{{out}} |
{{out}} |
||
<pre> |
<pre> |
||
Line 2,858: | Line 2,858: | ||
{{works with|PowerShell|4.0}} |
{{works with|PowerShell|4.0}} |
||
<syntaxhighlight lang="powershell"> |
|||
<lang PowerShell> |
|||
function Runge-Kutta (${function:F}, ${function:y}, $y0, $t0, $dt, $tEnd) { |
function Runge-Kutta (${function:F}, ${function:y}, $y0, $t0, $dt, $tEnd) { |
||
function RK ($tn,$yn) { |
function RK ($tn,$yn) { |
||
Line 2,895: | Line 2,895: | ||
$tEnd = 10 |
$tEnd = 10 |
||
Runge-Kutta F y $y0 $t0 $dt $tEnd |
Runge-Kutta F y $y0 $t0 $dt $tEnd |
||
</syntaxhighlight> |
|||
</lang> |
|||
<b>Output:</b> |
<b>Output:</b> |
||
<pre> |
<pre> |
||
Line 2,915: | Line 2,915: | ||
=={{header|PureBasic}}== |
=={{header|PureBasic}}== |
||
{{trans|BBC Basic}} |
{{trans|BBC Basic}} |
||
< |
<syntaxhighlight lang="purebasic">EnableExplicit |
||
Define.i i |
Define.i i |
||
Define.d y=1.0, k1=0.0, k2=0.0, k3=0.0, k4=0.0, t=0.0 |
Define.d y=1.0, k1=0.0, k2=0.0, k3=0.0, k4=0.0, t=0.0 |
||
Line 2,933: | Line 2,933: | ||
Print("Press return to exit...") : Input() |
Print("Press return to exit...") : Input() |
||
EndIf |
EndIf |
||
End</ |
End</syntaxhighlight> |
||
{{out}} |
{{out}} |
||
<pre>y( 0) = 1.0000 Error = 0.0000000000 |
<pre>y( 0) = 1.0000 Error = 0.0000000000 |
||
Line 2,949: | Line 2,949: | ||
=={{header|Python}}== |
=={{header|Python}}== |
||
< |
<syntaxhighlight lang="python">from math import sqrt |
||
def rk4(f, x0, y0, x1, n): |
def rk4(f, x0, y0, x1, n): |
||
Line 2,983: | Line 2,983: | ||
8.0 288.99997 -3.1565e-05 |
8.0 288.99997 -3.1565e-05 |
||
9.0 451.56246 -4.0723e-05 |
9.0 451.56246 -4.0723e-05 |
||
10.0 675.99995 -5.0983e-05</ |
10.0 675.99995 -5.0983e-05</syntaxhighlight> |
||
=={{header|R}}== |
=={{header|R}}== |
||
< |
<syntaxhighlight lang="r">rk4 <- function(f, x0, y0, x1, n) { |
||
vx <- double(n + 1) |
vx <- double(n + 1) |
||
vy <- double(n + 1) |
vy <- double(n + 1) |
||
Line 3,018: | Line 3,018: | ||
[9,] 8 288.999968 -3.156520e-05 |
[9,] 8 288.999968 -3.156520e-05 |
||
[10,] 9 451.562459 -4.072316e-05 |
[10,] 9 451.562459 -4.072316e-05 |
||
[11,] 10 675.999949 -5.098329e-05</ |
[11,] 10 675.999949 -5.098329e-05</syntaxhighlight> |
||
=={{header|Racket}}== |
=={{header|Racket}}== |
||
Line 3,025: | Line 3,025: | ||
The Runge-Kutta method |
The Runge-Kutta method |
||
< |
<syntaxhighlight lang="racket"> |
||
(define (RK4 F δt) |
(define (RK4 F δt) |
||
(λ (t y) |
(λ (t y) |
||
Line 3,034: | Line 3,034: | ||
(list (+ t δt) |
(list (+ t δt) |
||
(+ y (* 1/6 (+ δy1 (* 2 δy2) (* 2 δy3) δy4)))))) |
(+ y (* 1/6 (+ δy1 (* 2 δy2) (* 2 δy3) δy4)))))) |
||
</syntaxhighlight> |
|||
</lang> |
|||
The method modifier which divides each time-step into ''n'' sub-steps: |
The method modifier which divides each time-step into ''n'' sub-steps: |
||
< |
<syntaxhighlight lang="racket"> |
||
(define ((step-subdivision n method) F h) |
(define ((step-subdivision n method) F h) |
||
(λ (x . y) (last (ODE-solve F (cons x y) |
(λ (x . y) (last (ODE-solve F (cons x y) |
||
Line 3,043: | Line 3,043: | ||
#:step (/ h n) |
#:step (/ h n) |
||
#:method method)))) |
#:method method)))) |
||
</syntaxhighlight> |
|||
</lang> |
|||
Usage: |
Usage: |
||
< |
<syntaxhighlight lang="racket"> |
||
(define (F t y) (* t (sqrt y))) |
(define (F t y) (* t (sqrt y))) |
||
Line 3,057: | Line 3,057: | ||
(match-define (list t y) s) |
(match-define (list t y) s) |
||
(printf "t=~a\ty=~a\terror=~a\n" t y (- y (exact-solution t)))) |
(printf "t=~a\ty=~a\terror=~a\n" t y (- y (exact-solution t)))) |
||
</syntaxhighlight> |
|||
</lang> |
|||
{{out}} |
{{out}} |
||
<pre> |
<pre> |
||
Line 3,075: | Line 3,075: | ||
Graphical representation: |
Graphical representation: |
||
< |
<syntaxhighlight lang="racket"> |
||
> (require plot) |
> (require plot) |
||
> (plot (list (function exact-solution 0 10 #:label "Exact solution") |
> (plot (list (function exact-solution 0 10 #:label "Exact solution") |
||
(points numeric-solution #:label "Runge-Kutta method")) |
(points numeric-solution #:label "Runge-Kutta method")) |
||
#:x-label "t" #:y-label "y(t)") |
#:x-label "t" #:y-label "y(t)") |
||
</syntaxhighlight> |
|||
</lang> |
|||
[[File:runge-kutta.png]] |
[[File:runge-kutta.png]] |
||
Line 3,086: | Line 3,086: | ||
(formerly Perl 6) |
(formerly Perl 6) |
||
{{Works with|rakudo|2016.03}} |
{{Works with|rakudo|2016.03}} |
||
<lang |
<syntaxhighlight lang="raku" line>sub runge-kutta(&yp) { |
||
return -> \t, \y, \δt { |
return -> \t, \y, \δt { |
||
my $a = δt * yp( t, y ); |
my $a = δt * yp( t, y ); |
||
Line 3,106: | Line 3,106: | ||
printf "y(%2d) = %12f ± %e\n", $t, $y, abs($y - ($t**2 + 4)**2 / 16) |
printf "y(%2d) = %12f ± %e\n", $t, $y, abs($y - ($t**2 + 4)**2 / 16) |
||
if $t %% 1; |
if $t %% 1; |
||
}</ |
}</syntaxhighlight> |
||
{{out}} |
{{out}} |
||
<pre>y( 0) = 1.000000 ± 0.000000e+00 |
<pre>y( 0) = 1.000000 ± 0.000000e+00 |
||
Line 3,129: | Line 3,129: | ||
< |
<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*/ |
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 */ |
x0= 0; x1= 10; dx= .1 /*define variables: X0 X1 DX */ |
||
Line 3,156: | Line 3,156: | ||
numeric digits; parse value format(x,2,1,,0) 'E0' with g 'E' _ .; g=g * .5'e'_ % 2 |
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 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</ |
do k=j+5 to 0 by -1; numeric digits m.k; g=(g+x/g)*.5; end /*k*/; return g</syntaxhighlight> |
||
Programming note: the '''fmt''' function is used to |
Programming note: the '''fmt''' function is used to |
||
align the output with attention paid to the different ways some |
align the output with attention paid to the different ways some |
||
Line 3,194: | Line 3,194: | ||
=={{header|Ring}}== |
=={{header|Ring}}== |
||
< |
<syntaxhighlight lang="ring"> |
||
decimals(8) |
decimals(8) |
||
y = 1.0 |
y = 1.0 |
||
Line 3,208: | Line 3,208: | ||
y += 0.1 * (k1 + 2 * (k2 + k3) + k4) / 6 |
y += 0.1 * (k1 + 2 * (k2 + k3) + k4) / 6 |
||
next |
next |
||
</syntaxhighlight> |
|||
</lang> |
|||
Output: |
Output: |
||
Line 3,226: | Line 3,226: | ||
=={{header|Ruby}}== |
=={{header|Ruby}}== |
||
< |
<syntaxhighlight lang="ruby">def calc_rk4(f) |
||
return ->(t,y,dt){ |
return ->(t,y,dt){ |
||
->(dy1 ){ |
->(dy1 ){ |
||
Line 3,252: | Line 3,252: | ||
printf("y(%4.1f)\t= %12.6f \t error: %12.6e\n",t,y,find_error(t,y)) if is_whole?(t) |
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) |
t, y = t + DT, y + dy.call(t,y,DT) |
||
end</ |
end</syntaxhighlight> |
||
{{Out}} |
{{Out}} |
||
<pre> |
<pre> |
||
Line 3,269: | Line 3,269: | ||
=={{header|Run BASIC}}== |
=={{header|Run BASIC}}== |
||
< |
<syntaxhighlight lang="runbasic">y = 1 |
||
while t <= 10 |
while t <= 10 |
||
k1 = t * sqr(y) |
k1 = t * sqr(y) |
||
Line 3,280: | Line 3,280: | ||
t = t + .1 |
t = t + .1 |
||
wend |
wend |
||
end</ |
end</syntaxhighlight> |
||
{{out}} |
{{out}} |
||
<pre>y( 0) = 1.0000000 Error =0 |
<pre>y( 0) = 1.0000000 Error =0 |
||
Line 3,297: | Line 3,297: | ||
=={{header|Rust}}== |
=={{header|Rust}}== |
||
This is a translation of the javascript solution with some minor differences. |
This is a translation of the javascript solution with some minor differences. |
||
< |
<syntaxhighlight lang="rust">fn runge_kutta4(fx: &dyn Fn(f64, f64) -> f64, x: f64, y: f64, dx: f64) -> f64 { |
||
let k1 = dx * fx(x, y); |
let k1 = dx * fx(x, y); |
||
let k2 = dx * fx(x + dx / 2.0, y + k1 / 2.0); |
let k2 = dx * fx(x + dx / 2.0, y + k1 / 2.0); |
||
Line 3,330: | Line 3,330: | ||
x = ((x * 10.0) + (step * 10.0)) / 10.0; |
x = ((x * 10.0) + (step * 10.0)) / 10.0; |
||
} |
} |
||
}</ |
}</syntaxhighlight> |
||
<pre> |
<pre> |
||
y(0): 1.0000000000 0E0 |
y(0): 1.0000000000 0E0 |
||
Line 3,346: | Line 3,346: | ||
=={{header|Scala}}== |
=={{header|Scala}}== |
||
< |
<syntaxhighlight lang="scala">object Main extends App { |
||
val f = (t: Double, y: Double) => t * Math.sqrt(y) // Runge-Kutta solution |
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 |
val g = (t: Double) => Math.pow(t * t + 4, 2) / 16 // Exact solution |
||
Line 3,372: | Line 3,372: | ||
} |
} |
||
} |
} |
||
}</ |
}</syntaxhighlight> |
||
<pre> |
<pre> |
||
y( 0.0) = 1.00000000 Error: 0.00000e+00 |
y( 0.0) = 1.00000000 Error: 0.00000e+00 |
||
Line 3,389: | Line 3,389: | ||
=={{header|Sidef}}== |
=={{header|Sidef}}== |
||
{{trans|Raku}} |
{{trans|Raku}} |
||
< |
<syntaxhighlight lang="ruby">func runge_kutta(yp) { |
||
func (t, y, δt) { |
func (t, y, δt) { |
||
var a = (δt * yp(t, y)); |
var a = (δt * yp(t, y)); |
||
Line 3,409: | Line 3,409: | ||
y += δy(t, y, δt); |
y += δy(t, y, δt); |
||
t += δt; |
t += δt; |
||
}</ |
}</syntaxhighlight> |
||
{{out}} |
{{out}} |
||
<pre> |
<pre> |
||
Line 3,426: | Line 3,426: | ||
=={{header|Standard ML}}== |
=={{header|Standard ML}}== |
||
< |
<syntaxhighlight lang="sml">fun step y' (tn,yn) dt = |
||
let |
let |
||
val dy1 = dt * y'(tn,yn) |
val dy1 = dt * y'(tn,yn) |
||
Line 3,466: | Line 3,466: | ||
(* Run the suggested test case *) |
(* Run the suggested test case *) |
||
val () = test 0.0 1.0 0.1 101 10 testy testy'</ |
val () = test 0.0 1.0 0.1 101 10 testy testy'</syntaxhighlight> |
||
{{out}} |
{{out}} |
||
<pre>Time: 0.0 |
<pre>Time: 0.0 |
||
Line 3,524: | Line 3,524: | ||
=={{header|Stata}}== |
=={{header|Stata}}== |
||
< |
<syntaxhighlight lang="stata">function rk4(f, t0, y0, t1, n) { |
||
h = (t1-t0)/(n-1) |
h = (t1-t0)/(n-1) |
||
a = J(n, 2, 0) |
a = J(n, 2, 0) |
||
Line 3,564: | Line 3,564: | ||
10 | 9 451.5624593 -.0000407232 | |
10 | 9 451.5624593 -.0000407232 | |
||
11 | 10 675.999949 -.0000509833 | |
11 | 10 675.999949 -.0000509833 | |
||
+----------------------------------------------+</ |
+----------------------------------------------+</syntaxhighlight> |
||
=={{header|Swift}}== |
=={{header|Swift}}== |
||
{{trans|C}} |
{{trans|C}} |
||
< |
<syntaxhighlight lang="swift">import Foundation |
||
func rk4(dx: Double, x: Double, y: Double, f: (Double, Double) -> Double) -> Double { |
func rk4(dx: Double, x: Double, y: Double, f: (Double, Double) -> Double) -> Double { |
||
Line 3,605: | Line 3,605: | ||
print(String(format: "%2g %11.6g %11.5g", x, y[i], y[i]/y2 - 1)) |
print(String(format: "%2g %11.6g %11.5g", x, y[i], y[i]/y2 - 1)) |
||
}</ |
}</syntaxhighlight> |
||
{{out}} |
{{out}} |
||
Line 3,623: | Line 3,623: | ||
=={{header|Tcl}}== |
=={{header|Tcl}}== |
||
< |
<syntaxhighlight lang="tcl">package require Tcl 8.5 |
||
# Hack to bring argument function into expression |
# Hack to bring argument function into expression |
||
Line 3,655: | Line 3,655: | ||
printvals $t $y |
printvals $t $y |
||
} |
} |
||
}</ |
}</syntaxhighlight> |
||
{{out}} |
{{out}} |
||
<pre>y(0.0) = 1.00000000 Error: 0.00000000e+00 |
<pre>y(0.0) = 1.00000000 Error: 0.00000000e+00 |
||
Line 3,672: | Line 3,672: | ||
{{trans|Kotlin}} |
{{trans|Kotlin}} |
||
{{libheader|Wren-fmt}} |
{{libheader|Wren-fmt}} |
||
< |
<syntaxhighlight lang="ecmascript">import "/fmt" for Fmt |
||
var rungeKutta4 = Fn.new { |t0, tz, dt, y, yd| |
var rungeKutta4 = Fn.new { |t0, tz, dt, y, yd| |
||
Line 3,701: | Line 3,701: | ||
} |
} |
||
var yd = Fn.new { |t, yt| t * yt.sqrt } |
var yd = Fn.new { |t, yt| t * yt.sqrt } |
||
rungeKutta4.call(0, 10, 0.1, y, yd)</ |
rungeKutta4.call(0, 10, 0.1, y, yd)</syntaxhighlight> |
||
{{out}} |
{{out}} |
||
Line 3,722: | Line 3,722: | ||
=={{header|zkl}}== |
=={{header|zkl}}== |
||
{{trans|OCaml}} |
{{trans|OCaml}} |
||
< |
<syntaxhighlight lang="zkl">fcn yp(t,y) { t * y.sqrt() } |
||
fcn exact(t){ u:=0.25*t*t + 1.0; u*u } |
fcn exact(t){ u:=0.25*t*t + 1.0; u*u } |
||
Line 3,737: | Line 3,737: | ||
print("t = %f,\ty = %f,\terr = %g\n".fmt(t,y,(y - exact(t)).abs())); |
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 |
if(n < 102) return(loop(h,(n+1),rk4_step(T(y,t),h))) //tail recursion |
||
}</ |
}</syntaxhighlight> |
||
{{out}} |
{{out}} |
||
<pre> |
<pre> |