rungekutta2-1.0.3: Explicit Runge-Kutta methods of various orders (fork of 'rungekutta')
Copyright (c) Uwe Hollerbach uh@alumni.caltech.edu 2009Marco Zocca 2022 BSD-3 ocramz experimental POSIX Safe-Inferred Haskell2010

Numeric.RungeKutta

Description

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

# Forward Euler

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!

# Kutta 3rd order

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

# Kutta 4th order

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

# Kutta 4th order, more precise

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)

# Heun-Euler, order 2/1

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

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

# Bogacki-Shampine, order 3/2

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

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

# Runge-Kutta-Fehlberg, order 4/5

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

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

# Cash-Karp, order 4/5

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)

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

# Dormand-Prince, order 5/4

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

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

# Fehlberg, order 7/8

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

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

# Verner, order 6/5

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)

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