Runge-Kutta method: Difference between revisions
Content deleted Content added
m →{{header|Tcl}}: formatting glitch |
m Better integration of WP link |
||
Line 7: | Line 7: | ||
:<math>y(t) = \tfrac{1}{16}(t^2 +4)^2</math> |
:<math>y(t) = \tfrac{1}{16}(t^2 +4)^2</math> |
||
;Task |
;Task |
||
Demonstrate the commonly used explicit |
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) |
* 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. |
* 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. |
Revision as of 15:37, 15 March 2012
Runge-Kutta method is a draft programming task. It is not yet considered ready to be promoted as a complete task, for reasons that should be found in its talk page.
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
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