{-|
Module      : Numeric.RungeKutta
Description : Runge-Kutta integration of ODEs
Copyright   : (c) Uwe Hollerbach, uh@alumni.caltech.edu, 2009
                  Marco Zocca, 2022
License     : BSD-3
Maintainer  : ocramz
Stability   : experimental
Portability : POSIX


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).
-}
module Numeric.RungeKutta (
  -- * Forward Euler
  rkfe, show_rkfe,
  -- * Kutta 3rd order
  rk3, show_rk3,
  -- * Kutta 4th order
  rk4a, show_rk4a,
  -- * Kutta 4th order, more precise
  rk4b, show_rk4b,
  -- * Heun-Euler, order 2/1
  rkhe, show_rkhe, rkhe_aux, show_rkhe_aux,
  -- * Bogacki-Shampine, order 3/2
  rkbs, show_rkbs, rkbs_aux, show_rkbs_aux,
  -- * Runge-Kutta-Fehlberg, order 4/5
  rkf45, show_rkf45, rkf45_aux, show_rkf45_aux,
  -- * Cash-Karp, order 4/5
  rkck, show_rkck, rkck_aux, show_rkck_aux,
  -- * Dormand-Prince, order 5/4
  rkdp, show_rkdp, rkdp_aux, show_rkdp_aux,
  -- * Fehlberg, order 7/8
  rkf78, show_rkf78, rkf78_aux, show_rkf78_aux,
  -- * Verner, order 6/5
  rkv65, show_rkv65, rkv65_aux, show_rkv65_aux)
  where

import Data.Ratio
import Data.Maybe
import Data.List


-- Store a list of rational numbers as a common denominator, then a list
-- of numerators, all stored as doubles: this lets us write the values in
-- the source as exact rational numbers, yet we can take advantage of not
-- having to do too many conversions and also of not having to multiply
-- by 0 or 1 in some cases.

type RCL = (Double, [Double])

ratToRCL :: [Rational] -> RCL
ratToRCL :: [Rational] -> RCL
ratToRCL [] = (Double
1.0, [])
ratToRCL [Rational]
rs =
  let ds :: [Integer]
ds = (Rational -> Integer) -> [Rational] -> [Integer]
forall a b. (a -> b) -> [a] -> [b]
map Rational -> Integer
forall a. Ratio a -> a
denominator [Rational]
rs
      dp :: Integer
dp = (Integer -> Integer -> Integer) -> [Integer] -> Integer
forall a. (a -> a -> a) -> [a] -> a
foldl1' Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
(*) [Integer]
ds
      ns :: [Integer]
ns = (Rational -> Integer) -> [Rational] -> [Integer]
forall a b. (a -> b) -> [a] -> [b]
map (Rational -> Integer
forall a. Ratio a -> a
numerator (Rational -> Integer)
-> (Rational -> Rational) -> Rational -> Integer
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (Rational -> Rational -> Rational
forall a. Num a => a -> a -> a
*(Integer -> Rational
forall a. Num a => Integer -> a
fromInteger Integer
dp))) [Rational]
rs
      g :: Integer
g = (Integer -> Integer -> Integer) -> Integer -> [Integer] -> Integer
forall (t :: * -> *) b a.
Foldable t =>
(b -> a -> b) -> b -> t a -> b
foldl' Integer -> Integer -> Integer
forall a. Integral a => a -> a -> a
gcd Integer
dp [Integer]
ns
  in (Integer -> Double
forall a. Num a => Integer -> a
fromInteger (Integer -> Integer -> Integer
forall a. Integral a => a -> a -> a
quot Integer
dp Integer
g), (Integer -> Double) -> [Integer] -> [Double]
forall a b. (a -> b) -> [a] -> [b]
map (Integer -> Double
forall a. Num a => Integer -> a
fromInteger (Integer -> Double) -> (Integer -> Integer) -> Integer -> Double
forall b c a. (b -> c) -> (a -> b) -> a -> c
. ((Integer -> Integer -> Integer) -> Integer -> Integer -> Integer
forall a b c. (a -> b -> c) -> b -> a -> c
flip Integer -> Integer -> Integer
forall a. Integral a => a -> a -> a
quot Integer
g)) [Integer]
ns)

ratToRCLs :: [[Rational]] -> [RCL]
ratToRCLs :: [[Rational]] -> [RCL]
ratToRCLs = ([Rational] -> RCL) -> [[Rational]] -> [RCL]
forall a b. (a -> b) -> [a] -> [b]
map [Rational] -> RCL
ratToRCL

ratToDbls :: [Rational] -> [Double]
ratToDbls :: [Rational] -> [Double]
ratToDbls = (Rational -> Double) -> [Rational] -> [Double]
forall a b. (a -> b) -> [a] -> [b]
map Rational -> Double
forall a. Fractional a => Rational -> a
fromRational

-- Helper function to sum a list of K_i, skipping
-- un-necessary multiplications and additions

k_sum :: (t -> t -> t) -> (t -> t -> t) -> t -> (t, [t]) -> [t] -> t
k_sum t -> t -> t
sc_fn t -> t -> t
sum_fn t
h (t
d,[t]
ns) [t]
ks =
  t -> t -> t
sc_fn (t
ht -> t -> t
forall a. Fractional a => a -> a -> a
/t
d) ((t -> t -> t) -> [t] -> t
forall (t :: * -> *) a. Foldable t => (a -> a -> a) -> t a -> a
foldl1 t -> t -> t
sum_fn ([Maybe t] -> [t]
forall a. [Maybe a] -> [a]
catMaybes ((t -> t -> Maybe t) -> [t] -> [t] -> [Maybe t]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith t -> t -> Maybe t
m_scale [t]
ns [t]
ks)))
  where m_scale :: t -> t -> Maybe t
m_scale t
s t
v | t
s t -> t -> Bool
forall a. Eq a => a -> a -> Bool
== t
0.0    = Maybe t
forall a. Maybe a
Nothing
                    | t
s t -> t -> Bool
forall a. Eq a => a -> a -> Bool
== t
1.0    = t -> Maybe t
forall a. a -> Maybe a
Just t
v
                    | Bool
otherwise   = t -> Maybe t
forall a. a -> Maybe a
Just (t -> t -> t
sc_fn t
s t
v)

-- Helper function to generate a list of K_i

gen_ks :: (t -> [a] -> t)
-> (t -> t -> t)
-> (t -> t -> a)
-> t
-> (t, t)
-> [a]
-> [t]
-> [t]
-> [a]
gen_ks t -> [a] -> t
ksum_fn t -> t -> t
sum_fn t -> t -> a
der_fn t
h (t
tn,t
yn) [a]
ks [] [] = [a]
ks
gen_ks t -> [a] -> t
ksum_fn t -> t -> t
sum_fn t -> t -> a
der_fn t
h old :: (t, t)
old@(t
tn,t
yn) [a]
ks (t
c:[t]
cs) (t
a:[t]
as) =
  (t -> [a] -> t)
