Numeric.RungeKutta

Given an ODE of the form

\[ \partial_t y = f(t, y) \]

Runge-Kutta integrators approximate the evolution in time of the system with a summation

\[ y_{n+1} = y_n + h \sum_{i=1}^s b_i k_i \]

where

\[ k_i = f(t_n + c_i h, y_n + h \sum_{j=1}^s a_{ij} k_j) \]

And the \(a_{i, j}\) coefficients are taken from a Butcher Tableau of the form

c_1	a_11	a_12	...	a_1s
c_2	a_21	a_22	...	a_2s
...	...
c_s	a_s1	a_s2	...	a_ss
	b_1	b_2	...	b_s

This module implements a method that can use a generic tableau, then specializes with different tableaux to let the user pick a specific method. Adaptive step-size methods, see below, add a row of \(d_j\) coefficients and use that to report the error:

\[ e_{n+1} = h \sum_{i=1}^s d_i k_i \]

Non-adaptive solvers: rkfe, rk3, rk4a, rk4b

Adaptive solvers: rkhe, rkbs, rkf45, rkck, rkdp, rkf78, rkv65

Auxiliary non-adaptive solvers (error estimators from the adaptive ones): rkhe_aux, rkbs_aux, rkf45_aux, rkck_aux, rkdp_aux, rkf78_aux, rkv65_aux

Use rk4a or rk4b if you don't need an adaptive solver, rkdp or rkv65 if you do;

N.B. : DO NOT USE rkfe EXCEPT FOR DEMONSTRATIONS OF INSTABILITY!

Reference: E. Hairer, S. P. Norsett, G. Wanner, Solving Ordinary Differential Equations I: Nonstiff Problems (second revised edition, 1993).

Forward Euler

rkfe

Arguments

:: (Double -> a -> a)

scale function to scale a Y state vector

-> (a -> a -> a)

sum function to add two Y state vectors

-> (Double -> a -> a)

derivative function F

-> Double

step size H

-> (Double, a)

current state (t, Y)

-> (Double, a)

new state (t_new, Y_new)

forward Euler: unconditionally unstable: don't use this! If you do, your dangly bits will fall off!

show_rkfe :: String

Kutta 3rd order

rk3

Arguments

:: (Double -> a -> a)

scale function to scale a Y state vector

-> (a -> a -> a)

sum function to add two Y state vectors

-> (Double -> a -> a)

derivative function F

-> Double

step size H

-> (Double, a)

current state (t, Y)

-> (Double, a)

new state (t_new, Y_new)

Kutta's third-order method

show_rk3 :: String

Kutta 4th order

rk4a

Arguments

:: (Double -> a -> a)

scale function to scale a Y state vector

-> (a -> a -> a)

sum function to add two Y state vectors

-> (Double -> a -> a)

derivative function F

-> Double

step size H

-> (Double, a)

current state (t, Y)

-> (Double, a)

new state (t_new, Y_new)

Classic fourth-order method

show_rk4a :: String

Kutta 4th order, more precise

rk4b

Arguments

:: (Double -> a -> a)

scale function to scale a Y state vector

-> (a -> a -> a)

sum function to add two Y state vectors

-> (Double -> a -> a)

derivative function F

-> Double

step size H

-> (Double, a)

current state (t, Y)

-> (Double, a)

new state (t_new, Y_new)

Kutta's other fourth-order method... "The first [above] is more popular, the second is more precise." (Hairer, Norsett, Wanner)

show_rk4b :: String

Heun-Euler, order 2/1

rkhe

Arguments

:: (Double -> a -> a)

scale function to scale a Y state vector

-> (a -> a -> a)

sum function to add two Y state vectors

-> (Double -> a -> a)

derivative function F

-> Double

step size H

-> (Double, a)

current state (t, Y)

-> (Double, a, a)

new state (t_new, Y_new, e_new)

Heun-Euler, order 2/1

show_rkhe :: String

rkhe_aux :: (Double -> a -> a) -> (a -> a -> a) -> (Double -> a -> a) -> Double -> (Double, a) -> (Double, a)

show_rkhe_aux :: String

Bogacki-Shampine, order 3/2

rkbs

Arguments

:: (Double -> a -> a)

scale function to scale a Y state vector

-> (a -> a -> a)

sum function to add two Y state vectors

