Runge-Kutta method: Difference between revisions

This equation has an exact solution:

:<math>y(t) = \tfrac{1}{16}(t^2 +4)^2</math>

;Task

Demonstrate the commonly used explicit &nbsp; [[wp:Runge–Kutta_methods#Common_fourth-order_Runge.E2.80.93Kutta_method|fourth-order Runge–Kutta method]] &nbsp; 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:

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

 

=={{header|Ada}}==

with Ada.Numerics.Generic_Elementary_Functions;

procedure RungeKutta is

Runge (yprime'Access, t_arr, y_arr, dt);

Print (t_arr, y_arr, 10);

{{out}}

<pre>y(0.0) = 1.00000000 Error: 0.00000E+00

 

=={{header|ALGOL 68}}==

BEGIN

PROC rk4 = (PROC (REAL, REAL) REAL f, REAL y, x, dx) REAL :

OD

END

{{out}}

<pre>

9.0000000 451.5625000 451.5624593 -9.0183e-08

10.0000000 676.0000000 675.9999490 -7.5419e-08

</pre>

 

=={{header|APL}}==

∇RK4[⎕]∇

∇

[2] ⎕←'T' 'RK4 Y' 'ERROR'⍪TABLE,TABLE[;2]-{((4+⍵*2)*2)÷16}TABLE[;1]

∇

{{out}}

<pre>

 

=={{header|AWK}}==

# syntax: GAWK -f RUNGE-KUTTA_METHOD.AWK

# converted from BBC BASIC

exit(0)

}

{{out}}

<pre>

 

=={{header|BASIC}}==

==={{header|BBC BASIC}}===

FOR i% = 0 TO 100

t = i% / 10

k4 = (t + 0.10) * SQR(y + 0.10 * k3)

y += 0.1 * (k1 + 2 * (k2 + k3) + k4) / 6

{{out}}

<pre>y(0) = 1 Error = 0

 

==={{header|IS-BASIC}}===

110 LET Y=1

120 FOR T=0 TO 10 STEP .1

170 LET K4=(T+.1)*SQR(Y+.1*K3)

180 LET Y=Y+.1*(K1+2*(K2+K3)+K4)/6

 

=={{header|C}}==

#include <stdlib.h>

#include <math.h>

 

return 0;

{{out}} (errors are relative)

<pre>

 

=={{header|C sharp|C#}}==

using System;

 

}

}

 

=={{header|C++}}==

Using Lambdas

*

*/

# include <iostream>

# include <math.h>

auto rk4(double f(double, double))

{

}

 

int main(void)

{

 

 

=={{header|Common Lisp}}==

 

(let ((h (float (/ (- x-end x) n) 1d0))

k1 k2 k3 k4)

(7.999999999999988d0 288.9999684347983d0 -3.156520000402452d-5)

(8.999999999999984d0 451.56245927683887d0 -4.072315812209126d-5)

 

=={{header|Crystal}}==

{{trans|Run Basic and Ruby output}}

while t <= 10

k1 = t * Math.sqrt(y)

y += 0.1 * (k1 + 2 * (k2 + k3) + k4) / 6

t += 0.1

 

{{out}}

=={{header|D}}==

{{trans|Ada}}

 

alias FP = real;

t_arr[i], y_arr[i],

calc_err(t_arr[i], y_arr[i]));

{{out}}

<pre>y(0.0) = 1.00000000 Error: 0

 

=={{header|Dart}}==

 

num RungeKutta4(Function f, num t, num y, num dt){

t += dt;

}

{{out}}

<pre>

y(9.00) = 451.56245928 Error = 9.0182772312e-8

y(10.0) = 675.99994902 Error = 7.5419063100e-8

</pre>

 

=={{header|ERRE}}==

PROGRAM RUNGE_KUTTA

 

Y+=DELTA_T*(K1+2*(K2+K3)+K4)/6

END FOR

{{out}}

<pre>

</pre>

 

{{works with|F# interactive (fsi.exe)}}

open System

 

RungeKutta4 0.0 1.0 10.0 0.1

|> Seq.filter (fun (t,y) -> t % 1.0 = 0.0 )

 

{{out}}

 

=={{header|Fortran}}==

implicit none

{{out}}

<pre>

=={{header|FreeBASIC}}==

{{trans|BBC BASIC}}

' compile with: fbc -s console

' translation of BBC BASIC

Print : Print "hit any key to end program"

Sleep

{{out}}

<pre>y(0) = 1 Error = 0

 

=={{header|FutureBasic}}==

 

 

 

h = 0.1

y = 1

for i = 0 to 100

next

Output:

<pre>

=={{header|Go}}==

{{works with|Go1}}

 

import (

func printErr(t, y float64) {

fmt.Printf("y(%.1f) = %f Error: %e\n", t, y, math.Abs(actual(t)-y))

{{out}}

<pre>

 

=={{header|Groovy}}==

class Runge_Kutta{

static void main(String[] args){

static def dy(def x){return x*0.1;}

}

{{out}}

<pre>

y(9.0)=451.56245927683966 Error:4.07231603389846E-5

y(10.0)=675.9999490167097 Error:5.098329029351589E-5

</pre>

 

Using GHC 7.4.1.

 

:: Floating a

=> a -> a -> a

(print . (\(x, y) -> (truncate x, y, fy x - y)))

(filter (\(x, _) -> 0 == mod (truncate $ 10 * x) 10) $

 

Example executed in GHCi:

(0,1.0,0.0)

(1,1.5624998542781088,1.4572189122041834e-7)

(8,288.99996843479926,3.1565204153594095e-5)

(9,451.562459276841,4.0723166534917254e-5)

 

(See [[Euler method#Haskell]] for implementation of simple general ODE-solver)

 

Or, disaggregated a little, and expressed in terms of a single scanl:

rk4 y x dx =

let f x y = x * sqrt y

where

justifyLeft n c s = take n (s ++ replicate n c)

{{Out}}

<pre>y (0) = 1.0 ±0.0

=={{header|J}}==

'''Solution:'''

NB. y is: y(ta) , ta , tb , tstep

NB. u is: function to solve

end.

T ,. Y

'''Example:'''

fyp=: (* %:)/ NB. f'(t,y)

report_whole=: (10 * i. >:10)&{ NB. report at whole-numbered t values

8 289 _3.15652e_5

9 451.562 _4.07232e_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

ks=. (x * [: u y + (* x&,))/\. tableau

({:y) + 6 %~ +/ 1 2 2 1 * ks

 

Use:

Translation of [[Runge-Kutta_method#Ada|Ada]] via [[Runge-Kutta_method#D|D]]

{{works with|Java|8}}

import java.util.function.BiFunction;

 

calc_err(t_arr[i], y_arr[i]));

}

 

<pre>y(0,0) = 1,00000000 Error: 0,000000

=={{header|JavaScript}}==

===ES5===

function rk4(y, x, dx, f) {

var k1 = dx * f(x, y),

steps += 1;

}

{{out}}

<pre>

 

===ES6===

'use strict';

 

// MAIN ---

return main();

{{Out}}

<pre>y (0) = 1.0 ± 0e+0

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;

def while(cond; update):

def _while: if cond then ., (update | _while) else empty end;

 

===The Example Differential Equation and its Exact Solution===

def yprime: .[0] * (.[1] | sqrt);

. as $t

| (( $t*$t) + 4 )

 

===dy/dt===

 

'''Generic filters:'''

def round(n):

(if . < 0 then -1 else 1 end) as $s

 

# Is the input an integer?

 

'''dy(f)'''

 

# Input: [t, y]; yp is a filter that accepts [t,y] as input

 

# Input: [t,y]

''' Example''':

[0,1]

| while( .[0] <= 10;

"y(\($t|round(1))) = \($y|round(10000)) ± \( ($t|actual) - $y | abs)"

else empty

{{out}}

y(0) = 1 ± 0

y(1) = 1.5625 ± 1.4572189210859676e-07

real 0m0.048s

user 0m0.013s

 

===newRK4Step===

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

 

# main(t0; y0; tFinal; dtPrint)

{{out}}

y(0) = 1 with error: 0

y(1) = 1.562499854278108 with error: 1.4572189210859676e-07

real 0m0.023s

user 0m0.014s

 

=={{header|Julia}}==

{{works with|Julia|0.6}}

 

{{trans|Python}}

theoric(t) = (t ^ 2 + 4.0) ^ 2 / 16.0

 

y += δy(t, y, δt)

t += δt

 

{{out}}

y(10.0) = 675.999949 error: 5.098329e-05</pre>

 

{{trans|Python}}

vx = Vector{Float64}(undef, n + 1)

vy = Vector{Float64}(undef, n + 1)

for (x, y) in Iterators.take(zip(vx, vy), 10)

@printf("%4.1f %10.5f %+12.4e\n", x, y, y - theoric(x))

 

=={{header|Kotlin}}==

 

typealias Y = (Double) -> Double

val yd = fun(t: Double, yt: Double) = t * Math.sqrt(yt)

rungeKutta4(0.0, 10.0, 0.1, y, yd)

 

{{out}}

 

=={{header|Liberty BASIC}}==

'[RC] Runge-Kutta method

'initial conditions

exactY=(x^2 + 4)^2 / 16

end function

{{Out}}

<pre>

</pre>

 

DSolve[{y'[t] == t*Sqrt[y[t]], y[0] == 1}, y, t]

Table[{t, 1/16 (4 + t^2)^2}, {t, 0, 10}]

solution = NestList[phi, {0, 1}, 101];

Table[{y[[1]], y[[2]], Abs[y[[2]] - 1/16 (y[[1]]^2 + 4)^2]},

 

=={{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

y(k+1) = y(k)+(dy1+2*dy2+2*dy3+dy4)/6;

end

{{out}}

<pre>

 

=={{header|Maxima}}==

'diff(y, x) = x * sqrt(y);

ode2(%, y, x);

s: map(lambda([x], (x^2 + 4)^2 / 16), x)$

 

 

=={{header|МК-61/52}}==

 

=={{header|Nim}}==

 

proc fn(t, y: float): float =

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 dy4 = step * fn(cur_t + step, cur_y + dy3)

 

 

{{out}}

<pre>y(0.0) = 1.0, error = 0.0

 

=={{header|Objeck}}==

function : Main(args : String[]) ~ Nil {

x0 := 0.0; x1 := 10.0; dx := .1;

return x * y->SquareRoot();

}

 

Output:

 

=={{header|OCaml}}==

let exact t = let u = 0.25*.t*.t +. 1.0 in u*.u

 

if n < 102 then loop h (n+1) (rk4_step (y,t) h)

 

{{out}}

<pre>t = 0.000000, y = 1.000000, err = 0

 

=={{header|Octave}}==

#Applying the Runge-Kutta method (This code must be implement on a different file than the main one).

 

fprintf('%d \t %.5f \t %.5f \t %.4g \n',i,f(i),Yn(1+i*10),f(i)-Yn(1+i*10));

endfor

{{out}}

<pre>

=={{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

)

};

{{out}}

<pre>x y rel. err.

This code has been compiled using Free Pascal 2.6.2.

 

 

uses sysutils;

RungeKutta(@YPrime, tArr, yArr, dt);

Print(tArr, yArr, 10);

{{out}}

<pre>y( 0.0) = 1.00000000 Error: 0.00000E+000

 

Notice how we have to use sprintf to deal with floating point rounding. See perlfaq4.

my ($yp, $dt) = @_;

sub {

printf "y(%2.0f) = %12f ± %e\n", $t, $y, abs($y - ($t**2 + 4)**2 / 16)

if sprintf("%.4f", $t) =~ /0000$/;

 

{{out}}

=={{header|Phix}}==

{{trans|ERRE}}

{{out}}

<pre>

 

=={{header|PL/I}}==

Runge_Kutta: procedure options (main); /* 10 March 2014 */

declare (y, dy1, dy2, dy3, dy4) float (18);

 

end Runge_kutta;

{{out}}

<pre>

{{works with|PowerShell|4.0}}

 

function Runge-Kutta (${function:F}, ${function:y}, $y0, $t0, $dt, $tEnd) {

function RK ($tn,$yn) {

$tEnd = 10

Runge-Kutta F y $y0 $t0 $dt $tEnd

<b>Output:</b>

<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

Print("Press return to exit...") : Input()

EndIf

{{out}}

<pre>y( 0) = 1.0000 Error = 0.0000000000

 

=={{header|Python}}==

def rk4(f, x0, y0, x1, n):

8.0 288.99997 -3.1565e-05

9.0 451.56246 -4.0723e-05

 

=={{header|R}}==

 

vx <- double(n + 1)

vy <- double(n + 1)

[9,] 8 288.999968 -3.156520e-05

[10,] 9 451.562459 -4.072316e-05

 

=={{header|Racket}}==

 

The Runge-Kutta method

(define (RK4 F δt)

(λ (t y)

(list (+ t δt)

(+ y (* 1/6 (+ δy1 (* 2 δy2) (* 2 δy3) δy4))))))

 

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)

#:step (/ h n)

#:method method))))

 

Usage:

(define (F t y) (* t (sqrt y)))

 

(match-define (list t y) s)

(printf "t=~a\ty=~a\terror=~a\n" t y (- y (exact-solution t))))

{{out}}

<pre>

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)")

[[File:runge-kutta.png]]

 

(formerly Perl 6)

{{Works with|rakudo|2016.03}}

return -> \t, \y, \δt {

my $a = δt * yp( t, y );

printf "y(%2d) = %12f ± %e\n", $t, $y, abs($y - ($t**2 + 4)**2 / 16)

if $t %% 1;

{{out}}

<pre>y( 0) = 1.000000 ± 0.000000e+00

 

=={{header|REXX}}==

 

n=1 + (x1-x0) / dx

end /*m*/ /* [↑] use 4th order Runge─Kutta. */

say center('X', f, "═") center('Y', w+2, "═") center("relative error", w+8, '═') /*hdr*/

 

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

/*──────────────────────────────────────────────────────────────────────────────────────*/

/*──────────────────────────────────────────────────────────────────────────────────────*/

return y + (k1 + k2*2 + k3*2 + k4) / 6

/*──────────────────────────────────────────────────────────────────────────────────────*/

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

Programming note: &nbsp; the &nbsp; '''fmt''' &nbsp; function is used to

align the output with attention paid to the different ways some

 

=={{header|Ring}}==

decimals(8)

y = 1.0

y += 0.1 * (k1 + 2 * (k2 + k3) + k4) / 6

next

 

Output:

 

=={{header|Ruby}}==

return ->(t,y,dt){

->(dy1 ){

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)

{{Out}}

<pre>

 

=={{header|Run BASIC}}==

while t <= 10

k1 = t * sqr(y)

t = t + .1

wend

{{out}}

<pre>y( 0) = 1.0000000 Error =0

=={{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);

x = ((x * 10.0) + (step * 10.0)) / 10.0;

}

<pre>

y(0): 1.0000000000 0E0

 

=={{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

}

}

<pre>

y( 0.0) = 1.00000000 Error: 0.00000e+00

=={{header|Sidef}}==

{{trans|Raku}}

func (t, y, δt) {

var a = (δt * yp(t, y));

y += δy(t, y, δt);

t += δt;

{{out}}

<pre>

 

=={{header|Standard ML}}==

let

val dy1 = dt * y'(tn,yn)

 

(* Run the suggested test case *)

{{out}}

<pre>Time: 0.0

 

=={{header|Stata}}==

h = (t1-t0)/(n-1)

a = J(n, 2, 0)

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 {

 

print(String(format: "%2g %11.6g %11.5g", x, y[i], y[i]/y2 - 1))

 

{{out}}

 

=={{header|Tcl}}==

 

# Hack to bring argument function into expression

printvals $t $y

}

{{out}}

<pre>y(0.0) = 1.00000000 Error: 0.00000000e+00

y(9.0) = 451.56245928 Error: 4.07231581e-05

y(10.0) = 675.99994902 Error: 5.09832864e-05</pre>

 

=={{header|zkl}}==

{{trans|OCaml}}

fcn exact(t){ u:=0.25*t*t + 1.0; u*u }

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>