-> (t -> t -> t)
-> (t -> t -> a)
-> t
-> (t, t)
-> [a]
-> [t]
-> [t]
-> [a]
gen_ks t -> [a] -> t
ksum_fn t -> t -> t
sum_fn t -> t -> a
der_fn t
h (t, t)
old ([a]
ks [a] -> [a] -> [a]
forall a. [a] -> [a] -> [a]
++ [t -> t -> a
der_fn (t
tn t -> t -> t
forall a. Num a => a -> a -> a
+ t
ct -> t -> t
forall a. Num a => a -> a -> a
*t
h)
          (if [a] -> Bool
forall (t :: * -> *) a. Foldable t => t a -> Bool
null [a]
ks then t
yn else t -> t -> t
sum_fn t
yn (t -> [a] -> t
ksum_fn t
a [a]
ks))]) [t]
cs [t]
as

-- | This is the first core routine: it does not get exported,
-- only partial applications of it; see below.
--
-- Its arguments are:
--   c table (specified internally)
--   a table (specified internally)
--   b table (specified internally)
-- (user-specified arguments from here on:)
--   scale function to scale a Y state vector ::
--	(Double -> a -> a)
--   sum function to add two Y state vectors ::
--	(a -> a -> a)
--   derivative function F ::
--	(Double -> a -> a)
--   step size H ::
--	Double
--   current state (T,Y) ::
--	(Double, a)
--   and the return value is the new state (T_new,Y_new)
core1 ::
  [Double] -> [RCL] -> RCL -> (Double -> a -> a) ->
  (a -> a -> a) -> (Double -> a -> a) -> Double -> (Double, a) -> (Double, a)
core1 :: [Double]
-> [RCL]
-> RCL
-> (Double -> a -> a)
-> (a -> a -> a)
-> (Double -> a -> a)
-> Double
-> (Double, a)
-> (Double, a)
core1 [Double]
cl [RCL]
al RCL
bl Double -> a -> a
sc_fn a -> a -> a
sum_fn Double -> a -> a
der_fn Double
h old :: (Double, a)
old@(Double
tn,a
yn) =
  let ksum :: RCL -> [a] -> a
ksum = (Double -> a -> a) -> (a -> a -> a) -> Double -> RCL -> [a] -> a
forall t t.
(Eq t, Fractional t) =>
(t -> t -> t) -> (t -> t -> t) -> t -> (t, [t]) -> [t] -> t
k_sum Double -> a -> a
sc_fn a -> a -> a
sum_fn Double
h
  in (Double
tn Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
h, a -> a -> a
sum_fn a
yn (RCL -> [a] -> a
ksum RCL
bl ((RCL -> [a] -> a)
-> (a -> a -> a)
-> (Double -> a -> a)
-> Double
-> (Double, a)
-> [a]
-> [Double]
-> [RCL]
-> [a]
forall t t a t t.
Num t =>
(t -> [a] -> t)
-> (t -> t -> t)
-> (t -> t -> a)
-> t
-> (t, t)
-> [a]
-> [t]
-> [t]
-> [a]
gen_ks RCL -> [a] -> a
ksum a -> a -> a
sum_fn Double -> a -> a
der_fn Double
h (Double, a)
old [] [Double]
cl [RCL]
al)))

