Runge-Kutta method: Difference between revisions
+ D entry |
Stronger static typing in the D entry |
||
Line 88: | Line 88: | ||
=={{header|D}}== |
=={{header|D}}== |
||
{{trans|Ada}} |
{{trans|Ada}} |
||
<lang d>import std.stdio, std.math; |
<lang d>import std.stdio, std.math, std.typecons; |
||
alias real FP; |
alias real FP; |
||
alias Typedef!(FP[101]) FPs; |
|||
void runge(in FP function(in FP, in FP) pure nothrow yp_func, |
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) { |
foreach (n; 0 .. t.length - 1) { |
||
FP dy1 = dt * yp_func(t[n], y[n]); |
FP dy1 = dt * yp_func(t[n], y[n]); |
||
Line 111: | Line 112: | ||
void main() { |
void main() { |
||
enum FP dt = 0.10; |
enum FP dt = 0.10; |
||
FPs t_arr, y_arr; |
|||
FP[N] t_arr, y_arr; |
|||
t_arr[0] = 0.0; |
t_arr[0] = 0.0; |
Revision as of 18:28, 15 March 2012
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
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
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