-> (Double -> a -> a)

derivative function F

-> Double

step size H

-> (Double, a)

current state (t, Y)

-> (Double, a, a)

new state (t_new, Y_new, e_new)

Bogacki-Shampine, order 3/2

show_rkbs :: String

rkbs_aux :: (Double -> a -> a) -> (a -> a -> a) -> (Double -> a -> a) -> Double -> (Double, a) -> (Double, a)

show_rkbs_aux :: String

Runge-Kutta-Fehlberg, order 4/5

rkf45

Arguments

:: (Double -> a -> a)

scale function to scale a Y state vector

-> (a -> a -> a)

sum function to add two Y state vectors

-> (Double -> a -> a)

derivative function F

-> Double

step size H

-> (Double, a)

current state (t, Y)

-> (Double, a, a)

new state (t_new, Y_new, e_new)

Runge-Kutta-Fehlberg, order 4/5

show_rkf45 :: String

rkf45_aux :: (Double -> a -> a) -> (a -> a -> a) -> (Double -> a -> a) -> Double -> (Double, a) -> (Double, a)

show_rkf45_aux :: String

Cash-Karp, order 4/5

rkck

Arguments

:: (Double -> a -> a)

scale function to scale a Y state vector

-> (a -> a -> a)

sum function to add two Y state vectors

-> (Double -> a -> a)

derivative function F

-> Double

step size H

-> (Double, a)

current state (t, Y)

-> (Double, a, a)

new state (t_new, Y_new, e_new)

Cash-Karp, order 4/5 (use 4th-order sol'n, coeffs chosen to minimize error of 4th-order sol'n)

show_rkck :: String

rkck_aux :: (Double -> a -> a) -> (a -> a -> a) -> (Double -> a -> a) -> Double -> (Double, a) -> (Double, a)

show_rkck_aux :: String

Dormand-Prince, order 5/4

rkdp

Arguments

:: (Double -> a -> a)

scale function to scale a Y state vector

-> (a -> a -> a)

sum function to add two Y state vectors

-> (Double -> a -> a)

derivative function F

-> Double

step size H

-> (Double, a)

current state (t, Y)

-> (Double, a, a)

new state (t_new, Y_new, e_new)

Dormand-Prince, order 5/4 (use 5th-order sol'n, coeffs chosen to minimize error of 5th-order sol'n) This is DOPRI5 from Hairer, Norsett, Wanner

show_rkdp :: String

rkdp_aux :: (Double -> a -> a) -> (a -> a -> a) -> (Double -> a -> a) -> Double -> (Double, a) -> (Double, a)

show_rkdp_aux :: String

Fehlberg, order 7/8

rkf78

Arguments

:: (Double -> a -> a)

scale function to scale a Y state vector

-> (a -> a -> a)

sum function to add two Y state vectors

-> (Double -> a -> a)

derivative function F

-> Double

step size H

-> (Double, a)

current state (t, Y)

-> (Double, a, a)

new state (t_new, Y_new, e_new)

Fehlberg, order 7/8: "... of frequent use in all high precision computations, e.g., in astronomy." -- Hairer, Norsett, Wanner. But caveat: suffers from the drawback that error estimate is identically 0 for quadrature problems y' = f(x)

NOTE BUG in Hairer Norsett Wanner: the third-last A coefficient in the last row of the tableau is listed as "19/41" in the book. This is WRONG: that row does not then sum to 1, and the convergence of the auxiliary solver is then order 1 or 2, not 8

show_rkf78 :: String

rkf78_aux :: (Double -> a -> a) -> (a -> a -> a) -> (Double -> a -> a) -> Double -> (Double, a) -> (Double, a)

show_rkf78_aux :: String

Verner, order 6/5

rkv65

Arguments

:: (Double -> a -> a)

scale function to scale a Y state vector

-> (a -> a -> a)

sum function to add two Y state vectors

-> (Double -> a -> a)

derivative function F

-> Double

step size H

-> (Double, a)

current state (t, Y)

-> (Double, a, a)

new state (t_new, Y_new, e_new)

Verner, order 6/5 (DVERK)

show_rkv65 :: String

rkv65_aux :: (Double -> a -> a) -> (a -> a -> a) -> (Double -> a -> a) -> Double -> (Double, a) -> (Double, a)

show_rkv65_aux :: String