Runge-Kutta method: Difference between revisions
Content deleted Content added
m moved Runge-Kutta to Runge-Kutta method: Page names should be clearer |
m Improve formatting of problem (\delta t not dt, \times not *) |
||
Line 1: | Line 1: | ||
{{draft task}} |
{{draft task}} |
||
Given the example Differential equation: |
Given the example Differential equation: |
||
:<math>y' = t |
:<math>y' = t \times \sqrt y</math> |
||
With initial condition: |
With initial condition: |
||
:<math>t_0 = 0</math> and <math>y_0 = y(t_0) = y(0) = 1</math> |
:<math>t_0 = 0</math> and <math>y_0 = y(t_0) = y(0) = 1</math> |
||
This equation has an exact solution: |
This equation has an exact solution: |
||
:<math>y = \tfrac{1}{16}(t^2 +4)^2</math> |
:<math>y = \tfrac{1}{16}(t^2 +4)^2</math> |
||
;Task |
|||
Demonstrate the commonly used explicit fourth-order Runge–Kutta method as defined in the [[wp:Runge–Kutta_methods#Common_fourth-order_Runge.E2.80.93Kutta_method|Wikipedia article]] 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> |
:<math>\delta y_1 = \delta t\times y'(t_n, y_n)\quad</math> |
||
:<math> |
:<math>\delta y_2 = \delta t\times y'(t_n + \tfrac{1}{2}\delta t , y_n + \tfrac{1}{2}\delta y_1)</math> |
||
:<math> |
:<math>\delta y_3 = \delta t\times y'(t_n + \tfrac{1}{2}\delta t , y_n + \tfrac{1}{2}\delta y_2)</math> |
||
:<math> |
:<math>\delta y_4 = \delta t\times y'(t_n + \delta t , y_n + \delta y_3)\quad</math> |
||
then: |
|||
:<math>y_{n+1} = y_n + \tfrac{1}{6} ( |
:<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 + |
:<math>t_{n+1} = t_n + \delta t\quad</math> |
||
* The reference implementation is provided in Ada. |
|||
=={{header|Ada}}== |
=={{header|Ada}}== |
Revision as of 10:11, 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 as defined in the Wikipedia article 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