Runge-Kutta method
Given the example Differential equation:
With initial condition:
- and
This equation has an exact solution:
- Task
Demonstrate the commonly used explicit fourth-order Runge–Kutta method to solve the above differential equation.
- Solve the given differential equation over the range with a step value of (101 total points, the first being given)
- Print the calculated values of at whole numbered 's () along with error as compared to the exact solution.
- Method summary
Starting with a given and calculate:
then:
Ada
<lang Ada>with Ada.Text_IO; use Ada.Text_IO; with Ada.Numerics.Generic_Elementary_Functions; procedure RungeKutta is
type Floaty is digits 15;
type Floaty_Array is array (Natural range <>) of Floaty;
package FIO is new Ada.Text_IO.Float_IO(Floaty); use FIO;
type Derivative is access function(t, y : Floaty) return Floaty;
package Math is new Ada.Numerics.Generic_Elementary_Functions (Floaty);
function calc_err (t, calc : Floaty) return Floaty;
procedure Runge (yp_func : Derivative; t, y : in out Floaty_Array;
dt : Floaty) is
dy1, dy2, dy3, dy4 : Floaty;
begin
for n in t'First .. t'Last-1 loop
dy1 := dt * yp_func(t(n), y(n));
dy2 := dt * yp_func(t(n) + dt / 2.0, y(n) + dy1 / 2.0);
dy3 := dt * yp_func(t(n) + dt / 2.0, y(n) + dy2 / 2.0);
dy4 := dt * yp_func(t(n) + dt, y(n) + dy3);
t(n+1) := t(n) + dt;
y(n+1) := y(n) + (dy1 + 2.0 * (dy2 + dy3) + dy4) / 6.0;
end loop;
end Runge;
procedure Print (t, y : Floaty_Array; modnum : Positive) is begin
for i in t'Range loop
if i mod modnum = 0 then
Put("y("); Put (t(i), Exp=>0, Fore=>0, Aft=>1);
Put(") = "); Put (y(i), Exp=>0, Fore=>0, Aft=>8);
Put(" Error:"); Put (calc_err(t(i),y(i)), Aft=>5);
New_Line;
end if;
end loop;
end Print;
function yprime (t, y : Floaty) return Floaty is begin
return t * Math.Sqrt (y);
end yprime;
function calc_err (t, calc : Floaty) return Floaty is
actual : constant Floaty := (t**2 + 4.0)**2 / 16.0;
begin return abs(actual-calc);
end calc_err;
dt : constant Floaty := 0.10;
N : constant Positive := 100;
t_arr, y_arr : Floaty_Array(0 .. N);
begin
t_arr(0) := 0.0; y_arr(0) := 1.0; Runge (yprime'Access, t_arr, y_arr, dt); Print (t_arr, y_arr, 10);
end RungeKutta;</lang>
- Output:
y(0.0) = 1.00000000 Error: 0.00000E+00 y(1.0) = 1.56249985 Error: 1.45722E-07 y(2.0) = 3.99999908 Error: 9.19479E-07 y(3.0) = 10.56249709 Error: 2.90956E-06 y(4.0) = 24.99999377 Error: 6.23491E-06 y(5.0) = 52.56248918 Error: 1.08197E-05 y(6.0) = 99.99998341 Error: 1.65946E-05 y(7.0) = 175.56247648 Error: 2.35177E-05 y(8.0) = 288.99996843 Error: 3.15652E-05 y(9.0) = 451.56245928 Error: 4.07232E-05 y(10.0) = 675.99994902 Error: 5.09833E-05
C
<lang c>#include <stdio.h>
- include <stdlib.h>
- include <math.h>
double rk4(double(*f)(double, double), double dx, double x, double y) { double k1 = dx * f(x, y), k2 = dx * f(x + dx / 2, y + k1 / 2), k3 = dx * f(x + dx / 2, y + k2 / 2), k4 = dx * f(x + dx, y + k3); return y + (k1 + 2 * k2 + 2 * k3 + k4) / 6; }
double rate(double x, double y) { return x * sqrt(y); }
int main(void) { double *y, x, y2; double x0 = 0, x1 = 10, dx = .1; int i, n = 1 + (x1 - x0)/dx; y = malloc(sizeof(double) * n);
for (y[0] = 1, i = 1; i < n; i++) y[i] = rk4(rate, dx, x0 + dx * (i - 1), y[i-1]);
printf("x\ty\trel. err.\n------------\n"); for (i = 0; i < n; i += 10) { x = x0 + dx * i; y2 = pow(x * x / 4 + 1, 2); printf("%g\t%g\t%g\n", x, y[i], y[i]/y2 - 1); }
return 0;
}</lang>output (errors are relative)
x y rel. err. ------------ 0 1 0 1 1.5625 -9.3262e-08 2 4 -2.2987e-07 3 10.5625 -2.75462e-07 4 25 -2.49396e-07 5 52.5625 -2.05844e-07 6 100 -1.65946e-07 7 175.562 -1.33956e-07 8 289 -1.09222e-07 9 451.562 -9.01828e-08 10 676 -7.54191e-08
D
<lang d>import std.stdio, std.math, std.typecons;
alias real FP; alias Typedef!(FP[101]) FPs;
void runge(in FP function(in FP, in FP) pure nothrow yp_func,
ref FPs t, ref FPs y, in FP dt) pure nothrow {
foreach (n; 0 .. t.length - 1) {
FP dy1 = dt * yp_func(t[n], y[n]);
FP dy2 = dt * yp_func(t[n] + dt / 2.0, y[n] + dy1 / 2.0);
FP dy3 = dt * yp_func(t[n] + dt / 2.0, y[n] + dy2 / 2.0);
FP dy4 = dt * yp_func(t[n] + dt, y[n] + dy3);
t[n + 1] = t[n] + dt;
y[n + 1] = y[n] + (dy1 + 2.0 * (dy2 + dy3) + dy4) / 6.0;
}
}
FP calc_err(in FP t, in FP calc) pure nothrow {
immutable FP actual = (t ^^ 2 + 4.0) ^^ 2 / 16.0; return abs(actual - calc);
}
void main() {
enum FP dt = 0.10; FPs t_arr, y_arr;
t_arr[0] = 0.0; y_arr[0] = 1.0; runge((t, y) => t * sqrt(y), t_arr, y_arr, dt);
foreach (i; 0 .. t_arr.length)
if (i % 10 == 0)
writefln("y(%.1f) = %.8f Error: %.6g",
t_arr[i], y_arr[i],
calc_err(t_arr[i], y_arr[i]));
}</lang>
- Output:
y(0.0) = 1.00000000 Error: 0 y(1.0) = 1.56249985 Error: 1.45722e-07 y(2.0) = 3.99999908 Error: 9.19479e-07 y(3.0) = 10.56249709 Error: 2.90956e-06 y(4.0) = 24.99999377 Error: 6.23491e-06 y(5.0) = 52.56248918 Error: 1.08197e-05 y(6.0) = 99.99998341 Error: 1.65946e-05 y(7.0) = 175.56247648 Error: 2.35177e-05 y(8.0) = 288.99996843 Error: 3.15652e-05 y(9.0) = 451.56245928 Error: 4.07232e-05 y(10.0) = 675.99994902 Error: 5.09833e-05
Go
<lang go>package main
import (
"fmt" "math"
)
type ypFunc func(t, y float64) float64 type ypStepFunc func(t, y, dt float64) float64
// newRKStep takes a function representing a differential equation // and returns a function that performs a single step of the forth-order // Runge-Kutta method. func newRK4Step(yp ypFunc) ypStepFunc {
return func(t, y, dt float64) float64 {
dy1 := dt * yp(t, y)
dy2 := dt * yp(t+dt/2, y+dy1/2)
dy3 := dt * yp(t+dt/2, y+dy2/2)
dy4 := dt * yp(t+dt, y+dy3)
return y + (dy1+2*(dy2+dy3)+dy4)/6
}
}
// example differential equation func yprime(t, y float64) float64 {
return t * math.Sqrt(y)
}
// exact solution of example func actual(t float64) float64 {
t = t*t + 4 return t * t / 16
}
func main() {
t0, tFinal := 0, 10 // task specifies times as integers, dtPrint := 1 // and to print at whole numbers. y0 := 1. // initial y. dtStep := .1 // step value.
t, y := float64(t0), y0
ypStep := newRK4Step(yprime)
for t1 := t0 + dtPrint; t1 <= tFinal; t1 += dtPrint {
printErr(t, y) // print intermediate result
for steps := int(float64(dtPrint)/dtStep + .5); steps > 1; steps-- {
y = ypStep(t, y, dtStep)
t += dtStep
}
y = ypStep(t, y, float64(t1)-t) // adjust step to integer time
t = float64(t1)
}
printErr(t, y) // print final result
}
func printErr(t, y float64) {
fmt.Printf("y(%.1f) = %f Error: %e\n", t, y, math.Abs(actual(t)-y))
}</lang>
- Output:
y(0.0) = 1.000000 Error: 0.000000e+00 y(1.0) = 1.562500 Error: 1.457219e-07 y(2.0) = 3.999999 Error: 9.194792e-07 y(3.0) = 10.562497 Error: 2.909562e-06 y(4.0) = 24.999994 Error: 6.234909e-06 y(5.0) = 52.562489 Error: 1.081970e-05 y(6.0) = 99.999983 Error: 1.659460e-05 y(7.0) = 175.562476 Error: 2.351773e-05 y(8.0) = 288.999968 Error: 3.156520e-05 y(9.0) = 451.562459 Error: 4.072316e-05 y(10.0) = 675.999949 Error: 5.098329e-05
Haskell
Using GHC 7.4.1.
<lang haskell>import Data.List
dv :: Floating a => a -> a -> a dv = (. sqrt). (*)
fy t = 1/16 * (4+t^2)^2
rk4 :: (Enum a, Fractional a)=> (a -> a -> a) -> a -> a -> a -> [(a,a)] rk4 fd y0 a h = zip ts $ scanl (flip fc) y0 ts where
ts = [a,h ..] fc t y = sum. (y:). zipWith (*) [1/6,1/3,1/3,1/6] $ scanl (\k f -> h * fd (t+f*h) (y+f*k)) (h * fd t y) [1/2,1/2,1]
task = mapM_ print
$ map (\(x,y)-> (truncate x,y,fy x - y)) $ filter (\(x,_) -> 0== mod (truncate $ 10*x) 10) $ take 101 $ rk4 dv 1.0 0 0.1</lang>
Example executed in GHCi: <lang haskell>*Main> task (0,1.0,0.0) (1,1.5624998542781088,1.4572189122041834e-7) (2,3.9999990805208006,9.194792029987298e-7) (3,10.562497090437557,2.909562461184123e-6) (4,24.999993765090654,6.234909399438493e-6) (5,52.56248918030265,1.0819697635611192e-5) (6,99.99998340540378,1.6594596999652822e-5) (7,175.56247648227165,2.3517730085131916e-5) (8,288.99996843479926,3.1565204153594095e-5) (9,451.562459276841,4.0723166534917254e-5) (10,675.9999490167125,5.098330132113915e-5)</lang>
J
Solution: <lang j>NB.*rk4 a Solve function using Runge-Kutta method NB. y is: y(ta) , ta , tb , tstep NB. u is: function to solve NB. eg: fyp rk4 1 0 10 0.1 rk4=: adverb define
'Y0 a b h'=. 4{. y
T=. a + i.@>:&.(%&h) b - a
Y=. Yt=. Y0
for_t. }: T do.
ty=. t,Yt
k1=. h * u ty
k2=. h * u ty + -: h,k1
k3=. h * u ty + -: h,k2
k4=. h * u ty + h,k3
Y=. Y, Yt=. Yt + (%6) * 1 2 2 1 +/@:* k1, k2, k3, k4
end.
T ,. Y )</lang> Example: <lang j> fy=: (%16) * [: *: 4 + *: NB. f(t,y)
fyp=: (* %:)/ NB. f'(t,y)
report_whole=: (10 * i. >:10)&{ NB. report at whole-numbered t values
report_err=: (, {: - [: fy {.)"1 NB. report errors
report_err report_whole fyp rk4 1 0 10 0.1 0 1 0 1 1.5625 _1.45722e_7 2 4 _9.19479e_7 3 10.5625 _2.90956e_6 4 25 _6.23491e_6 5 52.5625 _1.08197e_5 6 100 _1.65946e_5 7 175.562 _2.35177e_5 8 289 _3.15652e_5 9 451.562 _4.07232e_5
10 676 _5.09833e_5</lang>
Alternative solution:
The following solution replaces the for loop as well as the calculation of the increments (ks) with an accumulating suffix. <lang j>rk4=: adverb define
'Y0 a b h'=. 4{. y
T=. a + i.@>:&.(%&h) b-a
(,. [: h&(u nextY)@,/\. Y0 ,~ }.)&.|. T
)
NB. nextY a Calculate Yn+1 of a function using Runge-Kutta method NB. y is: 2-item numeric list of time t and y(t) NB. u is: function to use NB. x is: step size NB. eg: 0.001 fyp nextY 0 1 nextY=: adverb define
tableau=. 1 0.5 0.5, x * u y
ks=. (x * [: u y + (* x&,))/\. tableau
({:y) + 6 %~ +/ 1 2 2 1 * ks
)</lang>
Use:
report_err report_whole fyp rk4 1 0 10 0.1
Mathematica
<lang Mathematica>DSolve[{y'[t] == t * Sqrt[y[t]], y[0] == 1, y'[0] == 0}, y[t], t]//Simplify -> {{y[t] -> 1/16 (4+t^2)^2}} (* Symbolic expression found : no delta from expected solution *) y[t] -> 1/16*(4+t^2)^2 /. t -> Range[10] ->y[{1,2,3,4,5,6,7,8,9,10}]->{25/16, 4, 169/16, 25, 841/16, 100, 2809/16, 289, 7225/16, 676}</lang>
Maxima
<lang maxima>/* Here is how to solve a differential equation */ 'diff(y, x) = x * sqrt(y); ode2(%, y, x); ic1(%, x = 0, y = 1); factor(solve(%, y)); /* [y = (x^2 + 4)^2 / 16] */
/* The Runge-Kutta solver is builtin */
load(dynamics)$ sol: rk(t * sqrt(y), y, 1, [t, 0, 10, 1.0])$ plot2d([discrete, sol])$
/* An implementation of RK4 for one equation */
rk4(f, x0, y0, x1, n) := block([h, x, y, vx, vy, k1, k2, k3, k4],
h: bfloat((x1 - x0) / (n - 1)),
x: x0,
y: y0,
vx: makelist(0, n + 1),
vy: makelist(0, n + 1),
vx[1]: x0,
vy[1]: y0,
for i from 1 thru n do (
k1: bfloat(h * f(x, y)),
k2: bfloat(h * f(x + h / 2, y + k1 / 2)),
k3: bfloat(h * f(x + h / 2, y + k2 / 2)),
k4: bfloat(h * f(x + h, y + k3)),
vy[i + 1]: y: y + (k1 + 2 * k2 + 2 * k3 + k4) / 6,
vx[i + 1]: x: x + h
),
[vx, vy]
)$
[x, y]: rk4(lambda([x, y], x * sqrt(y)), 0, 1, 10, 101)$
plot2d([discrete, x, y])$
s: map(lambda([x], (x^2 + 4)^2 / 16), x)$
for i from 1 step 10 thru 101 do print(x[i], " ", y[i], " ", y[i] - s[i]);</lang>
МК-61/52
ПП 38 П1 ПП 30 П2 ПП 35 П3 2 * ПП 30 ИП2 ИП3 + 2 * + ИП1 + 3 / ИП7 + П7 П8 С/П БП 00 ИП6 ИП5 + П6 <-> ИП7 + П8 ИП8 КвКор ИП6 * ИП5 * В/О
Input: 1/2 (h/2) - Р5, 1 (y0) - Р8 and Р7, 0 (t0) - Р6.
OCaml
<lang ocaml>let y' t y = t *. sqrt y let exact t = let u = 0.25*.t*.t +. 1.0 in u*.u
let rk4_step (y,t) h =
let k1 = h *. y' t y in let k2 = h *. y' (t +. 0.5*.h) (y +. 0.5*.k1) in let k3 = h *. y' (t +. 0.5*.h) (y +. 0.5*.k2) in let k4 = h *. y' (t +. h) (y +. k3) in (y +. (k1+.k4)/.6.0 +. (k2+.k3)/.3.0, t +. h)
let rec loop h n (y,t) =
if n mod 10 = 1 then Printf.printf "t = %f,\ty = %f,\terr = %g\n" t y (abs_float (y -. exact t)); if n < 102 then loop h (n+1) (rk4_step (y,t) h)
let _ = loop 0.1 1 (1.0, 0.0)</lang>
Output:
t = 0.000000, y = 1.000000, err = 0 t = 1.000000, y = 1.562500, err = 1.45722e-07 t = 2.000000, y = 3.999999, err = 9.19479e-07 t = 3.000000, y = 10.562497, err = 2.90956e-06 t = 4.000000, y = 24.999994, err = 6.23491e-06 t = 5.000000, y = 52.562489, err = 1.08197e-05 t = 6.000000, y = 99.999983, err = 1.65946e-05 t = 7.000000, y = 175.562476, err = 2.35177e-05 t = 8.000000, y = 288.999968, err = 3.15652e-05 t = 9.000000, y = 451.562459, err = 4.07232e-05 t = 10.000000, y = 675.999949, err = 5.09833e-05
Perl 6
<lang perl6>sub runge-kutta (&yp, @t, $y is copy, \𝛿t) {
$y, gather for @t -> \t {
my \𝛿y1 = 𝛿t * yp t, $y; my \𝛿y2 = 𝛿t * yp t + 𝛿t / 2, $y + 𝛿y1 / 2; my \𝛿y3 = 𝛿t * yp t + 𝛿t / 2, $y + 𝛿y2 / 2; my \𝛿y4 = 𝛿t * yp t + 𝛿t, $y + 𝛿y3; take $y = $y + (𝛿y1 + 2 * (𝛿y2 + 𝛿y3) + 𝛿y4) / 6;
}
}
sub yprime ($t, $y) { $t * sqrt $y }
constant 𝛿t = 0.10; constant n = 10.0;
constant @t = 0, 𝛿t ... n; my @y = runge-kutta(&yprime, @t, 1, 𝛿t);
for @t Z @y -> $t, $y {
say "y($t) = $y".fmt("%-30s"), "±", abs($y - (($t ** 2 + 4) ** 2 / 16))
if $t == $t.floor; }</lang>
- Output:
y(0) = 1 ±0 y(1) = 1.5624998542781079 ±1.4572189210859676E-07 y(2) = 3.9999990805207992 ±9.1947920077828371E-07 y(3) = 10.562497090437551 ±2.9095624487496252E-06 y(4) = 24.999993765090636 ±6.2349093639113562E-06 y(5) = 52.562489180302585 ±1.0819697415342944E-05 y(6) = 99.999983405403583 ±1.6594596417007779E-05 y(7) = 175.56247648227125 ±2.3517728749311573E-05 y(8) = 288.99996843479857 ±3.156520142510999E-05 y(9) = 451.56245927683966 ±4.07231603389846E-05 y(10) = 675.99994901670971 ±5.0983290293515893E-05
Python
<lang Python>
- RHS of the ODE to be solved
def yp(t,y): from math import sqrt return t*sqrt(y)
- The RK4 routine
- INPUT ARGUMENTS
- yp function object providing the RHS of the 1st order ODE
- tmin, tmax the domain of the solution
- y0, t0 initial values
- dt time step over which the solution is sought
- OUTPUT
- ti the list of points where the equation has been solved
- yi the values of the solution found
def RK4(yp,tmin,tmax,y0,t0,dt): yi=[y0] ti=[t0] nmax=int( (tmax-tmin)/dt +1) for n in range(1,nmax,1): tn=ti[n-1] yn=yi[n-1] dy1=dt*yp( tn, yn) dy2=dt*yp( tn+dt/2.0, yn+dy1/2.0) dy3=dt*yp( tn+dt/2.0, yn+dy2/2.0) dy4=dt*yp( tn+dt, yn+dy3) yi.append( yn+(1.0/6.0)*(dy1+2.0*dy2+2.0*dy3+dy4) ) ti.append(tn+dt) return [ti,yi]
[tmin, tmax, dt] = [0.0, 10.0, 0.1] [y0, t0] = [1.0, 0.0] [t,y]=RK4(yp,tmin,tmax,y0,t0,dt) for i in range(0,len(t),10): print ("y(%2.1f)\t= %4.6f \t error:%4.6g")%(t[i], y[i], abs( y[i]- ((t[i]**2 + 4.0)**2)/16.0 ) )
</lang> Output
y(0.0) = 1.000000 error: 0 y(1.0) = 1.562500 error:1.45722e-07 y(2.0) = 3.999999 error:9.19479e-07 y(3.0) = 10.562497 error:2.90956e-06 y(4.0) = 24.999994 error:6.23491e-06 y(5.0) = 52.562489 error:1.08197e-05 y(6.0) = 99.999983 error:1.65946e-05 y(7.0) = 175.562476 error:2.35177e-05 y(8.0) = 288.999968 error:3.15652e-05 y(9.0) = 451.562459 error:4.07232e-05 y(10.0) = 675.999949 error:5.09833e-05
REXX
<lang rexx>/*REXX program uses the Runge-Kutta method to solve the differential */ /* ____ */ /*equation: y'(t)=t²√y(t) which has the exact solution: y(t)=(t²+4)²/16*/
numeric digits 40; d=digits()%2 /*use forty digits, show ½ that. */ x0=0; x1=10; dx=.1; n=1 + (x1-x0) / dx; y.=1
do m=1 for n-1; mm=m-1
y.m=Runge_Kutta(dx, x0+dx*mm, y.mm)
end /*m*/
say center(x,13,'─') center(y,d,'─') ' ' center('relative error',d,'─')
do i=0 to n-1 by 10; x=(x0+dx*i)/1; y2=(x*x/4+1)**2
relE=format(y.i/y2-1,,13)/1; if relE=0 then relE=' 0'
say center(x,13) right(format(y.i,,12),d) ' ' left(relE,d)
end /*i*/
exit /*stick a fork in it, we're done.*/ /*──────────────────────────────────RATE subroutine─────────────────────*/ rate: return arg(1)*sqrt(arg(2)) /*──────────────────────────────────Runge_Kutta subroutine──────────────*/ Runge_Kutta: procedure; parse arg dx,x,y
k1 = dx * rate(x , y )
k2 = dx * rate(x+dx/2 , y+k1/2 )
k3 = dx * rate(x+dx/2 , y+k2/2 )
k4 = dx * rate(x+dx , y+k3 )
return y + (k1 + 2*k2 + 2*k3 + k4) / 6 /*──────────────────────────────────SQRT subroutine─────────────────────*/ sqrt: procedure;parse arg x; if x=0 then return 0; d=digits()
numeric digits 11; g=.sqrtG()
do j=0 while p>9; m.j=p; p=p%2+1; end; do k=j+5 to 0 by -1
if m.k>11 then numeric digits m.k
g=.5*(g+x/g); end; numeric digits d; return g/1
.sqrtG: numeric form; m.=11; p=d+d%4+2
parse value format(x,2,1,,0) 'E0' with g 'E' _ .; return g*.5'E'_%2</lang>
output
──────X────── ─────────Y────────── ───relative error───
0 1.000000000000 0
1 1.562499854278 -9.3262010934636E-8
2 3.999999080521 -2.2986980018794E-7
3 10.562497090438 -2.7546153355698E-7
4 24.999993765091 -2.4939637458826E-7
5 52.562489180303 -2.058444217357E-7
6 99.999983405404 -1.6594596403363E-7
7 175.562476482271 -1.3395644712797E-7
8 288.999968434799 -1.092221504E-7
9 451.562459276840 -9.018277747572E-8
10 675.999949016709 -7.5419068845528E-8
Tcl
<lang tcl>package require Tcl 8.5
- Hack to bring argument function into expression
proc tcl::mathfunc::dy {t y} {upvar 1 dyFn dyFn; $dyFn $t $y}
proc rk4step {dyFn y* t* dt} {
upvar 1 ${y*} y ${t*} t
set dy1 [expr {$dt * dy($t, $y)}]
set dy2 [expr {$dt * dy($t+$dt/2, $y+$dy1/2)}]
set dy3 [expr {$dt * dy($t+$dt/2, $y+$dy2/2)}]
set dy4 [expr {$dt * dy($t+$dt, $y+$dy3)}]
set y [expr {$y + ($dy1 + 2*$dy2 + 2*$dy3 + $dy4)/6.0}]
set t [expr {$t + $dt}]
}
proc y {t} {expr {($t**2 + 4)**2 / 16}} proc δy {t y} {expr {$t * sqrt($y)}}
proc printvals {t y} {
set err [expr {abs($y - [y $t])}]
puts [format "y(%.1f) = %.8f\tError: %.8e" $t $y $err]
}
set t 0.0 set y 1.0 set dt 0.1 printvals $t $y for {set i 1} {$i <= 101} {incr i} {
rk4step δy y t $dt
if {$i%10 == 0} {
printvals $t $y
}
}</lang>
- Output:
y(0.0) = 1.00000000 Error: 0.00000000e+00 y(1.0) = 1.56249985 Error: 1.45721892e-07 y(2.0) = 3.99999908 Error: 9.19479203e-07 y(3.0) = 10.56249709 Error: 2.90956245e-06 y(4.0) = 24.99999377 Error: 6.23490939e-06 y(5.0) = 52.56248918 Error: 1.08196973e-05 y(6.0) = 99.99998341 Error: 1.65945961e-05 y(7.0) = 175.56247648 Error: 2.35177280e-05 y(8.0) = 288.99996843 Error: 3.15652000e-05 y(9.0) = 451.56245928 Error: 4.07231581e-05 y(10.0) = 675.99994902 Error: 5.09832864e-05