Copyright | (c) Uwe Hollerbach uh@alumni.caltech.edu 2009 Marco Zocca 2022 |
---|---|

License | BSD-3 |

Maintainer | ocramz |

Stability | experimental |

Portability | POSIX |

Safe Haskell | Safe-Inferred |

Language | Haskell2010 |

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).

## Synopsis

- rkfe :: (Double -> a -> a) -> (a -> a -> a) -> (Double -> a -> a) -> Double -> (Double, a) -> (Double, a)
- show_rkfe :: String
- rk3 :: (Double -> a -> a) -> (a -> a -> a) -> (Double -> a -> a) -> Double -> (Double, a) -> (Double, a)
- show_rk3 :: String
- rk4a :: (Double -> a -> a) -> (a -> a -> a) -> (Double -> a -> a) -> Double -> (Double, a) -> (Double, a)
- show_rk4a :: String
- rk4b :: (Double -> a -> a) -> (a -> a -> a) -> (Double -> a -> a) -> Double -> (Double, a) -> (Double, a)
- show_rk4b :: String
- rkhe :: (Double -> a -> a) -> (a -> a -> a) -> (Double -> a -> a) -> Double -> (Double, a) -> (Double, a, a)
- show_rkhe :: String
- rkhe_aux :: (Double -> a -> a) -> (a -> a -> a) -> (Double -> a -> a) -> Double -> (Double, a) -> (Double, a)
- show_rkhe_aux :: String
- rkbs :: (Double -> a -> a) -> (a -> a -> a) -> (Double -> a -> a) -> Double -> (Double, a) -> (Double, a, a)
- show_rkbs :: String
- rkbs_aux :: (Double -> a -> a) -> (a -> a -> a) -> (Double -> a -> a) -> Double -> (Double, a) -> (Double, a)
- show_rkbs_aux :: String
- rkf45 :: (Double -> a -> a) -> (a -> a -> a) -> (Double -> a -> a) -> Double -> (Double, a) -> (Double, a, a)
- show_rkf45 :: String
- rkf45_aux :: (Double -> a -> a) -> (a -> a -> a) -> (Double -> a -> a) -> Double -> (Double, a) -> (Double, a)
- show_rkf45_aux :: String
- rkck :: (Double -> a -> a) -> (a -> a -> a) -> (Double -> a -> a) -> Double -> (Double, a) -> (Double, a, a)
- show_rkck :: String
- rkck_aux :: (Double -> a -> a) -> (a -> a -> a) -> (Double -> a -> a) -> Double -> (Double, a) -> (Double, a)
- show_rkck_aux :: String
- rkdp :: (Double -> a -> a) -> (a -> a -> a) -> (Double -> a -> a) -> Double -> (Double, a) -> (Double, a, a)
- show_rkdp :: String
- rkdp_aux :: (Double -> a -> a) -> (a -> a -> a) -> (Double -> a -> a) -> Double -> (Double, a) -> (Double, a)
- show_rkdp_aux :: String
- rkf78 :: (Double -> a -> a) -> (a -> a -> a) -> (Double -> a -> a) -> Double -> (Double, a) -> (Double, a, a)
- show_rkf78 :: String
- rkf78_aux :: (Double -> a -> a) -> (a -> a -> a) -> (Double -> a -> a) -> Double -> (Double, a) -> (Double, a)
- show_rkf78_aux :: String
- rkv65 :: (Double -> a -> a) -> (a -> a -> a) -> (Double -> a -> a) -> Double -> (Double, a) -> (Double, a, a)
- show_rkv65 :: String
- rkv65_aux :: (Double -> a -> a) -> (a -> a -> a) -> (Double -> a -> a) -> Double -> (Double, a) -> (Double, a)
- show_rkv65_aux :: String

# Forward Euler

:: (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!

# Kutta 3rd order

:: (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

# Kutta 4th order

:: (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

# Kutta 4th order, more precise

:: (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)

# Heun-Euler, order 2/1

:: (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

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

# Bogacki-Shampine, order 3/2

:: (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

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

# Runge-Kutta-Fehlberg, order 4/5

:: (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 Source #

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

# Cash-Karp, order 4/5

:: (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)

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

# Dormand-Prince, order 5/4

:: (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

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

# Fehlberg, order 7/8

:: (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 Source #

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

# Verner, order 6/5

:: (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 Source #