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

Numeric.RungeKutta

Contents

Forward Euler
Kutta 3rd order
Kutta 4th order
Kutta 4th order, more precise
Heun-Euler, order 2/1
Bogacki-Shampine, order 3/2
Runge-Kutta-Fehlberg, order 4/5
Cash-Karp, order 4/5
Dormand-Prince, order 5/4
Fehlberg, order 7/8
Verner, order 6/5

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

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

rkfe Source #

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

Kutta 3rd order

rk3 Source #

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

Kutta 4th order

rk4a Source #

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

Kutta 4th order, more precise

rk4b Source #

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

Heun-Euler, order 2/1

rkhe Source #

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

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

show_rkhe_aux :: String Source #

Bogacki-Shampine, order 3/2

rkbs Source #

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

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

show_rkbs_aux :: String Source #

Runge-Kutta-Fehlberg, order 4/5

rkf45 Source #

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

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

show_rkf45_aux :: String Source #

Cash-Karp, order 4/5

rkck Source #

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

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

show_rkck_aux :: String Source #

Dormand-Prince, order 5/4

rkdp Source #

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

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

show_rkdp_aux :: String Source #

Fehlberg, order 7/8

rkf78 Source #

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

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

show_rkf78_aux :: String Source #

Verner, order 6/5

rkv65 Source #

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

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

show_rkv65_aux :: String Source #