-- | This is the second core routine, analogous to the previous one.
-- The difference is that this gets an additional internal table arg,
-- and it returns a 3-tuple instead of a 2-tuple: (tnew,ynew,enew),
-- where enew is the error state vector
-- e_{n+1} = h sum_{i=1}^s (b_i - b'_i) k_i
core2 ::
  [Double] -> [RCL] -> RCL -> RCL -> (Double -> a -> a) ->
  (a -> a -> a) -> (Double -> a -> a) -> Double -> (Double, a) ->
  (Double, a, a)
core2 :: [Double]
-> [RCL]
-> RCL
-> RCL
-> (Double -> a -> a)
-> (a -> a -> a)
-> (Double -> a -> a)
-> Double
-> (Double, a)
-> (Double, a, a)
core2 [Double]
cl [RCL]
al RCL
bl RCL
dl Double -> a -> a
sc_fn a -> a -> a
sum_fn Double -> a -> a
der_fn Double
h old :: (Double, a)
old@(Double
tn,a
yn) =
  let ksum :: RCL -> [a] -> a
ksum = (Double -> a -> a) -> (a -> a -> a) -> Double -> RCL -> [a] -> a
forall t t.
(Eq t, Fractional t) =>
(t -> t -> t) -> (t -> t -> t) -> t -> (t, [t]) -> [t] -> t
k_sum Double -> a -> a
sc_fn a -> a -> a
sum_fn Double
h
      ks :: [a]
ks = (RCL -> [a] -> a)
-> (a -> a -> a)
-> (Double -> a -> a)
-> Double
-> (Double, a)
-> [a]
-> [Double]
-> [RCL]
-> [a]
forall t t a t t.
Num t =>
(t -> [a] -> t)
-> (t -> t -> t)
-> (t -> t -> a)
-> t
-> (t, t)
-> [a]
-> [t]
-> [t]
-> [a]
gen_ks RCL -> [a] -> a
ksum a -> a -> a
sum_fn Double -> a -> a
der_fn Double
h (Double, a)
old [] [Double]
cl [RCL]
al
  in (Double
tn Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
h, a -> a -> a
sum_fn a
yn (RCL -> [a] -> a
ksum RCL
bl [a]
ks), RCL -> [a] -> a
ksum RCL
dl [a]
ks)

-- A pair of helper routines to show the internal tables

rk_show1 :: String -> [Double] -> [RCL] -> RCL -> String
rk_show1 :: String -> [Double] -> [RCL] -> RCL -> String
rk_show1 String
title [Double]
cs [RCL]
as RCL
bs =
  String
title String -> String -> String
forall a. [a] -> [a] -> [a]
++ String
":\ncs:\t" String -> String -> String
forall a. [a] -> [a] -> [a]
++ ([Double] -> String
forall a. Show a => a -> String
show [Double]
cs) String -> String -> String
forall a. [a] -> [a] -> [a]
++ String
"\nas:\t" String -> String -> String
forall a. [a] -> [a] -> [a]
++
    (String -> String -> String) -> [String] -> String
forall a. (a -> a -> a) -> [a] -> a
foldl1' String -> String -> String
gl ((RCL -> String) -> [RCL] -> [String]
forall a b. (a -> b) -> [a] -> [b]
map RCL -> String
forall a. Show a => a -> String
show [RCL]
as) String -> String -> String
forall a. [a] -> [a] -> [a]
++ String
"\nbs:\t" String -> String -> String
forall a. [a] -> [a] -> [a]
++ (RCL -> String
forall a. Show a => a -> String
show RCL
bs)
  where gl :: String -> String -> String
gl String
a String
b = String
a String -> String -> String
forall a. [a] -> [a] -> [a]
++ String
"\n\t" String -> String -> String
forall a. [a] -> [a] -> [a]
++ String
b

rk_show2 :: String -> [Double] -> [RCL] -> RCL -> RCL -> String
rk_show2 :: String -> [Double] -> [RCL] -> RCL -> RCL -> String
rk_show2 String
title [Double]
cs [RCL]
as RCL
bs RCL
ds =
  String
title String -> String -> String
forall a. [a] -> [a] -> [a]
++ String
":\ncs:\t" String -> String -> String
forall a. [a] -> [a] -> [a]
++ ([Double] -> String
forall a. Show a => a -> String
show [Double]
cs) String -> String -> String
forall a. [a] -> [a] -> [a]
++ String
"\nas:\t" String -> String -> String
forall a. [a] -> [a] -> [a]
++
    (String -> String -> String) -> [String] -> String
forall a. (a -> a -> a) -> [a] -> a
foldl1' String -> String -> String
gl ((RCL -> String) -> [RCL] -> [String]
forall a b. (a -> b) -> [a] -> [b]
map RCL -> String
forall a. Show a => a -> String
show [RCL]
as) String -> String -> String
forall a. [a] -> [a] -> [a]
++ String
"\nbs:\t" String -> String -> String
forall a. [a] -> [a] -> [a]
++ (RCL -> String
forall a. Show a => a -> String
show RCL
bs) String -> String -> String
forall a. [a] -> [a] -> [a]
++
    String
"\nds:\t" String -> String -> String
forall a. [a] -> [a] -> [a]
++ (RCL -> String
forall a. Show a => a -> String
show RCL
ds)
  where gl :: String -> String -> String
gl String
a String
b = String
a String -> String -> String
forall a. [a] -> [a] -> [a]
++ String
"\n\t" String -> String -> String
forall a. [a] -> [a] -> [a]
++ String
b

-- Some specific explicit methods, taken from
-- "List of Runge-Kutta methods" at Wikipedia




-- | forward Euler: unconditionally unstable: don't use this!
-- If you do, your dangly bits will fall off!
rkfe
  :: (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)
rkfe :: (Double -> a -> a)
-> (a -> a -> a)
-> (Double -> a -> a)
-> Double
-> (Double, a)
-> (Double, a)
rkfe = [Double]
-> [RCL]
-> RCL
-> (Double -> a -> a)
-> (a -> a -> a)
-> (Double -> a -> a)
-> Double
-> (Double, a)
-> (Double, a)
forall a.
[Double]
-> [RCL]
-> RCL
-> (Double -> a -> a)
-> (a -> a -> a)
-> (Double -> a -> a)
-> Double
-> (Double, a)
-> (Double, a)
core1 [Double]
cs_fe [RCL]
as_fe RCL
bs_fe
show_rkfe :: String
show_rkfe :: String
show_rkfe = String -> [Double] -> [RCL] -> RCL -> String
rk_show1 String
"Forward Euler" [Double]
cs_fe [RCL]
as_fe RCL
bs_fe
cs_fe :: [Double]
cs_fe = [Rational] -> [Double]
ratToDbls [Rational
0]
as_fe :: [RCL]
as_fe = [[Rational]] -> [RCL]
ratToRCLs [[]]
bs_fe :: RCL
bs_fe = [Rational] -> RCL
ratToRCL  [Rational
1]

-- | Kutta's third-order method
rk3
  :: (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)
rk3 :: (Double -> a -> a)
-> (a -> a -> a)
-> (Double -> a -> a)
-> Double
-> (Double, a)
-> (Double, a)
rk3 = [Double]
-> [RCL]
-> RCL
-> (Double -> a -> a)
-> (a -> a -> a)
-> (Double -> a -> a)
-> Double
-> (Double, a)
-> (Double, a)
forall a.
[Double]
-> [RCL]
-> RCL
-> (Double -> a -> a)
-> (a -> a -> a)
-> (Double -> a -> a)
-> Double
-> (Double, a)
-> (Double, a)
core1 [Double]
cs_rk3 [RCL]
as_rk3 RCL
bs_rk3
show_rk3 :: String
show_rk3 :: String
show_rk3 = String -> [Double] -> [RCL] -> RCL -> String
rk_show1 String
"Kutta\'s third-order method" [Double]
cs_rk3 [RCL]
as_rk3 RCL
bs_rk3
cs_rk3 :: [Double]
cs_rk3 = [Rational] -> [Double]
ratToDbls [Rational
0, Integer
1Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
%Integer
2, Rational
1]
as_rk3 :: [RCL]
as_rk3 = [[Rational]] -> [RCL]
ratToRCLs [[], [Integer
1Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
%Integer
2], [-Rational
1, Rational
2]]
bs_rk3 :: RCL
bs_rk3 = [Rational] -> RCL
ratToRCL  [Integer
1Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
%Integer
6, Integer
2Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
%Integer
3, Integer
1Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
%Integer
6]

-- | Classic fourth-order method
rk4a
  :: (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)
rk4a :: (Double -> a -> a)
-> (a -> a -> a)
-> (Double -> a -> a)
-> Double
-> (Double, a)
-> (Double, a)
rk4a = [Double]
-> [RCL]
-> RCL
-> (Double -> a -> a)
-> (a -> a -> a)
-> (Double -> a -> a)
-> Double
-> (Double, a)
-> (Double, a)
forall a.
[Double]
-> [RCL]
-> RCL
-> (Double -> a -> a)
-> (a -> a -> a)
-> (Double -> a -> a)
-> Double
-> (Double, a)
-> (Double, a)
core1 [Double]
cs_rk4a [RCL]
as_rk4a RCL
bs_rk4a
show_rk4a :: String
show_rk4a :: String
show_rk4a = String -> [Double] -> [RCL] -> RCL -> String
rk_show1 String
"Classic fourth-order method" [Double]
cs_rk4a [RCL]
as_rk4a RCL
bs_rk4a
cs_rk4a :: [Double]
cs_rk4a = [Rational] -> [Double]
ratToDbls [Rational
0, Integer
1Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
%Integer
2, Integer
1Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
%Integer
2, Rational
1]
as_rk4a :: [RCL]
as_rk4a = [[Rational]] -> [RCL]
ratToRCLs [[], [Integer
1Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
%Integer
2], [Rational
0, Integer
1Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
%Integer
2], [Rational
0, Rational
0, Rational
1]]
bs_rk4a :: RCL
bs_rk4a = [Rational] -> RCL
ratToRCL  [Integer
1Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
%Integer
6, Integer
1Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
%Integer
3, Integer
1Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
%Integer
3, Integer
1Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
%Integer
6]


-- | Kutta's other fourth-order method... "The first [above] is more popular,
-- the second is more precise." (Hairer, Norsett, Wanner)
rk4b
  :: (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)
rk4b :: (Double -> a -> a)
-> (a -> a -> a)
-> (Double -> a -> a)
-> Double
-> (Double, a)
-> (Double, a)
rk4b = [Double]
-> [RCL]
-> RCL
-> (Double -> a -> a)
-> (a -> a -> a)
-> (Double -> a -> a)
-> Double
-> (Double, a)
-> (Double, a)
forall a.
[Double]
-> [RCL]
-> RCL
-> (Double -> a -> a)
-> (a -> a -> a)
-> (Double -> a -> a)
-> Double
-> (Double, a)
-> (Double, a)
core1 [Double]
cs_rk4b [RCL]
as_rk4b RCL
bs_rk4b
show_rk4b :: String
show_rk4b :: String
show_rk4b = String -> [Double] -> [RCL] -> RCL -> String
rk_show1 String
"Kutta's other classic fourth-order method" [Double]
cs_rk4b [RCL]
as_rk4b RCL
bs_rk4b
cs_rk4b :: [Double]
cs_rk4b = [Rational] -> [Double]
ratToDbls [Rational
0, Integer
1Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
%Integer
3, Integer
2Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
%Integer
3, Rational
1]
as_rk4b :: [RCL]
as_rk4b = [[Rational]] -> [RCL]
ratToRCLs [[], [Integer
1Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
%Integer
3], [-Integer
1Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
%Integer
3, Rational
1], [Rational
1, -Rational
1, Rational
1]]
bs_rk4b :: RCL
bs_rk4b = [Rational] -> RCL
ratToRCL  [Integer
1Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
%Integer
8, Integer
3Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
%Integer
8, Integer
3Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
%Integer
8, Integer
1Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
%Integer
8]

-- Some adaptive-stepsize methods, also from Wikipedia; more from HNW.
-- These don't auto-adapt, but they do allow the user to make a somewhat
-- intelligent decision about what the step size ought to be at each step.

-- Helper function to take the difference of two lists of rationals:
-- we don't want to use zipWith (-) because that gives only the head where
-- both lists have entries; we want implicit zeros at the end, as far as
-- is necessary.

diffs :: [a] -> [a] -> [a]
diffs [] [] = []
diffs [a]
xs [] = [a]
xs
diffs [] [a]
xs = (a -> a) -> [a] -> [a]
forall a b. (a -> b) -> [a] -> [b]
map a -> a
forall a. Num a => a -> a
negate [a]
xs
diffs (a
x:[a]
xs) (a
y:[a]
ys) = (a
x a -> a -> a
forall a. Num a => a -> a -> a
- a
y) a -> [a] -> [a]
forall a. a -> [a] -> [a]
: [a] -> [a] -> [a]
diffs [a]
xs [a]
ys



-- | Heun-Euler, order 2/1
rkhe
  :: (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)
rkhe :: (Double -> a -> a)
-> (a -> a -> a)
-> (Double -> a -> a)
-> Double
-> (Double, a)
-> (Double, a, a)
rkhe = [Double]
-> [RCL]
-> RCL
-> RCL
-> (Double -> a -> a)
-> (a -> a -> a)
-> (Double -> a -> a)
-> Double
-> (Double, a)
-> (Double, a, a)
forall a.
[Double]
-> [RCL]
-> RCL
-> RCL
-> (Double -> a -> a)
-> (a -> a -> a)
-> (Double -> a -> a)
-> Double
-> (Double, a)
-> (Double, a, a)
core2 [Double]
cs_he [RCL]
as_he RCL
bs_he RCL
ds_he
show_rkhe :: String
show_rkhe :: String
show_rkhe = String -> [Double] -> [RCL] -> RCL -> RCL -> String
rk_show2 String
"Heun-Euler 2(1)" [Double]
cs_he [RCL]
as_he RCL
bs_he RCL
ds_he
cs_he :: [Double]
cs_he = [Rational] -> [Double]
ratToDbls [Rational
0, Rational
1]
as_he :: [RCL]
as_he = [[Rational]] -> [RCL]
ratToRCLs [[], [Rational
1]]
r1_he :: [Rational]
r1_he = [Integer
1Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
%Integer
2, Integer
1Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
%Integer
2] -- second-order coeffs
r2_he :: [Rational]
r2_he = [Rational
1] -- first-order coeffs
bs_he :: RCL
bs_he = [Rational] -> RCL
ratToRCL [Rational]
r1_he
ds_he :: RCL
ds_he = [Rational] -> RCL
ratToRCL ([Rational] -> [Rational] -> [Rational]
forall a. Num a => [a] -> [a] -> [a]
diffs [Rational]
r1_he [Rational]
r2_he)

bs_he_aux :: RCL
bs_he_aux = [Rational] -> RCL
ratToRCL [Rational]
r2_he
rkhe_aux :: (Double -> a -> a)
-> (a -> a -> a)
-> (Double -> a -> a)
-> Double
-> (Double, a)
-> (Double, a)
rkhe_aux = [Double]
-> [RCL]
-> RCL
-> (Double -> a -> a)
-> (a -> a -> a)
-> (Double -> a -> a)
-> Double
-> (Double, a)
-> (Double, a)
forall a.
[Double]
-> [RCL]
-> RCL
-> (Double -> a -> a)
-> (a -> a -> a)
-> (Double -> a -> a)
-> Double
-> (Double, a)
-> (Double, a)
core1 [Double]
cs_he [RCL]
as_he RCL
bs_he_aux
show_rkhe_aux :: String
show_rkhe_aux = String -> [Double] -> [RCL] -> RCL -> String
rk_show1 String
"Heun-Euler (1)" [Double]
cs_he [RCL]
as_he RCL
bs_he_aux




-- | Bogacki-Shampine, order 3/2
rkbs
  :: (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)
rkbs :: (Double -> a -> a)
-> (a -> a -> a)
-> (Double -> a -> a)
-> Double
-> (Double, a)
-> (Double, a, a)
rkbs = [Double]
-> [RCL]
-> RCL
-> RCL
-> (Double -> a -> a)
-> (a -> a -> a)
-> (Double -> a -> a)
-> Double
-> (Double, a)
-> (Double, a, a)
forall a.
[Double]
-> [RCL]
-> RCL
-> RCL
-> (Double -> a -> a)
-> (a -> a -> a)
-> (Double -> a -> a)
-> Double
-> (Double, a)
-> (Double, a, a)
core2 [Double]
cs_bs [RCL]
as_bs RCL
bs_bs RCL
ds_bs
show_rkbs :: String
show_rkbs :: String
show_rkbs = String -> [Double] -> [RCL] -> RCL -> RCL -> String
rk_show2 String
"Bogacki-Shampine 3(2)" [Double]
cs_bs [RCL]
as_bs RCL
bs_bs RCL
ds_bs
cs_bs :: [Double]
cs_bs = [Rational] -> [Double]
ratToDbls [Rational
0, Integer
1Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
%Integer
2, Integer
3Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
%Integer
4, Rational
1]
as_bs :: [RCL]
as_bs = [[Rational]] -> [RCL]
ratToRCLs [[], [Integer
1Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
%Integer
2], [Rational
0, Integer
3Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
%Integer
4], [Integer
2Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
%Integer
9, Integer
1Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
%Integer
3, Integer
4Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
%Integer
9]]
r1_bs :: [Rational]
r1_bs = [Integer
2Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
%Integer
9, Integer
1Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
%Integer
3, Integer
4Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
%Integer
9] -- third-order coeffs
r2_bs :: [Rational]
r2_bs = [Integer
7Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
%Integer
24, Integer
1Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
%Integer
4, Integer
1Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
%Integer
3, Integer
1Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
%Integer
8]-- second-order coeffs
bs_bs :: RCL
bs_bs = [Rational] -> RCL
ratToRCL [Rational]
r1_bs
ds_bs :: RCL
ds_bs = [Rational] -> RCL
ratToRCL ([Rational] -> [Rational] -> [Rational]
forall a. Num a => [a] -> [a] -> [a]
diffs [Rational]
r1_bs [Rational]
r2_bs)

bs_bs_aux :: RCL
bs_bs_aux = [Rational] -> RCL
ratToRCL [Rational]
r2_bs
rkbs_aux :: (Double -> a -> a)
-> (a -> a -> a)
-> (Double -> a -> a)
-> Double
-> (Double, a)
-> (Double, a)
rkbs_aux = [Double]
-> [RCL]
-> RCL
-> (Double -> a -> a)
-> (a -> a -> a)
-> (Double -> a -> a)
-> Double
-> (Double, a)
-> (Double, a)
forall a.
[Double]
-> [RCL]
-> RCL
-> (Double -> a -> a)
-> (a -> a -> a)
-> (Double -> a -> a)
-> Double
-> (Double, a)
-> (Double, a)
core1 [Double]
cs_bs [RCL]
as_bs RCL
bs_bs_aux
show_rkbs_aux :: String
show_rkbs_aux = String -> [Double] -> [RCL] -> RCL -> String
rk_show1 String
"Bogacki-Shampine (2)" [Double]
cs_bs [RCL]
as_bs RCL
bs_bs_aux



-- | Runge-Kutta-Fehlberg, order 4/5
rkf45
  :: (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)
rkf45 :: (Double -> a -> a)
-> (a -> a -> a)
-> (Double -> a -> a)
-> Double
-> (Double, a)
-> (Double, a, a)
rkf45 = [Double]
-> [RCL]
-> RCL
-> RCL
-> (Double -> a -> a)
-> (a -> a -> a)
-> (Double -> a -> a)
-> Double
-> (Double, a)
-> (Double, a, a)
forall a.
[Double]
-> [RCL]
-> RCL
-> RCL
-> (Double -> a -> a)
-> (a -> a -> a)
-> (Double -> a -> a)
-> Double
-> (Double, a)
-> (Double, a, a)
core2 [Double]
cs_rkf [RCL]
as_rkf RCL
bs_rkf RCL
ds_rkf
show_rkf45 :: String
show_rkf45 :: String
show_rkf45 = String -> [Double] -> [RCL] -> RCL -> RCL -> String
rk_show2 String
"Runge-Kutta-Fehlberg 4(5)" [Double]
cs_rkf [RCL]
as_rkf RCL
bs_rkf RCL
ds_rkf
cs_rkf :: [Double]
cs_rkf = [Rational] -> [Double]
ratToDbls [Rational
0, Integer
1Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
%Integer
4, Integer
3Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
%Integer
8, Integer
12Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
%Integer
13, Rational
1, Integer
1Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
%Integer
2]
as_rkf :: [RCL]
as_rkf = [[Rational]] -> [RCL]
ratToRCLs [[],
                    [Integer
1Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
%Integer
4],
                    [Integer
3Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
%Integer
32, Integer
9Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
%Integer
32],
                    [Integer
1932Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
%Integer
2197, -Integer
7200Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
%Integer
2197, Integer
7296Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
%Integer
2197],
                    [Integer
439Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
%Integer
216, -Rational
8, Integer
3680Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
%Integer
513, -Integer
845Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
%Integer
4104],
                    [-Integer
8Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
%Integer
27, Rational
2, -Integer
3544Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
%Integer
2565, Integer
1859Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
%Integer
4104, -Integer
11Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
%Integer
40]]
-- fourth-order coeffs
r1_rkf :: [Rational]
r1_rkf = [Integer
25Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
%Integer
216, Rational
0, Integer
1408Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
%Integer
2565, Integer
2197Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
%Integer
4104, -Integer
1Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
%Integer
5]
-- fifth-order coeffs
r2_rkf :: [Rational]
r2_rkf = [Integer
16Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
%Integer
135, Rational
0, Integer
6656Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
%Integer
12825, Integer
28561Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
%Integer
56430, -Integer
9Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
%Integer
50, Integer
2Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
%Integer
55]
bs_rkf :: RCL
bs_rkf = [Rational] -> RCL
ratToRCL [Rational]
r1_rkf
ds_rkf :: RCL
ds_rkf = [Rational] -> RCL
ratToRCL ([Rational] -> [Rational] -> [Rational]
forall a. Num a => [a] -> [a] -> [a]
diffs [Rational]
r1_rkf [Rational]
r2_rkf)

bs_rkf_aux :: RCL
bs_rkf_aux = [Rational] -> RCL
ratToRCL [Rational]
r2_rkf
rkf45_aux :: (Double -> a -> a)
-> (a -> a -> a)
-> (Double -> a -> a)
-> Double
-> (Double, a)
-> (Double, a)
rkf45_aux = [Double]
-> [RCL]
-> RCL
-> (Double -> a -> a)
-> (a -> a -> a)
-> (Double -> a -> a)
-> Double
-> (Double, a)
-> (Double, a)
forall a.
[Double]
-> [RCL]
-> RCL
-> (Double -> a -> a)
-> (a -> a -> a)
-> (Double -> a -> a)
-> Double
-> (Double, a)
-> (Double, a)
core1 [Double]
cs_rkf [RCL]
as_rkf RCL
bs_rkf_aux
show_rkf45_aux :: String
show_rkf45_aux = String -> [Double] -> [RCL] -> RCL -> String
rk_show1 String
"Runge-Kutta-Fehlberg (5)" [Double]
cs_rkf [RCL]
as_rkf RCL
bs_rkf_aux




-- | Cash-Karp, order 4/5 (use 4th-order sol'n,
-- coeffs chosen to minimize error of 4th-order sol'n)
rkck
  :: (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)
rkck :: (Double -> a -> a)
-> (a -> a -> a)
-> (Double -> a -> a)
-> Double
-> (Double, a)
-> (Double, a, a)
rkck = [Double]
-> [RCL]
-> RCL
-> RCL
-> (Double -> a -> a)
-> (a -> a -> a)
-> (Double -> a -> a)
-> Double
-> (Double, a)
-> (Double, a, a)
forall a.
[Double]
-> [RCL]
-> RCL
-> RCL
-> (Double -> a -> a)
-> (a -> a -> a)
-> (Double -> a -> a)
-> Double
-> (Double, a)
-> (Double, a, a)
core2 [Double]
cs_ck [RCL]
as_ck RCL
bs_ck RCL
ds_ck
show_rkck :: String
show_rkck :: String
show_rkck = String -> [Double] -> [RCL] -> RCL -> RCL -> String
rk_show2 String
"Cash-Karp 4(5)" [Double]
cs_ck [RCL]
as_ck RCL
bs_ck RCL
ds_ck
cs_ck :: [Double]
cs_ck = [Rational] -> [Double]
ratToDbls [Rational
0, Integer
1Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
%Integer
5, Integer
3Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
%Integer
10, Integer
3Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
%Integer
5, Rational
1, Integer
7Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
%Integer
8]
as_ck :: [RCL]
as_ck = [[Rational]] -> [RCL]
ratToRCLs [[],
                   [Integer
1Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
%Integer
5],
                   [Integer
3Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
%Integer
40, Integer
9Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
%Integer
40],
                   [Integer
3Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
%Integer
10, -Integer
9Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
%Integer
10, Integer
6Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
%Integer
5],
                   [-Integer
11Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
%Integer
54, Integer
5Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
%Integer
2, -Integer
70Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
%Integer
27, Integer
35Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
%Integer
27],
                   [Integer
1631Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
%Integer
55296, Integer
175Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
%Integer
512, Integer
575Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
%Integer
13824, Integer
44275Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
%Integer
110592, Integer
253Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
%Integer
4096]]
-- fourth-order coeffs
r1_ck :: [Rational]
r1_ck = [Integer
2825Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
%Integer
27648, Rational
0, Integer
18575Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
%Integer
48384, Integer
13525Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
%Integer
55296, Integer
277Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
%Integer
14336, Integer
1Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
%Integer
4]
-- fifth-order coeffs
r2_ck :: [Rational]
r2_ck = [Integer
37Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
%Integer
378, Rational
0, Integer
250Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
%Integer
621, Integer
125Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
%Integer
594, Rational
0, Integer
512Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
%Integer
1771]
bs_ck :: RCL
bs_ck = [Rational] -> RCL
ratToRCL [Rational]
r1_ck
ds_ck :: RCL
ds_ck = [Rational] -> RCL
ratToRCL ([Rational] -> [Rational] -> [Rational]
forall a. Num a => [a] -> [a] -> [a]
diffs [Rational]
r1_ck [Rational]
r2_ck)

bs_ck_aux :: RCL
bs_ck_aux = [Rational] -> RCL
ratToRCL [Rational]
r2_ck
rkck_aux :: (Double -> a -> a)
-> (a -> a -> a)
-> (Double -> a -> a)
-> Double
-> (Double, a)
-> (Double, a)
rkck_aux = [Double]
-> [RCL]
-> RCL
-> (Double -> a -> a)
-> (a -> a -> a)
-> (Double -> a -> a)
-> Double
-> (Double, a)
-> (Double, a)
forall a.
[Double]
-> [RCL]
-> RCL
-> (Double -> a -> a)
-> (a -> a -> a)
-> (Double -> a -> a)
-> Double
-> (Double, a)
-> (Double, a)
core1 [Double]
cs_ck [RCL]
as_ck RCL
bs_ck_aux
show_rkck_aux :: String
show_rkck_aux = String -> [Double] -> [RCL] -> RCL -> String
rk_show1 String
"Cash-Karp (5)" [Double]
cs_ck [RCL]
as_ck RCL
bs_ck_aux




-- | 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
  :: (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)
rkdp :: (Double -> a -> a)
-> (a -> a -> a)
-> (Double -> a -> a)
-> Double
-> (Double, a)
-> (Double, a, a)
rkdp = [Double]
-> [RCL]
-> RCL
-> RCL
-> (Double -> a -> a)
-> (a -> a -> a)
-> (Double -> a -> a)
-> Double
-> (Double, a)
-> (Double, a, a)
forall a.
[Double]
-> [RCL]
-> RCL
-> RCL
-> (Double -> a -> a)
-> (a -> a -> a)
-> (Double -> a -> a)
-> Double
-> (Double, a)
-> (Double, a, a)
core2 [Double]
cs_dp [RCL]
as_dp RCL
bs_dp RCL
ds_dp
show_rkdp :: String
show_rkdp :: String
show_rkdp = String -> [Double] -> [RCL] -> RCL -> RCL -> String
rk_show2 String
"Dormand-Prince 5(4) \"DOPRI5\"" [Double]
cs_dp [RCL]
as_dp RCL
bs_dp RCL
ds_dp
cs_dp :: [Double]
cs_dp = [Rational] -> [Double]
ratToDbls [Rational
0, Integer
1Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
%Integer
5, Integer
3Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
%Integer
10, Integer
4Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
%Integer
5, Integer
8Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
%Integer
9, Rational
1, Rational
1]
as_dp :: [RCL]
as_dp = [[Rational]] -> [RCL]
ratToRCLs [[],
                   [Integer
1Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
%Integer
5],
                   [Integer
3Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
%Integer
40, Integer
9Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
%Integer
40],
                   [Integer
44Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
%Integer
45, -Integer
56Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
%Integer
15, Integer
32Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
%Integer
9],
                   [Integer
19372Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
%Integer
6561, -Integer
25360Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
%Integer
2187, Integer
64448Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
%Integer
6561, -Integer
212Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
%Integer
729],
                   [Integer
9017Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
%Integer
3168, -Integer
355Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
%Integer
33, Integer
46732Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
%Integer
5247, Integer
49Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
%Integer
176, -Integer
5103Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
%Integer
18656],
                   [Integer
35Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
%Integer
384, Rational
0, Integer
500Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
%Integer
1113, Integer
125Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
%Integer
192, -Integer
2187Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
%Integer
6784, Integer
11Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
%Integer
84]]
-- fifth-order coeffs
r1_dp :: [Rational]
r1_dp = [Integer
35Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
%Integer
384, Rational
0, Integer
500Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
%Integer
1113, Integer
125Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
%Integer
192, -Integer
2187Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
%Integer
6784, Integer
11Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
%Integer
84]
-- fourth-order coeffs
r2_dp :: [Rational]
r2_dp = [Integer
5179Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
%Integer
57600, Rational
0, Integer
7571Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
%Integer
16695, Integer
393Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
%Integer
640, -Integer
92097Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
%Integer
339200, Integer
187Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
%Integer
2100, Integer
1Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
%Integer
40]
bs_dp :: RCL
bs_dp = [Rational] -> RCL
ratToRCL [Rational]
r1_dp
ds_dp :: RCL
ds_dp = [Rational] -> RCL
ratToRCL ([Rational] -> [Rational] -> [Rational]
forall a. Num a => [a] -> [a] -> [a]
diffs [Rational]
r1_dp [Rational]
r2_dp)

bs_dp_aux :: RCL
bs_dp_aux = [Rational] -> RCL
ratToRCL [Rational]
r2_dp
rkdp_aux :: (Double -> a -> a)
-> (a -> a -> a)
-> (Double -> a -> a)
-> Double
-> (Double, a)
-> (Double, a)
rkdp_aux = [Double]
-> [RCL]
-> RCL
-> (Double -> a -> a)
-> (a -> a -> a)
-> (Double -> a -> a)
-> Double
-> (Double, a)
-> (Double, a)
forall a.
[Double]
-> [RCL]
-> RCL
-> (Double -> a -> a)
-> (a -> a -> a)
-> (Double -> a -> a)
-> Double
-> (Double, a)
-> (Double, a)
core1 [Double]
cs_dp [RCL]
as_dp RCL
bs_dp_aux
show_rkdp_aux :: String
show_rkdp_aux = String -> [Double] -> [RCL] -> RCL -> String
rk_show1 String
"Dormand-Prince (4)" [Double]
cs_dp [RCL]
as_dp RCL
bs_dp_aux




-- | 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
  :: (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)
rkf78 :: (Double -> a -> a)
-> (a -> a -> a)
-> (Double -> a -> a)
-> Double
-> (Double, a)
-> (Double, a, a)
rkf78 = [Double]
-> [RCL]
-> RCL
-> RCL
-> (Double -> a -> a)
-> (a -> a -> a)
-> (Double -> a -> a)
-> Double
-> (Double, a)
-> (Double, a, a)
forall a.
[Double]
-> [RCL]
-> RCL
-> RCL
-> (Double -> a -> a)
-> (a -> a -> a)
-> (Double -> a -> a)
-> Double
-> (Double, a)
-> (Double, a, a)
core2 [Double]
cs_f78 [RCL]
as_f78 RCL
bs_f78 RCL
ds_f78
show_rkf78 :: String
show_rkf78 :: String
show_rkf78 = String -> [Double] -> [RCL] -> RCL -> RCL -> String
rk_show2 String
"Fehlberg 7(8)" [Double]
cs_f78 [RCL]
as_f78 RCL
bs_f78 RCL
ds_f78
cs_f78 :: [Double]
cs_f78 = [Rational] -> [Double]
ratToDbls [Rational
0, Integer
2Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
%Integer
27, Integer
1Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
%Integer
9, Integer
1Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
%Integer
6, Integer
5Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
%Integer
12, Integer
1Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
%Integer
2, Integer
5Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
%Integer
6, Integer
1Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
%Integer
6, Integer
2Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
%Integer
3, Integer
1Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
%Integer
3, Rational
1, Rational
0, Rational
1]
as_f78 :: [RCL]
as_f78 =
  [[Rational]] -> [RCL]
ratToRCLs
    [[],
    [Integer
2Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
%Integer
27],
    [Integer
1Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
%Integer
36, Integer
1Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
%Integer
12],
    [Integer
1Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
%Integer
24, Rational
0, Integer
1Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
%Integer
8],
    [Integer
5Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
%Integer
12, Rational
0, -Integer
25Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
%Integer
16, Integer
25Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
%Integer
16],
    [Integer
1Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
%Integer
20, Rational
0, Rational
0, Integer
1Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
%Integer
4, Integer
1Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
%Integer
5],
    [-Integer
25Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
%Integer
108, Rational
0, Rational
0, Integer
125Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
%Integer
108, -Integer
65Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
%Integer
27, Integer
125Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
%Integer
54],
    [Integer
31Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
%Integer
300, Rational
0, Rational
0, Rational
0, Integer
61Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
%Integer
225, -Integer
2Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
%Integer
9, Integer
13Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
%Integer
900],
    [Rational
2, Rational
0, Rational
0, -Integer
53Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
%Integer
6, Integer
704Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
%Integer
45, -Integer
107Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
%Integer
9, Integer
67Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
%Integer
90, Rational
3],
    [-Integer
91Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
%Integer
108, Rational
0, Rational
0, Integer
23Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
%Integer
108, -Integer
976Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
%Integer
135, Integer
311Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
%Integer
54, -Integer
19Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
%Integer
60, Integer
17Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
%Integer
6, -Integer
1Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
%Integer
12],
    [Integer
2383Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
%Integer
4100, Rational
0, Rational
0, -Integer
341Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
%Integer
164, Integer
4496Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
%Integer
1025, -Integer
301Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
%Integer
82, Integer
2133Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
%Integer
4100, Integer
45Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
%Integer
82, Integer
45Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
%Integer
164, Integer
18Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
%Integer
41],
    [Integer
3Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
%Integer
205, Rational
0, Rational
0, Rational
0, Rational
0, -Integer
6Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
%Integer
41, -Integer
3Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
%Integer
205, -Integer
3Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
%Integer
41, Integer
3Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
%Integer
41, Integer
6Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
%Integer
41, Rational
0],
    [-Integer
1777Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
%Integer
4100, Rational
0, Rational
0, -Integer
341Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
%Integer
164, Integer
4496Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
%Integer
1025, -Integer
289Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
%Integer
82, Integer
2193Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
%Integer
4100, Integer
51Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
%Integer
82, Integer
33Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
%Integer
164, Integer
12Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
%Integer
41, Rational
0, Rational
1]]
-- seventh-order coeffs
r1_f78 :: [Rational]
r1_f78 = [Integer
41Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
%Integer
840, Rational
0, Rational
0, Rational
0, Rational
0, Integer
34Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
%Integer
105, Integer
9Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
%Integer
35, Integer
9Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
%Integer
35, Integer
9Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
%Integer
280, Integer
9Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
%Integer
280, Integer
41Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
%Integer
840]
-- eighth-order coeffs
r2_f78 :: [Rational]
r2_f78 = [Rational
0, Rational
0, Rational
0, Rational
0, Rational
0, Integer
34Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
%Integer
105, Integer
9Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
%Integer
35, Integer
9Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
%Integer
35, Integer
9Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
%Integer
280, Integer
9Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
%Integer
280, Rational
0, Integer
41Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
%Integer
840, Integer
41Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
%Integer
840]
bs_f78 :: RCL
bs_f78 = [Rational] -> RCL
ratToRCL [Rational]
r1_f78
ds_f78 :: RCL
ds_f78 = [Rational] -> RCL
ratToRCL ([Rational] -> [Rational] -> [Rational]
forall a. Num a => [a] -> [a] -> [a]
diffs [Rational]
r1_f78 [Rational]
r2_f78)

bs_f78_aux :: RCL
bs_f78_aux = [Rational] -> RCL
ratToRCL [Rational]
r2_f78
rkf78_aux :: (Double -> a -> a)
-> (a -> a -> a)
-> (Double -> a -> a)
-> Double
-> (Double, a)
-> (Double, a)
rkf78_aux = [Double]
-> [RCL]
-> RCL
-> (Double -> a -> a)
-> (a -> a -> a)
-> (Double -> a -> a)
-> Double
-> (Double, a)
-> (Double, a)
forall a.
[Double]
-> [RCL]
-> RCL
-> (Double -> a -> a)
-> (a -> a -> a)
-> (Double -> a -> a)
-> Double
-> (Double, a)
-> (Double, a)
core1 [Double]
cs_f78 [RCL]
as_f78 RCL
bs_f78_aux
show_rkf78_aux :: String
show_rkf78_aux = String -> [Double] -> [RCL] -> RCL -> String
rk_show1 String
"Fehlberg (8)" [Double]
cs_f78 [RCL]
as_f78 RCL
bs_f78_aux




-- | Verner, order 6/5 (DVERK)
rkv65
  :: (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)
rkv65 :: (Double -> a -> a)
-> (a -> a -> a)
-> (Double -> a -> a)
-> Double
-> (Double, a)
-> (Double, a, a)
rkv65 = [Double]
-> [RCL]
-> RCL
-> RCL
-> (Double -> a -> a)
-> (a -> a -> a)
-> (Double -> a -> a)
-> Double
-> (Double, a)
-> (Double, a, a)
forall a.
[Double]
-> [RCL]
-> RCL
-> RCL
-> (Double -> a -> a)
-> (a -> a -> a)
-> (Double -> a -> a)
-> Double
-> (Double, a)
-> (Double, a, a)
core2 [Double]
cs_v65 [RCL]
as_v65 RCL
bs_v65 RCL
ds_v65
show_rkv65 :: String
show_rkv65 = String -> [Double] -> [RCL] -> RCL -> RCL -> String
rk_show2 String
"Verner 6(5) \"DVERK\"" [Double]
cs_v65 [RCL]
as_v65 RCL
bs_v65 RCL
ds_v65
cs_v65 :: [Double]
cs_v65 = [Rational] -> [Double]
ratToDbls [Rational
0, Integer
1Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
%Integer
6, Integer
4Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
%Integer
15, Integer
2Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
%Integer
3, Integer
5Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
%Integer
6, Rational
1, Integer
1Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
%Integer
15, Rational
1]
as_v65 :: [RCL]
as_v65 =
  [[Rational]] -> [RCL]
ratToRCLs
    [[],
    [Integer
1Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
%Integer
6],
    [Integer
4Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
%Integer
75, Integer
16Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
%Integer
75],
    [Integer
5Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
%Integer
6, -Integer
8Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
%Integer
3, Integer
5Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
%Integer
2],
    [-Integer
165Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
%Integer
64, Integer
55Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
%Integer
6, -Integer
425Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
%Integer
64, Integer
85Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
%Integer
96],
    [Integer
12Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
%Integer
5, -Rational
8, Integer
4015Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
%Integer
612, -Integer
11Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
%Integer
36, Integer
88Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
%Integer
255],
    [-Integer
8263Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
%Integer
15000, Integer
124Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
%Integer
75, -Integer
643Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
%Integer
680, -Integer
81Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
%Integer
250, Integer
2484Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
%Integer
10625],
    [Integer
3501Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
%Integer
1720, -Integer
300Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
%Integer
43, Integer
297275Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
%Integer
52632, -Integer
319Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
%Integer
2322, Integer
24068Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
%Integer
84065, Rational
0, Integer
3850Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
%Integer
26703]]
-- sixth-order coeffs
r1_v65 :: [Rational]
r1_v65 = [Integer
3Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
%Integer
40, Rational
0, Integer
875Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
%Integer
2244, Integer
23Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
%Integer
72, Integer
264Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
%Integer
1955, Rational
0, Integer
125Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
%Integer
11592, Integer
43Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
%Integer
616]
-- fifth-order coeffs
r2_v65 :: [Rational]
r2_v65 = [Integer
13Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
%Integer
160, Rational
0, Integer
2375Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
%Integer
5984, Integer
5Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
%Integer
16, Integer
12Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
%Integer
85, Integer
3Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
%Integer
44]
bs_v65 :: RCL
bs_v65 = [Rational] -> RCL
ratToRCL [Rational]
r1_v65
ds_v65 :: RCL
ds_v65 = [Rational] -> RCL
ratToRCL ([Rational] -> [Rational] -> [Rational]
forall a. Num a => [a] -> [a] -> [a]
diffs [Rational]
r1_v65 [Rational]
r2_v65)

bs_v65_aux :: RCL
bs_v65_aux = [Rational] -> RCL
ratToRCL [Rational]
r2_v65
rkv65_aux :: (Double -> a -> a)
-> (a -> a -> a)
-> (Double -> a -> a)
-> Double
-> (Double, a)
-> (Double, a)
rkv65_aux = [Double]
-> [RCL]
-> RCL
-> (Double -> a -> a)
-> (a -> a -> a)
-> (Double -> a -> a)
-> Double
-> (Double, a)
-> (Double, a)
forall a.
[Double]
-> [RCL]
-> RCL
-> (Double -> a -> a)
-> (a -> a -> a)
-> (Double -> a -> a)
-> Double
-> (Double, a)
-> (Double, a)
core1 [Double]
cs_v65 [RCL]
as_v65 RCL
bs_v65_aux
show_rkv65_aux :: String
show_rkv65_aux = String -> [Double] -> [RCL] -> RCL -> String
rk_show1 String
"Verner (5)" [Double]
cs_v65 [RCL]
as_v65 RCL
bs_v65_aux