-- SPDX-License-Identifier: MPL-2.0

{- | Checked real ODE/DAE integration. Fixed-step methods use the requested
grid. Adaptive methods have explicit tolerances, step bounds and work budgets.
Callbacks must be pure, total and retain the state dimension. Dense output
uses a cubic Hermite interpolant, independently of the integration method.
-}
module HMath.ODE (
    ODE,
    Method (..),
    ODEOptions (..),
    defaultODEOptions,
    Trajectory (..),
    solveIVP,
    solveAdaptive,
    interpolateTrajectory,
    Event (..),
    EventDirection (..),
    locateEvents,
    solveBoundary,
    solveDAE,
    solveHamiltonian,
    HamiltonianMethod (..),
    solveNystrom,
    sensitivities,
    equilibrium,
    phaseTrajectories,
    lyapunovExponents,
) where

import Control.Monad (foldM)
import Data.List (transpose)
import HMath.Calculus (jacobian)
import HMath.Error (MathError (..))
import HMath.Matrix (Matrix, fromListChecked, toList)
import HMath.Matrix.LU (solve)
import HMath.ODE.Tableau

-- | Time and state callback, returning a derivative of exactly the same length.
type ODE = Double -> [Double] -> [Double]

{- | Explicit, multistep, implicit, and linearly implicit methods. RK45 is the
Cash-Karp 5(4) pair; DormandPrince is the Dormand-Prince 5(4) pair. Radau is
two-stage Radau IIA (order 3), BDF is order 2 with backward-Euler startup,
and Rosenbrock is linearly implicit Rosenbrock-Euler (order 1).
-}
data Method
    = Euler
    | Heun
    | Midpoint
    | RK2
    | RK4
    | RK23
    | RK45
    | DormandPrince
    | DOP853
    | AdamsBashforth
    | AdamsMoulton
    | PredictorCorrector
    | BackwardEuler
    | BDF
    | Radau
    | Rosenbrock
    deriving (Method -> Method -> Bool
(Method -> Method -> Bool)
-> (Method -> Method -> Bool) -> Eq Method
forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
$c== :: Method -> Method -> Bool
== :: Method -> Method -> Bool
$c/= :: Method -> Method -> Bool
/= :: Method -> Method -> Bool
Eq, Int -> Method -> ShowS
[Method] -> ShowS
Method -> String
(Int -> Method -> ShowS)
-> (Method -> String) -> ([Method] -> ShowS) -> Show Method
forall a.
(Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
$cshowsPrec :: Int -> Method -> ShowS
showsPrec :: Int -> Method -> ShowS
$cshow :: Method -> String
show :: Method -> String
$cshowList :: [Method] -> ShowS
showList :: [Method] -> ShowS
Show)

{- | Positive step bounds and tolerances; maximumSteps bounds accepted and
rejected attempts, and newtonIterations bounds each implicit correction.
-}
data ODEOptions = ODEOptions
    { ODEOptions -> Double
odeAbsoluteTolerance :: Double
    -- ^ Finite positive absolute local error tolerance.
    , ODEOptions -> Double
odeRelativeTolerance :: Double
    -- ^ Finite nonnegative local error tolerance relative to state scale.
    , ODEOptions -> Double
initialStep :: Double
    -- ^ Requested positive initial step magnitude; integration direction follows time bounds.
    , ODEOptions -> Double
minimumStep :: Double
    -- ^ Positive minimum permitted step magnitude before reporting failure.
    , ODEOptions -> Double
maximumStep :: Double
    -- ^ Positive maximum permitted step magnitude.
    , ODEOptions -> Int
maximumSteps :: Int
    -- ^ Maximum integration-step work budget.
    , ODEOptions -> Int
newtonIterations :: Int
    -- ^ Iteration budget for implicit-method nonlinear solves.
    }
    deriving (ODEOptions -> ODEOptions -> Bool
(ODEOptions -> ODEOptions -> Bool)
-> (ODEOptions -> ODEOptions -> Bool) -> Eq ODEOptions
forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
$c== :: ODEOptions -> ODEOptions -> Bool
== :: ODEOptions -> ODEOptions -> Bool
$c/= :: ODEOptions -> ODEOptions -> Bool
/= :: ODEOptions -> ODEOptions -> Bool
Eq, Int -> ODEOptions -> ShowS
[ODEOptions] -> ShowS
ODEOptions -> String
(Int -> ODEOptions -> ShowS)
-> (ODEOptions -> String)
-> ([ODEOptions] -> ShowS)
-> Show ODEOptions
forall a.
(Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
$cshowsPrec :: Int -> ODEOptions -> ShowS
showsPrec :: Int -> ODEOptions -> ShowS
$cshow :: ODEOptions -> String
show :: ODEOptions -> String
$cshowList :: [ODEOptions] -> ShowS
showList :: [ODEOptions] -> ShowS
Show)

-- | Conservative default adaptive settings for Double callbacks.
defaultODEOptions :: ODEOptions
defaultODEOptions :: ODEOptions
defaultODEOptions = Double
-> Double -> Double -> Double -> Double -> Int -> Int -> ODEOptions
ODEOptions Double
1e-9 Double
1e-7 Double
1e-2 Double
1e-12 Double
0.1 Int
100000 Int
30

{- | Strictly monotonic time samples, states and endpoint derivatives. The
constructor is public for inspection; interpolation validates all dimensions.
-}
data Trajectory = Trajectory
    { Trajectory -> [Double]
trajectoryTimes :: [Double]
    -- ^ Accepted sample times in integration order.
    , Trajectory -> [[Double]]
trajectoryStates :: [[Double]]
    -- ^ State vector at each trajectory time.
    , Trajectory -> [[Double]]
trajectoryDerivatives :: [[Double]]
    -- ^ ODE derivative vector at each trajectory sample, used by dense output.
    }
    deriving (Trajectory -> Trajectory -> Bool
(Trajectory -> Trajectory -> Bool)
-> (Trajectory -> Trajectory -> Bool) -> Eq Trajectory
forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
$c== :: Trajectory -> Trajectory -> Bool
== :: Trajectory -> Trajectory -> Bool
$c/= :: Trajectory -> Trajectory -> Bool
/= :: Trajectory -> Trajectory -> Bool
Eq, Int -> Trajectory -> ShowS
[Trajectory] -> ShowS
Trajectory -> String
(Int -> Trajectory -> ShowS)
-> (Trajectory -> String)
-> ([Trajectory] -> ShowS)
-> Show Trajectory
forall a.
(Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
$cshowsPrec :: Int -> Trajectory -> ShowS
showsPrec :: Int -> Trajectory -> ShowS
$cshow :: Trajectory -> String
show :: Trajectory -> String
$cshowList :: [Trajectory] -> ShowS
showList :: [Trajectory] -> ShowS
Show)

finite :: Double -> Bool
finite :: Double -> Bool
finite Double
x = Bool -> Bool
not (Double -> Bool
forall a. RealFloat a => a -> Bool
isNaN Double
x Bool -> Bool -> Bool
|| Double -> Bool
forall a. RealFloat a => a -> Bool
isInfinite Double
x)
require :: Bool -> String -> Either MathError ()
require :: Bool -> String -> Either MathError ()
require Bool
condition String
message = if Bool
condition then () -> Either MathError ()
forall a b. b -> Either a b
Right () else MathError -> Either MathError ()
forall a b. a -> Either a b
Left (String -> MathError
InvalidArgument String
message)
checked :: [Double] -> Either MathError [Double]
checked :: [Double] -> Either MathError [Double]
checked [Double]
xs = if (Double -> Bool) -> [Double] -> Bool
forall (t :: * -> *) a. Foldable t => (a -> Bool) -> t a -> Bool
all Double -> Bool
finite [Double]
xs then [Double] -> Either MathError [Double]
forall a b. b -> Either a b
Right [Double]
xs else MathError -> Either MathError [Double]
forall a b. a -> Either a b
Left MathError
NonFiniteResult
add :: [Double] -> [Double] -> [Double]
add :: [Double] -> [Double] -> [Double]
add = (Double -> Double -> Double) -> [Double] -> [Double] -> [Double]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith Double -> Double -> Double
forall a. Num a => a -> a -> a
(+)
scale :: Double -> [Double] -> [Double]
scale :: Double -> [Double] -> [Double]
scale Double
a = (Double -> Double) -> [Double] -> [Double]
forall a b. (a -> b) -> [a] -> [b]
map (Double
a Double -> Double -> Double
forall a. Num a => a -> a -> a
*)
norm :: [Double] -> Double
norm :: [Double] -> Double
norm [Double]
xs =
    let m :: Double
m = [Double] -> Double
forall a. Ord a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Ord a) => t a -> a
maximum (Double
0 Double -> [Double] -> [Double]
forall a. a -> [a] -> [a]
: (Double -> Double) -> [Double] -> [Double]
forall a b. (a -> b) -> [a] -> [b]
map Double -> Double
forall a. Num a => a -> a
abs [Double]
xs)
     in if Double
m Double -> Double -> Bool
forall a. Eq a => a -> a -> Bool
== Double
0 then Double
0 else Double
m Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double -> Double
forall a. Floating a => a -> a
sqrt ([Double] -> Double
forall a. Num a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum [(Double
x Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
m) Double -> Int -> Double
forall a b. (Num a, Integral b) => a -> b -> a
^ (Int
2 :: Int) | Double
x <- [Double]
xs])
evaluate :: ODE -> Double -> [Double] -> Either MathError [Double]
evaluate :: (Double -> [Double] -> [Double])
-> Double -> [Double] -> Either MathError [Double]
evaluate Double -> [Double] -> [Double]
f Double
t [Double]
y = do
    Bool -> String -> Either MathError ()
require (Double -> Bool
finite Double
t Bool -> Bool -> Bool
&& (Double -> Bool) -> [Double] -> Bool
forall (t :: * -> *) a. Foldable t => (a -> Bool) -> t a -> Bool
all Double -> Bool
finite [Double]
y) String
"nonfinite ODE point"
    let result :: [Double]
result = Double -> [Double] -> [Double]
f Double
t [Double]
y
    Bool -> String -> Either MathError ()
require ([Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
result Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== [Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
y) String
"ODE callback changed state dimension"
    [Double] -> Either MathError [Double]
checked [Double]
result
validateOptions :: ODEOptions -> Either MathError ()
validateOptions :: ODEOptions -> Either MathError ()
validateOptions ODEOptions
o =
    Bool -> String -> Either MathError ()
require
        ( (Double -> Bool) -> [Double] -> Bool
forall (t :: * -> *) a. Foldable t => (a -> Bool) -> t a -> Bool
all
            Double -> Bool
finite
            [ODEOptions -> Double
odeAbsoluteTolerance ODEOptions
o, ODEOptions -> Double
odeRelativeTolerance ODEOptions
o, ODEOptions -> Double
initialStep ODEOptions
o, ODEOptions -> Double
minimumStep ODEOptions
o, ODEOptions -> Double
maximumStep ODEOptions
o]
            Bool -> Bool -> Bool
&& ODEOptions -> Double
odeAbsoluteTolerance ODEOptions
o Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
> Double
0
            Bool -> Bool -> Bool
&& ODEOptions -> Double
odeRelativeTolerance ODEOptions
o Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
>= Double
0
            Bool -> Bool -> Bool
&& ODEOptions -> Double
minimumStep ODEOptions
o Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
> Double
0
            Bool -> Bool -> Bool
&& ODEOptions -> Double
initialStep ODEOptions
o Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
> Double
0
            Bool -> Bool -> Bool
&& ODEOptions -> Double
maximumStep ODEOptions
o Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
>= ODEOptions -> Double
minimumStep ODEOptions
o
            Bool -> Bool -> Bool
&& ODEOptions -> Int
maximumSteps ODEOptions
o Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
> Int
0
            Bool -> Bool -> Bool
&& ODEOptions -> Int
newtonIterations ODEOptions
o Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
> Int
0
        )
        String
"invalid ODE options"

{- | Integrate on n equal panels. Multistep methods use genuine derivative
history; implicit methods use checked finite-difference Newton iterations.
-}
solveIVP ::
    Method -> ODEOptions -> Int -> ODE -> Double -> Double -> [Double] -> Either MathError Trajectory
solveIVP :: Method
-> ODEOptions
-> Int
-> (Double -> [Double] -> [Double])
-> Double
-> Double
-> [Double]
-> Either MathError Trajectory
solveIVP Method
method ODEOptions
options Int
panels Double -> [Double] -> [Double]
f Double
start Double
end [Double]
initial = do
    ODEOptions -> Either MathError ()
validateOptions ODEOptions
options
    Bool -> String -> Either MathError ()
require
        (Int
panels Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
> Int
0 Bool -> Bool -> Bool
&& Double -> Bool
finite Double
start Bool -> Bool -> Bool
&& Double -> Bool
finite Double
end Bool -> Bool -> Bool
&& (Double -> Bool) -> [Double] -> Bool
forall (t :: * -> *) a. Foldable t => (a -> Bool) -> t a -> Bool
all Double -> Bool
finite [Double]
initial)
        String
"invalid IVP grid or initial state"
    Bool -> String -> Either MathError ()
require (Int
panels Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
<= ODEOptions -> Int
maximumSteps ODEOptions
options) String
"IVP grid exceeds step budget"
    d0 <- (Double -> [Double] -> [Double])
-> Double -> [Double] -> Either MathError [Double]
evaluate Double -> [Double] -> [Double]
f Double
start [Double]
initial
    if start == end
        then Right (Trajectory [start] [initial] [d0])
        else do
            let h = (Double
end Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
start) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
panels
                go Int
i Double
t [Double]
y Maybe [Double]
previous Maybe [Double]
previousDerivative [Double]
times [[Double]]
states [[Double]]
derivatives
                    | Int
i Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
panels = Trajectory -> Either MathError Trajectory
forall a b. b -> Either a b
Right ([Double] -> [[Double]] -> [[Double]] -> Trajectory
Trajectory ([Double] -> [Double]
forall a. [a] -> [a]
reverse [Double]
times) ([[Double]] -> [[Double]]
forall a. [a] -> [a]
reverse [[Double]]
states) ([[Double]] -> [[Double]]
forall a. [a] -> [a]
reverse [[Double]]
derivatives))
                    | Bool
otherwise = do
                        let nextTime :: Double
nextTime = if Int
i Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1 Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
panels then Double
end else Double
start Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral (Int
i Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1) Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
h
                            dt :: Double
dt = Double
nextTime Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
t
                        Bool -> String -> Either MathError ()
require (Double
nextTime Double -> Double -> Bool
forall a. Eq a => a -> a -> Bool
/= Double
t) String
"ODE grid step rounded away"
                        next <- case (Method
method, Maybe [Double]
previous, Maybe [Double]
previousDerivative) of
                            (Method
AdamsBashforth, Maybe [Double]
_, Just [Double]
oldD) -> [Double] -> Either MathError [Double]
checked ([Double] -> [Double] -> [Double]
add [Double]
y (Double -> [Double] -> [Double]
scale Double
dt ((Double -> Double -> Double) -> [Double] -> [Double] -> [Double]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith (\Double
a Double
b -> Double
1.5 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
a Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
0.5 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
b) ([[Double]]
derivatives [[Double]] -> Int -> [Double]
forall a. HasCallStack => [a] -> Int -> a
!! Int
0) [Double]
oldD)))
                            (Method
AdamsMoulton, Maybe [Double]
_, Maybe [Double]
_) -> ODEOptions
-> (Double -> [Double] -> [Double])
-> Double
-> Double
-> [Double]
-> Maybe [Double]
-> Bool
-> Either MathError [Double]
implicitStep ODEOptions
options Double -> [Double] -> [Double]
f Double
t Double
dt [Double]
y Maybe [Double]
forall a. Maybe a
Nothing Bool
True
                            (Method
PredictorCorrector, Maybe [Double]
_, Just [Double]
oldD) -> do
                                predicted <-
                                    [Double] -> Either MathError [Double]
checked ([Double] -> [Double] -> [Double]
add [Double]
y (Double -> [Double] -> [Double]
scale Double
dt ((Double -> Double -> Double) -> [Double] -> [Double] -> [Double]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith (\Double
a Double
b -> Double
1.5 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
a Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
0.5 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
b) ([[Double]]
derivatives [[Double]] -> Int -> [Double]
forall a. HasCallStack => [a] -> Int -> a
!! Int
0) [Double]
oldD)))
                                dn <- evaluate f nextTime predicted
                                checked (add y (scale (dt / 2) (add (derivatives !! 0) dn)))
                            (Method
BDF, Just [Double]
oldY, Maybe [Double]
_) -> ODEOptions
-> (Double -> [Double] -> [Double])
-> Double
-> Double
-> [Double]
-> Maybe [Double]
-> Bool
-> Either MathError [Double]
implicitStep ODEOptions
options Double -> [Double] -> [Double]
f Double
t Double
dt [Double]
y ([Double] -> Maybe [Double]
forall a. a -> Maybe a
Just [Double]
oldY) Bool
False
                            (Method
BDF, Maybe [Double]
Nothing, Maybe [Double]
_) -> ODEOptions
-> (Double -> [Double] -> [Double])
-> Double
-> Double
-> [Double]
-> Maybe [Double]
-> Bool
-> Either MathError [Double]
implicitStep ODEOptions
options Double -> [Double] -> [Double]
f Double
t Double
dt [Double]
y Maybe [Double]
forall a. Maybe a
Nothing Bool
False
                            (Method, Maybe [Double], Maybe [Double])
_ -> Method
-> ODEOptions
-> (Double -> [Double] -> [Double])
-> Double
-> Double
-> [Double]
-> Either MathError [Double]
fixedStep Method
method ODEOptions
options Double -> [Double] -> [Double]
f Double
t Double
dt [Double]
y
                        dn <- evaluate f nextTime next
                        go
                            (i + 1)
                            nextTime
                            next
                            (Just y)
                            (Just (derivatives !! 0))
                            (nextTime : times)
                            (next : states)
                            (dn : derivatives)
            go 0 start initial Nothing Nothing [start] [initial] [d0]

{- | Adaptive explicit RK integration. Euler/Heun/Midpoint/RK4 use step-doubling
error estimates. RK23/RK45/DormandPrince/DOP853 use their embedded estimates.
Multistep and implicit methods require 'solveIVP' to preserve their grid/history.
-}
solveAdaptive ::
    Method -> ODEOptions -> ODE -> Double -> Double -> [Double] -> Either MathError Trajectory
solveAdaptive :: Method
-> ODEOptions
-> (Double -> [Double] -> [Double])
-> Double
-> Double
-> [Double]
-> Either MathError Trajectory
solveAdaptive Method
method ODEOptions
options Double -> [Double] -> [Double]
f Double
start Double
end [Double]
initial = do
    ODEOptions -> Either MathError ()
validateOptions ODEOptions
options
    Bool -> String -> Either MathError ()
require
        (Method
method Method -> [Method] -> Bool
forall a. Eq a => a -> [a] -> Bool
forall (t :: * -> *) a. (Foldable t, Eq a) => a -> t a -> Bool
`elem` [Method
Euler, Method
Heun, Method
Midpoint, Method
RK2, Method
RK4, Method
RK23, Method
RK45, Method
DormandPrince, Method
DOP853])
        String
"adaptive driver requires an explicit RK method"
    Bool -> String -> Either MathError ()
require (Double -> Bool
finite Double
start Bool -> Bool -> Bool
&& Double -> Bool
finite Double
end Bool -> Bool -> Bool
&& (Double -> Bool) -> [Double] -> Bool
forall (t :: * -> *) a. Foldable t => (a -> Bool) -> t a -> Bool
all Double -> Bool
finite [Double]
initial) String
"invalid adaptive IVP"
    d0 <- (Double -> [Double] -> [Double])
-> Double -> [Double] -> Either MathError [Double]
evaluate Double -> [Double] -> [Double]
f Double
start [Double]
initial
    if start == end
        then Right (Trajectory [start] [initial] [d0])
        else do
            let direction = if Double
end Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
> Double
start then Double
1 else -Double
1
                initialH = Double
direction Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double -> Double -> Double
forall a. Ord a => a -> a -> a
min (ODEOptions -> Double
maximumStep ODEOptions
options) (Double -> Double -> Double
forall a. Ord a => a -> a -> a
max (ODEOptions -> Double
minimumStep ODEOptions
options) (ODEOptions -> Double
initialStep ODEOptions
options))
                go Int
attempts Double
t [Double]
y Double
h [Double]
times [[Double]]
states [[Double]]
derivatives
                    | Double
direction Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Double
end Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
t) Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
<= Double
0 =
                        Trajectory -> Either MathError Trajectory
forall a b. b -> Either a b
Right ([Double] -> [[Double]] -> [[Double]] -> Trajectory
Trajectory ([Double] -> [Double]
forall a. [a] -> [a]
reverse [Double]
times) ([[Double]] -> [[Double]]
forall a. [a] -> [a]
reverse [[Double]]
states) ([[Double]] -> [[Double]]
forall a. [a] -> [a]
reverse [[Double]]
derivatives))
                    | Int
attempts Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
>= ODEOptions -> Int
maximumSteps ODEOptions
options = MathError -> Either MathError Trajectory
forall a b. a -> Either a b
Left (String -> Int -> MathError
NonConvergence String
"adaptive ODE" Int
attempts)
                    | Bool
otherwise = do
                        let dt :: Double
dt = Double
direction Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double -> Double -> Double
forall a. Ord a => a -> a -> a
min (Double -> Double
forall a. Num a => a -> a
abs Double
h) (Double -> Double
forall a. Num a => a -> a
abs (Double
end Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
t))
                            nextTime :: Double
nextTime = if Double -> Double
forall a. Num a => a -> a
abs Double
dt Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
>= Double -> Double
forall a. Num a => a -> a
abs (Double
end Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
t) then Double
end else Double
t Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
dt
                        Bool -> String -> Either MathError ()
require (Double
nextTime Double -> Double -> Bool
forall a. Eq a => a -> a -> Bool
/= Double
t) String
"adaptive ODE step rounded away"
                        (candidate, errVector) <- Method
-> ODEOptions
-> (Double -> [Double] -> [Double])
-> Double
-> Double
-> [Double]
-> Either MathError ([Double], [Double])
oneStep Method
method ODEOptions
options Double -> [Double] -> [Double]
f Double
t Double
dt [Double]
y
                        let scaled =
                                (Double -> Double -> Double -> Double)
-> [Double] -> [Double] -> [Double] -> [Double]
forall a b c d. (a -> b -> c -> d) -> [a] -> [b] -> [c] -> [d]
zipWith3
                                    ( \Double
e Double
a Double
b -> Double -> Double
forall a. Num a => a -> a
abs Double
e Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (ODEOptions -> Double
odeAbsoluteTolerance ODEOptions
options Double -> Double -> Double
forall a. Num a => a -> a -> a
+ ODEOptions -> Double
odeRelativeTolerance ODEOptions
options Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double -> Double -> Double
forall a. Ord a => a -> a -> a
max (Double -> Double
forall a. Num a => a -> a
abs Double
a) (Double -> Double
forall a. Num a => a -> a
abs Double
b))
                                    )
                                    [Double]
errVector
                                    [Double]
y
                                    [Double]
candidate
                            err = [Double] -> Double
forall a. Ord a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Ord a) => t a -> a
maximum (Double
0 Double -> [Double] -> [Double]
forall a. a -> [a] -> [a]
: [Double]
scaled)
                            factor = if Double
err Double -> Double -> Bool
forall a. Eq a => a -> a -> Bool
== Double
0 then Double
5 else Double -> Double -> Double
forall a. Ord a => a -> a -> a
max Double
0.1 (Double -> Double -> Double
forall a. Ord a => a -> a -> a
min Double
5 (Double
0.9 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
err Double -> Double -> Double
forall a. Floating a => a -> a -> a
** (-Double
1 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral (Method -> Int
methodOrder Method
method))))
                            nextH = Double
direction Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double -> Double -> Double
forall a. Ord a => a -> a -> a
min (ODEOptions -> Double
maximumStep ODEOptions
options) (Double -> Double -> Double
forall a. Ord a => a -> a -> a
max (ODEOptions -> Double
minimumStep ODEOptions
options) (Double -> Double
forall a. Num a => a -> a
abs Double
dt Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
factor))
                        if not (finite err)
                            then Left NonFiniteResult
                            else
                                if err <= 1
                                    then do
                                        dn <- evaluate f nextTime candidate
                                        go
                                            (attempts + 1)
                                            nextTime
                                            candidate
                                            nextH
                                            (nextTime : times)
                                            (candidate : states)
                                            (dn : derivatives)
                                    else
                                        if abs dt <= minimumStep options
                                            then Left (NonConvergence "minimum ODE step" attempts)
                                            else go (attempts + 1) t y nextH times states derivatives
            go 0 start initial initialH [start] [initial] [d0]

methodOrder :: Method -> Int
methodOrder :: Method -> Int
methodOrder Method
method = case Method
method of
    Method
Euler -> Int
2
    Method
Heun -> Int
3
    Method
Midpoint -> Int
3
    Method
RK2 -> Int
3
    Method
RK4 -> Int
5
    Method
RK23 -> Int
3
    Method
RK45 -> Int
5
    Method
DormandPrince -> Int
5
    Method
DOP853 -> Int
8
    Method
_ -> Int
2

-- Explicit Runge-Kutta stages share the same checked callback and accumulation.
stages ::
    ODE -> Double -> Double -> [Double] -> [Double] -> [[Double]] -> Either MathError [[Double]]
stages :: (Double -> [Double] -> [Double])
-> Double
-> Double
-> [Double]
-> [Double]
-> [[Double]]
-> Either MathError [[Double]]
stages Double -> [Double] -> [Double]
f Double
t Double
h [Double]
y [Double]
cs [[Double]]
rows = ([[Double]] -> (Double, [Double]) -> Either MathError [[Double]])
-> [[Double]]
-> [(Double, [Double])]
-> Either MathError [[Double]]
forall (t :: * -> *) (m :: * -> *) b a.
(Foldable t, Monad m) =>
(b -> a -> m b) -> b -> t a -> m b
foldM [[Double]] -> (Double, [Double]) -> Either MathError [[Double]]
stage [] ([Double] -> [[Double]] -> [(Double, [Double])]
forall a b. [a] -> [b] -> [(a, b)]
zip [Double]
cs [[Double]]
rows)
  where
    stage :: [[Double]] -> (Double, [Double]) -> Either MathError [[Double]]
stage [[Double]]
known (Double
c, [Double]
coefficients) = do
        state <- [Double] -> Either MathError [Double]
checked (([Double] -> [Double] -> [Double])
-> [Double] -> [[Double]] -> [Double]
forall b a. (b -> a -> b) -> b -> [a] -> b
forall (t :: * -> *) b a.
Foldable t =>
(b -> a -> b) -> b -> t a -> b
foldl [Double] -> [Double] -> [Double]
add [Double]
y [Double -> [Double] -> [Double]
scale (Double
h Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
a) [Double]
k | (Double
a, [Double]
k) <- [Double] -> [[Double]] -> [(Double, [Double])]
forall a b. [a] -> [b] -> [(a, b)]
zip [Double]
coefficients [[Double]]
known])
        k <- evaluate f (t + c * h) state
        Right (known ++ [k])
combine :: Double -> [Double] -> [Double] -> [[Double]] -> Either MathError [Double]
combine :: Double
-> [Double] -> [Double] -> [[Double]] -> Either MathError [Double]
combine Double
h [Double]
y [Double]
weights [[Double]]
ks = [Double] -> Either MathError [Double]
checked (([Double] -> [Double] -> [Double])
-> [Double] -> [[Double]] -> [Double]
forall b a. (b -> a -> b) -> b -> [a] -> b
forall (t :: * -> *) b a.
Foldable t =>
(b -> a -> b) -> b -> t a -> b
foldl [Double] -> [Double] -> [Double]
add [Double]
y [Double -> [Double] -> [Double]
scale (Double
h Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
w) [Double]
k | (Double
w, [Double]
k) <- [Double] -> [[Double]] -> [(Double, [Double])]
forall a b. [a] -> [b] -> [(a, b)]
zip [Double]
weights [[Double]]
ks])

oneStep ::
    Method -> ODEOptions -> ODE -> Double -> Double -> [Double] -> Either MathError ([Double], [Double])
oneStep :: Method
-> ODEOptions
-> (Double -> [Double] -> [Double])
-> Double
-> Double
-> [Double]
-> Either MathError ([Double], [Double])
oneStep Method
method ODEOptions
options Double -> [Double] -> [Double]
f Double
t Double
h [Double]
y
    | Method
method Method -> [Method] -> Bool
forall a. Eq a => a -> [a] -> Bool
forall (t :: * -> *) a. (Foldable t, Eq a) => a -> t a -> Bool
`elem` [Method
BackwardEuler, Method
BDF, Method
AdamsMoulton] = do
        z <- ODEOptions
-> (Double -> [Double] -> [Double])
-> Double
-> Double
-> [Double]
-> Maybe [Double]
-> Bool
-> Either MathError [Double]
implicitStep ODEOptions
options Double -> [Double] -> [Double]
f Double
t Double
h [Double]
y Maybe [Double]
forall a. Maybe a
Nothing (Method
method Method -> Method -> Bool
forall a. Eq a => a -> a -> Bool
== Method
AdamsMoulton)
        Right (z, replicate (length y) 0)
    | Method
method Method -> Method -> Bool
forall a. Eq a => a -> a -> Bool
== Method
Rosenbrock = do
        d <- (Double -> [Double] -> [Double])
-> Double -> [Double] -> Either MathError [Double]
evaluate Double -> [Double] -> [Double]
f Double
t [Double]
y
        j <- jacobian (sqrt (2.220446049250313e-16) * max 1 (norm y)) (f t) y
        let n = [Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
y; values = Matrix Double -> [Double]
forall a. Matrix a -> [a]
toList Matrix Double
j
        matrix <-
            fromListChecked
                n
                n
                [(if i == k then 1 else 0) - h * values !! (i * n + k) | i <- [0 .. n - 1], k <- [0 .. n - 1]]
        delta <- solve matrix (scale h d)
        z <- checked (add y delta)
        Right (z, replicate n 0)
    | Method
method Method -> Method -> Bool
forall a. Eq a => a -> a -> Bool
== Method
Radau = do
        d <- (Double -> [Double] -> [Double])
-> Double -> [Double] -> Either MathError [Double]
evaluate Double -> [Double] -> [Double]
f Double
t [Double]
y
        let n = [Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
y
            residual [Double]
zs =
                let ([Double]
k1, [Double]
k2) = Int -> [Double] -> ([Double], [Double])
forall a. Int -> [a] -> ([a], [a])
splitAt Int
n [Double]
zs
                    y1 :: [Double]
y1 = [Double] -> [Double] -> [Double]
add [Double]
y (Double -> [Double] -> [Double]
scale Double
h ((Double -> Double -> Double) -> [Double] -> [Double] -> [Double]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith (\Double
a Double
b -> Double
5 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
a Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
12 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
b Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
12) [Double]
k1 [Double]
k2))
                    y2 :: [Double]
y2 = [Double] -> [Double] -> [Double]
add [Double]
y (Double -> [Double] -> [Double]
scale Double
h ((Double -> Double -> Double) -> [Double] -> [Double] -> [Double]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith (\Double
a Double
b -> Double
3 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
a Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
4 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
b Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
4) [Double]
k1 [Double]
k2))
                    v1 :: [Double]
v1 = Double -> [Double] -> [Double]
f (Double
t Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
h Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
3) [Double]
y1
                    v2 :: [Double]
v2 = Double -> [Double] -> [Double]
f (Double
t Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
h) [Double]
y2
                 in if [Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
v1 Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
/= Int
n Bool -> Bool -> Bool
|| [Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
v2 Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
/= Int
n
                        then Int -> Double -> [Double]
forall a. Int -> a -> [a]
replicate (Int
2 Int -> Int -> Int
forall a. Num a => a -> a -> a
* Int
n) (Double
0 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
0)
                        else (Double -> Double -> Double) -> [Double] -> [Double] -> [Double]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith (-) ([Double]
k1 [Double] -> [Double] -> [Double]
forall a. [a] -> [a] -> [a]
++ [Double]
k2) ([Double]
v1 [Double] -> [Double] -> [Double]
forall a. [a] -> [a] -> [a]
++ [Double]
v2)
        ks <- newton options residual (d ++ d)
        let (k1, k2) = splitAt n ks
        z <- checked (add y (scale h (zipWith (\Double
a Double
b -> Double
3 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
a Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
4 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
b Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
4) k1 k2)))
        Right (z, replicate n 0)
    | Bool
otherwise = do
        let ([Double]
cs, [[Double]]
rows, [Double]
weights, Maybe [Double]
low) = Method -> ([Double], [[Double]], [Double], Maybe [Double])
tableau Method
method
        ks <- (Double -> [Double] -> [Double])
-> Double
-> Double
-> [Double]
-> [Double]
-> [[Double]]
-> Either MathError [[Double]]
stages Double -> [Double] -> [Double]
f Double
t Double
h [Double]
y [Double]
cs [[Double]]
rows
        high <- combine h y weights ks
        if method == DOP853
            then do
                lastStage <- evaluate f (t + h) high
                let allStages = [[Double]]
ks [[Double]] -> [[Double]] -> [[Double]]
forall a. [a] -> [a] -> [a]
++ [[Double]
lastStage]
                    errorPart [Double]
coefficients = ([Double] -> [Double] -> [Double])
-> [Double] -> [[Double]] -> [Double]
forall b a. (b -> a -> b) -> b -> [a] -> b
forall (t :: * -> *) b a.
Foldable t =>
(b -> a -> b) -> b -> t a -> b
foldl [Double] -> [Double] -> [Double]
add (Int -> Double -> [Double]
forall a. Int -> a -> [a]
replicate ([Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
y) Double
0) [Double -> [Double] -> [Double]
scale Double
w [Double]
k | (Double
w, [Double]
k) <- [Double] -> [[Double]] -> [(Double, [Double])]
forall a b. [a] -> [b] -> [(a, b)]
zip [Double]
coefficients [[Double]]
allStages]
                    e5 = [Double] -> [Double]
errorPart [Double]
dop853E5
                    e3 = [Double] -> [Double]
errorPart [Double]
dop853E3
                    correction Double
a Double
b =
                        let m :: Double
m = Double -> Double -> Double
forall a. Ord a => a -> a -> a
max (Double -> Double
forall a. Num a => a -> a
abs Double
a) (Double
0.1 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double -> Double
forall a. Num a => a -> a
abs Double
b)
                            denominator :: Double
denominator = if Double
m Double -> Double -> Bool
forall a. Eq a => a -> a -> Bool
== Double
0 then Double
1 else Double -> Double
forall a. Floating a => a -> a
sqrt ((Double
a Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
m) Double -> Int -> Double
forall a b. (Num a, Integral b) => a -> b -> a
^ (Int
2 :: Int) Double -> Double -> Double
forall a. Num a => a -> a -> a
+ (Double
0.1 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
b Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
m) Double -> Int -> Double
forall a b. (Num a, Integral b) => a -> b -> a
^ (Int
2 :: Int))
                         in if Double
m Double -> Double -> Bool
forall a. Eq a => a -> a -> Bool
== Double
0 then Double
0 else Double
h Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
a Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Double -> Double
forall a. Num a => a -> a
abs Double
a Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
m) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
denominator
                err <- checked (zipWith correction e5 e3)
                Right (high, err)
            else case low of
                Just [Double]
lowWeights -> do lower <- Double
-> [Double] -> [Double] -> [[Double]] -> Either MathError [Double]
combine Double
h [Double]
y [Double]
lowWeights [[Double]]
ks; err <- checked (zipWith (-) high lower); Right (high, err)
                Maybe [Double]
Nothing -> do
                    halfKs <- (Double -> [Double] -> [Double])
-> Double
-> Double
-> [Double]
-> [Double]
-> [[Double]]
-> Either MathError [[Double]]
stages Double -> [Double] -> [Double]
f Double
t (Double
h Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
2) [Double]
y [Double]
cs [[Double]]
rows
                    middle <- combine (h / 2) y weights halfKs
                    nextKs <- stages f (t + h / 2) (h / 2) middle cs rows
                    refined <- combine (h / 2) middle weights nextKs
                    let denominator = Double
2 Double -> Double -> Double
forall a. Floating a => a -> a -> a
** Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral (Method -> Int
methodOrder Method
method Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1) Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
1
                    err <- checked (map (/ denominator) (zipWith (-) refined high))
                    Right (refined, err)

tableau :: Method -> ([Double], [[Double]], [Double], Maybe [Double])
tableau :: Method -> ([Double], [[Double]], [Double], Maybe [Double])
tableau Method
method = case Method
method of
    Method
Euler -> ([Double
0], [[]], [Double
1], Maybe [Double]
forall a. Maybe a
Nothing)
    Method
Heun -> ([Double
0, Double
1], [[], [Double
1]], [Double
1 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
2, Double
1 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
2], Maybe [Double]
forall a. Maybe a
Nothing)
    Method
Midpoint -> ([Double], [[Double]], [Double], Maybe [Double])
forall {a}. ([Double], [[Double]], [Double], Maybe a)
midpoint
    Method
RK2 -> ([Double], [[Double]], [Double], Maybe [Double])
forall {a}. ([Double], [[Double]], [Double], Maybe a)
midpoint
    Method
AdamsBashforth -> ([Double], [[Double]], [Double], Maybe [Double])
forall {a}. ([Double], [[Double]], [Double], Maybe a)
rk4
    Method
PredictorCorrector -> ([Double], [[Double]], [Double], Maybe [Double])
forall {a}. ([Double], [[Double]], [Double], Maybe a)
rk4
    Method
RK4 -> ([Double], [[Double]], [Double], Maybe [Double])
forall {a}. ([Double], [[Double]], [Double], Maybe a)
rk4
    Method
RK23 ->
        ( [Double
0, Double
1 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
2, Double
3 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
4, Double
1]
        , [[], [Double
1 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
2], [Double
0, Double
3 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
4], [Double
2 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
9, Double
1 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
3, Double
4 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
9]]
        , [Double
2 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
9, Double
1 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
3, Double
4 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
9, Double
0]
        , [Double] -> Maybe [Double]
forall a. a -> Maybe a
Just [Double
7 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
24, Double
1 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
4, Double
1 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
3, Double
1 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
8]
        )
    Method
RK45 ->
        ( [Double
0, Double
1 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
5, Double
3 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
10, Double
3 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
5, Double
1, Double
7 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
8]
        ,
            [ []
            , [Double
1 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
5]
            , [Double
3 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
40, Double
9 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
40]
            , [Double
3 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
10, -Double
9 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
10, Double
6 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
5]
            , [-Double
11 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
54, Double
5 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
2, -Double
70 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
27, Double
35 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
27]
            , [Double
1631 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
55296, Double
175 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
512, Double
575 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
13824, Double
44275 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
110592, Double
253 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
4096]
            ]
        , [Double
37 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
378, Double
0, Double
250 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
621, Double
125 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
594, Double
0, Double
512 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
1771]
        , [Double] -> Maybe [Double]
forall a. a -> Maybe a
Just [Double
2825 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
27648, Double
0, Double
18575 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
48384, Double
13525 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
55296, Double
277 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
14336, Double
1 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
4]
        )
    Method
DormandPrince ->
        ( [Double
0, Double
1 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
5, Double
3 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
10, Double
4 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
5, Double
8 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
9, Double
1, Double
1]
        ,
            [ []
            , [Double
1 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
5]
            , [Double
3 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
40, Double
9 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
40]
            , [Double
44 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
45, -Double
56 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
15, Double
32 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
9]
            , [Double
19372 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
6561, -Double
25360 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
2187, Double
64448 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
6561, -Double
212 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
729]
            , [Double
9017 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
3168, -Double
355 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
33, Double
46732 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
5247, Double
49 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
176, -Double
5103 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
18656]
            , [Double
35 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
384, Double
0, Double
500 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
1113, Double
125 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
192, -Double
2187 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
6784, Double
11 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
84]
            ]
        , [Double
35 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
384, Double
0, Double
500 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
1113, Double
125 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
192, -Double
2187 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
6784, Double
11 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
84, Double
0]
        , [Double] -> Maybe [Double]
forall a. a -> Maybe a
Just [Double
5179 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
57600, Double
0, Double
7571 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
16695, Double
393 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
640, -Double
92097 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
339200, Double
187 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
2100, Double
1 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
40]
        )
    Method
DOP853 -> ([Double]
dop853C, [[Double]]
dop853A, [Double]
dop853B, Maybe [Double]
forall a. Maybe a
Nothing)
    Method
_ -> ([Double], [[Double]], [Double], Maybe [Double])
forall {a}. ([Double], [[Double]], [Double], Maybe a)
rk4
  where
    midpoint :: ([Double], [[Double]], [Double], Maybe a)
midpoint = ([Double
0, Double
1 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
2], [[], [Double
1 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
2]], [Double
0, Double
1], Maybe a
forall a. Maybe a
Nothing)
    rk4 :: ([Double], [[Double]], [Double], Maybe a)
rk4 =
        ([Double
0, Double
1 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
2, Double
1 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
2, Double
1], [[], [Double
1 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
2], [Double
0, Double
1 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
2], [Double
0, Double
0, Double
1]], [Double
1 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
6, Double
1 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
3, Double
1 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
3, Double
1 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
6], Maybe a
forall a. Maybe a
Nothing)

fixedStep ::
    Method -> ODEOptions -> ODE -> Double -> Double -> [Double] -> Either MathError [Double]
fixedStep :: Method
-> ODEOptions
-> (Double -> [Double] -> [Double])
-> Double
-> Double
-> [Double]
-> Either MathError [Double]
fixedStep Method
method ODEOptions
options Double -> [Double] -> [Double]
f Double
t Double
h [Double]
y
    | Method
method Method -> [Method] -> Bool
forall a. Eq a => a -> [a] -> Bool
forall (t :: * -> *) a. (Foldable t, Eq a) => a -> t a -> Bool
`elem` [Method
Euler, Method
Heun, Method
Midpoint, Method
RK2, Method
RK4, Method
AdamsBashforth, Method
PredictorCorrector] =
        let ([Double]
cs, [[Double]]
rows, [Double]
weights, Maybe [Double]
_) = Method -> ([Double], [[Double]], [Double], Maybe [Double])
tableau Method
method in (Double -> [Double] -> [Double])
-> Double
-> Double
-> [Double]
-> [Double]
-> [[Double]]
-> Either MathError [[Double]]
stages Double -> [Double] -> [Double]
f Double
t Double
h [Double]
y [Double]
cs [[Double]]
rows Either MathError [[Double]]
-> ([[Double]] -> Either MathError [Double])
-> Either MathError [Double]
forall a b.
Either MathError a
-> (a -> Either MathError b) -> Either MathError b
forall (m :: * -> *) a b. Monad m => m a -> (a -> m b) -> m b
>>= Double
-> [Double] -> [Double] -> [[Double]] -> Either MathError [Double]
combine Double
h [Double]
y [Double]
weights
    | Bool
otherwise = ([Double], [Double]) -> [Double]
forall a b. (a, b) -> a
fst (([Double], [Double]) -> [Double])
-> Either MathError ([Double], [Double])
-> Either MathError [Double]
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
<$> Method
-> ODEOptions
-> (Double -> [Double] -> [Double])
-> Double
-> Double
-> [Double]
-> Either MathError ([Double], [Double])
oneStep Method
method ODEOptions
options Double -> [Double] -> [Double]
f Double
t Double
h [Double]
y

implicitStep ::
    ODEOptions ->
    ODE ->
    Double ->
    Double ->
    [Double] ->
    Maybe [Double] ->
    Bool ->
    Either MathError [Double]
implicitStep :: ODEOptions
-> (Double -> [Double] -> [Double])
-> Double
-> Double
-> [Double]
-> Maybe [Double]
-> Bool
-> Either MathError [Double]
implicitStep ODEOptions
options Double -> [Double] -> [Double]
f Double
t Double
h [Double]
y Maybe [Double]
old Bool
trapezoid = do
    d <- (Double -> [Double] -> [Double])
-> Double -> [Double] -> Either MathError [Double]
evaluate Double -> [Double] -> [Double]
f Double
t [Double]
y
    predicted <- checked (add y (scale h d))
    let residual [Double]
z =
            let dz :: [Double]
dz = Double -> [Double] -> [Double]
f (Double
t Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
h) [Double]
z
             in if [Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
dz Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
/= [Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
y
                    then Int -> Double -> [Double]
forall a. Int -> a -> [a]
replicate ([Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
y) (Double
0 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
0)
                    else case Maybe [Double]
old of
                        Just [Double]
previous -> (Double -> Double -> Double -> Double)
-> [Double] -> [Double] -> [Double] -> [Double]
forall a b c d. (a -> b -> c -> d) -> [a] -> [b] -> [c] -> [d]
zipWith3 (\Double
a Double
b Double
c -> Double
a Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
4 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
b Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
3 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
c Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
3) [Double]
z [Double]
y [Double]
previous [Double] -> [Double] -> [Double]
`add` Double -> [Double] -> [Double]
scale (-Double
2 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
h Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
3) [Double]
dz
                        Maybe [Double]
Nothing ->
                            if Bool
trapezoid
                                then (Double -> Double -> Double) -> [Double] -> [Double] -> [Double]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith (-) [Double]
z [Double]
y [Double] -> [Double] -> [Double]
`add` Double -> [Double] -> [Double]
scale (-Double
h Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
2) ([Double] -> [Double] -> [Double]
add [Double]
d [Double]
dz)
                                else (Double -> Double -> Double) -> [Double] -> [Double] -> [Double]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith (-) [Double]
z [Double]
y [Double] -> [Double] -> [Double]
`add` Double -> [Double] -> [Double]
scale (-Double
h) [Double]
dz
    newton options residual predicted

newton :: ODEOptions -> ([Double] -> [Double]) -> [Double] -> Either MathError [Double]
newton :: ODEOptions
-> ([Double] -> [Double]) -> [Double] -> Either MathError [Double]
newton ODEOptions
options [Double] -> [Double]
residual [Double]
initial = [Double] -> Either MathError [Double]
checked [Double]
initial Either MathError [Double]
-> Either MathError [Double] -> Either MathError [Double]
forall a b.
Either MathError a -> Either MathError b -> Either MathError b
forall (m :: * -> *) a b. Monad m => m a -> m b -> m b
>> Int -> [Double] -> Either MathError [Double]
go Int
0 [Double]
initial
  where
    go :: Int -> [Double] -> Either MathError [Double]
go Int
iteration [Double]
x = do
        let r :: [Double]
r = [Double] -> [Double]
residual [Double]
x
        Bool -> String -> Either MathError ()
require ([Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
r Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== [Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
x) String
"Newton residual dimension mismatch"
        _ <- [Double] -> Either MathError [Double]
checked [Double]
r
        if converged options x r
            then Right x
            else
                if iteration >= newtonIterations options
                    then Left (NonConvergence "implicit Newton" iteration)
                    else do
                        matrix <- jacobian (sqrt (2.220446049250313e-16) * maximum (1 : map abs x)) residual x
                        delta <- solve matrix (map negate r)
                        next <- checked (add x delta)
                        if next == x
                            then Left (NonConvergence "Newton floating-point resolution" iteration)
                            else go (iteration + 1) next

-- Compare componentwise after scaling, avoiding overflowing vector norms.
converged :: ODEOptions -> [Double] -> [Double] -> Bool
converged :: ODEOptions -> [Double] -> [Double] -> Bool
converged ODEOptions
options [Double]
x [Double]
residual = [Bool] -> Bool
forall (t :: * -> *). Foldable t => t Bool -> Bool
and ((Double -> Double -> Bool) -> [Double] -> [Double] -> [Bool]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith Double -> Double -> Bool
acceptable [Double]
x [Double]
residual)
  where
    acceptable :: Double -> Double -> Bool
acceptable Double
value Double
r =
        let magnitude :: Double
magnitude = Double -> Double -> Double
forall a. Ord a => a -> a -> a
max (ODEOptions -> Double
odeAbsoluteTolerance ODEOptions
options) (Double -> Double
forall a. Num a => a -> a
abs Double
value)
         in Double -> Double
forall a. Num a => a -> a
abs Double
r Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
magnitude
                Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
<= ODEOptions -> Double
odeAbsoluteTolerance ODEOptions
options Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
magnitude Double -> Double -> Double
forall a. Num a => a -> a -> a
+ ODEOptions -> Double
odeRelativeTolerance ODEOptions
options Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Double -> Double
forall a. Num a => a -> a
abs Double
value Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
magnitude)

{- | Cubic Hermite interpolation inside the trajectory's closed time range.
Extrapolation, malformed public trajectories and repeated sample times fail.
-}
interpolateTrajectory :: Trajectory -> Double -> Either MathError [Double]
interpolateTrajectory :: Trajectory -> Double -> Either MathError [Double]
interpolateTrajectory Trajectory
trajectory Double
t = do
    let ts :: [Double]
ts = Trajectory -> [Double]
trajectoryTimes Trajectory
trajectory
        ys :: [[Double]]
ys = Trajectory -> [[Double]]
trajectoryStates Trajectory
trajectory
        ds :: [[Double]]
ds = Trajectory -> [[Double]]
trajectoryDerivatives Trajectory
trajectory
        n :: Int
n = case [[Double]]
ys of [] -> Int
0; [Double]
first : [[Double]]
_ -> [Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
first
        increasing :: Bool
increasing = [Bool] -> Bool
forall (t :: * -> *). Foldable t => t Bool -> Bool
and ((Double -> Double -> Bool) -> [Double] -> [Double] -> [Bool]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
(<) [Double]
ts (Int -> [Double] -> [Double]
forall a. Int -> [a] -> [a]
drop Int
1 [Double]
ts))
        decreasing :: Bool
decreasing = [Bool] -> Bool
forall (t :: * -> *). Foldable t => t Bool -> Bool
and ((Double -> Double -> Bool) -> [Double] -> [Double] -> [Bool]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
(>) [Double]
ts (Int -> [Double] -> [Double]
forall a. Int -> [a] -> [a]
drop Int
1 [Double]
ts))
    Bool -> String -> Either MathError ()
require
        ( Bool -> Bool
not ([Double] -> Bool
forall a. [a] -> Bool
forall (t :: * -> *) a. Foldable t => t a -> Bool
null [Double]
ts)
            Bool -> Bool -> Bool
&& [Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
ts Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== [[Double]] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [[Double]]
ys
            Bool -> Bool -> Bool
&& [Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
ts Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== [[Double]] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [[Double]]
ds
            Bool -> Bool -> Bool
&& ([Double] -> Bool) -> [[Double]] -> Bool
forall (t :: * -> *) a. Foldable t => (a -> Bool) -> t a -> Bool
all (\[Double]
v -> [Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
v Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
n Bool -> Bool -> Bool
&& (Double -> Bool) -> [Double] -> Bool
forall (t :: * -> *) a. Foldable t => (a -> Bool) -> t a -> Bool
all Double -> Bool
finite [Double]
v) ([[Double]]
ys [[Double]] -> [[Double]] -> [[Double]]
forall a. [a] -> [a] -> [a]
++ [[Double]]
ds)
            Bool -> Bool -> Bool
&& (Double -> Bool) -> [Double] -> Bool
forall (t :: * -> *) a. Foldable t => (a -> Bool) -> t a -> Bool
all Double -> Bool
finite [Double]
ts
            Bool -> Bool -> Bool
&& Double -> Bool
finite Double
t
            Bool -> Bool -> Bool
&& (Bool
increasing Bool -> Bool -> Bool
|| Bool
decreasing)
        )
        String
"malformed dense trajectory"
    Bool -> String -> Either MathError ()
require
        (Double
t Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
>= Double -> Double -> Double
forall a. Ord a => a -> a -> a
min ([Double]
ts [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
0) ([Double] -> Double
forall a. HasCallStack => [a] -> a
last [Double]
ts) Bool -> Bool -> Bool
&& Double
t Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
<= Double -> Double -> Double
forall a. Ord a => a -> a -> a
max ([Double]
ts [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
0) ([Double] -> Double
forall a. HasCallStack => [a] -> a
last [Double]
ts))
        String
"dense-output time is outside trajectory"
    if Double
t Double -> Double -> Bool
forall a. Eq a => a -> a -> Bool
== [Double] -> Double
forall a. HasCallStack => [a] -> a
last [Double]
ts
        then [Double] -> Either MathError [Double]
forall a b. b -> Either a b
Right ([[Double]] -> [Double]
forall a. HasCallStack => [a] -> a
last [[Double]]
ys)
        else
            if [Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
ts Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
1
                then [Double] -> Either MathError [Double]
forall a b. b -> Either a b
Right ([[Double]]
ys [[Double]] -> Int -> [Double]
forall a. HasCallStack => [a] -> Int -> a
!! Int
0)
                else do
                    let intervals :: [(Double, Double, ([Double], [Double]), ([Double], [Double]))]
intervals = [Double]
-> [Double]
-> [([Double], [Double])]
-> [([Double], [Double])]
-> [(Double, Double, ([Double], [Double]), ([Double], [Double]))]
forall a b c d. [a] -> [b] -> [c] -> [d] -> [(a, b, c, d)]
zip4 [Double]
ts (Int -> [Double] -> [Double]
forall a. Int -> [a] -> [a]
drop Int
1 [Double]
ts) ([[Double]] -> [[Double]] -> [([Double], [Double])]
forall a b. [a] -> [b] -> [(a, b)]
zip [[Double]]
ys [[Double]]
ds) ([[Double]] -> [[Double]] -> [([Double], [Double])]
forall a b. [a] -> [b] -> [(a, b)]
zip (Int -> [[Double]] -> [[Double]]
forall a. Int -> [a] -> [a]
drop Int
1 [[Double]]
ys) (Int -> [[Double]] -> [[Double]]
forall a. Int -> [a] -> [a]
drop Int
1 [[Double]]
ds))
                        selected :: [(Double, Double, [Double], [Double], [Double], [Double])]
selected = [(Double
a, Double
b, [Double]
y0, [Double]
d0, [Double]
y1, [Double]
d1) | (Double
a, Double
b, ([Double]
y0, [Double]
d0), ([Double]
y1, [Double]
d1)) <- [(Double, Double, ([Double], [Double]), ([Double], [Double]))]
intervals, Double
t Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
>= Double -> Double -> Double
forall a. Ord a => a -> a -> a
min Double
a Double
b Bool -> Bool -> Bool
&& Double
t Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
<= Double -> Double -> Double
forall a. Ord a => a -> a -> a
max Double
a Double
b]
                    case [(Double, Double, [Double], [Double], [Double], [Double])]
selected of
                        [] -> MathError -> Either MathError [Double]
forall a b. a -> Either a b
Left (String -> MathError
InvalidArgument String
"dense output has no containing interval")
                        (Double
a, Double
b, [Double]
y0, [Double]
d0, [Double]
y1, [Double]
d1) : [(Double, Double, [Double], [Double], [Double], [Double])]
_ -> do
                            let h :: Double
h = Double
b Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
a
                                s :: Double
s = (Double
t Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
a) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
h
                                w0 :: Double
w0 = Double
2 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
s Double -> Int -> Double
forall a b. (Num a, Integral b) => a -> b -> a
^ (Int
3 :: Int) Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
3 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
s Double -> Int -> Double
forall a b. (Num a, Integral b) => a -> b -> a
^ (Int
2 :: Int) Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
1
                                v0 :: Double
v0 = Double
s Double -> Int -> Double
forall a b. (Num a, Integral b) => a -> b -> a
^ (Int
3 :: Int) Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
2 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
s Double -> Int -> Double
forall a b. (Num a, Integral b) => a -> b -> a
^ (Int
2 :: Int) Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
s
                                w1 :: Double
w1 = -Double
2 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
s Double -> Int -> Double
forall a b. (Num a, Integral b) => a -> b -> a
^ (Int
3 :: Int) Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
3 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
s Double -> Int -> Double
forall a b. (Num a, Integral b) => a -> b -> a
^ (Int
2 :: Int)
                                v1 :: Double
v1 = Double
s Double -> Int -> Double
forall a b. (Num a, Integral b) => a -> b -> a
^ (Int
3 :: Int) Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
s Double -> Int -> Double
forall a b. (Num a, Integral b) => a -> b -> a
^ (Int
2 :: Int)
                            [Double] -> Either MathError [Double]
checked (([Double] -> [Double] -> [Double])
-> [Double] -> [[Double]] -> [Double]
forall b a. (b -> a -> b) -> b -> [a] -> b
forall (t :: * -> *) b a.
Foldable t =>
(b -> a -> b) -> b -> t a -> b
foldl [Double] -> [Double] -> [Double]
add (Int -> Double -> [Double]
forall a. Int -> a -> [a]
replicate Int
n Double
0) [Double -> [Double] -> [Double]
scale Double
w0 [Double]
y0, Double -> [Double] -> [Double]
scale (Double
h Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
v0) [Double]
d0, Double -> [Double] -> [Double]
scale Double
w1 [Double]
y1, Double -> [Double] -> [Double]
scale (Double
h Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
v1) [Double]
d1])

zip4 :: [a] -> [b] -> [c] -> [d] -> [(a, b, c, d)]
zip4 :: forall a b c d. [a] -> [b] -> [c] -> [d] -> [(a, b, c, d)]
zip4 (a
a : [a]
as) (b
b : [b]
bs) (c
c : [c]
cs) (d
d : [d]
ds) = (a
a, b
b, c
c, d
d) (a, b, c, d) -> [(a, b, c, d)] -> [(a, b, c, d)]
forall a. a -> [a] -> [a]
: [a] -> [b] -> [c] -> [d] -> [(a, b, c, d)]
forall a b c d. [a] -> [b] -> [c] -> [d] -> [(a, b, c, d)]
zip4 [a]
as [b]
bs [c]
cs [d]
ds
zip4 [a]
_ [b]
_ [c]
_ [d]
_ = []

-- | Event direction in increasing physical time, even for backward integration.
data EventDirection = AnyDirection | Rising | Falling deriving (EventDirection -> EventDirection -> Bool
(EventDirection -> EventDirection -> Bool)
-> (EventDirection -> EventDirection -> Bool) -> Eq EventDirection
forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
$c== :: EventDirection -> EventDirection -> Bool
== :: EventDirection -> EventDirection -> Bool
$c/= :: EventDirection -> EventDirection -> Bool
/= :: EventDirection -> EventDirection -> Bool
Eq, Int -> EventDirection -> ShowS
[EventDirection] -> ShowS
EventDirection -> String
(Int -> EventDirection -> ShowS)
-> (EventDirection -> String)
-> ([EventDirection] -> ShowS)
-> Show EventDirection
forall a.
(Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
$cshowsPrec :: Int -> EventDirection -> ShowS
showsPrec :: Int -> EventDirection -> ShowS
$cshow :: EventDirection -> String
show :: EventDirection -> String
$cshowList :: [EventDirection] -> ShowS
showList :: [EventDirection] -> ShowS
Show)

-- | A located event time and its dense-interpolated state.
data Event = Event
    { Event -> Double
eventTime :: Double
    -- ^ Located event crossing time.
    , Event -> [Double]
eventState :: [Double]
    -- ^ Interpolated state vector at eventTime.
    }
    deriving (Event -> Event -> Bool
(Event -> Event -> Bool) -> (Event -> Event -> Bool) -> Eq Event
forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
$c== :: Event -> Event -> Bool
== :: Event -> Event -> Bool
$c/= :: Event -> Event -> Bool
/= :: Event -> Event -> Bool
Eq, Int -> Event -> ShowS
[Event] -> ShowS
Event -> String
(Int -> Event -> ShowS)
-> (Event -> String) -> ([Event] -> ShowS) -> Show Event
forall a.
(Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
$cshowsPrec :: Int -> Event -> ShowS
showsPrec :: Int -> Event -> ShowS
$cshow :: Event -> String
show :: Event -> String
$cshowList :: [Event] -> ShowS
showList :: [Event] -> ShowS
Show)

{- | Locate bracketed sign-changing events by dense-output bisection. Tangential
roots and multiple roots inside one sample interval are not guaranteed.
-}
locateEvents ::
    EventDirection -> Double -> (Double -> [Double] -> Double) -> Trajectory -> Either MathError [Event]
locateEvents :: EventDirection
-> Double
-> (Double -> [Double] -> Double)
-> Trajectory
-> Either MathError [Event]
locateEvents EventDirection
direction Double
tolerance Double -> [Double] -> Double
function Trajectory
trajectory = do
    Bool -> String -> Either MathError ()
require (Double -> Bool
finite Double
tolerance Bool -> Bool -> Bool
&& Double
tolerance Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
> Double
0) String
"event time tolerance must be positive"
    let ts :: [Double]
ts = Trajectory -> [Double]
trajectoryTimes Trajectory
trajectory
    Bool -> String -> Either MathError ()
require (Bool -> Bool
not ([Double] -> Bool
forall a. [a] -> Bool
forall (t :: * -> *) a. Foldable t => t a -> Bool
null [Double]
ts)) String
"empty event trajectory"
    _ <- Trajectory -> Double -> Either MathError [Double]
interpolateTrajectory Trajectory
trajectory ([Double]
ts [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
0)
    pairs <-
        traverse
            (\Double
t -> do y <- Trajectory -> Double -> Either MathError [Double]
interpolateTrajectory Trajectory
trajectory Double
t; value <- scalar (function t y); Right (t, value))
            ts
    events <- traverse locate (zip pairs (drop 1 pairs))
    Right (deduplicate (concat events))
  where
    scalar :: Double -> Either MathError Double
scalar Double
x = if Double -> Bool
finite Double
x then Double -> Either MathError Double
forall a b. b -> Either a b
Right Double
x else MathError -> Either MathError Double
forall a b. a -> Either a b
Left MathError
NonFiniteResult
    accepts :: a -> a -> Bool
accepts a
a a
b = case EventDirection
direction of EventDirection
AnyDirection -> Bool
True; EventDirection
Rising -> a
b a -> a -> Bool
forall a. Ord a => a -> a -> Bool
> a
a; EventDirection
Falling -> a
b a -> a -> Bool
forall a. Ord a => a -> a -> Bool
< a
a
    locate :: ((Double, Double), (Double, Double)) -> Either MathError [Event]
locate ((Double
a, Double
fa), (Double
b, Double
fb)) = do
        let (Double
left, Double
right, Double
fl, Double
fr) = if Double
a Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< Double
b then (Double
a, Double
b, Double
fa, Double
fb) else (Double
b, Double
a, Double
fb, Double
fa)
        if Bool -> Bool
not (Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
accepts Double
fl Double
fr) Bool -> Bool -> Bool
|| (Double
fl Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
> Double
0 Bool -> Bool -> Bool
&& Double
fr Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
> Double
0) Bool -> Bool -> Bool
|| (Double
fl Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< Double
0 Bool -> Bool -> Bool
&& Double
fr Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< Double
0)
            then [Event] -> Either MathError [Event]
forall a b. b -> Either a b
Right []
            else
                if Double
fl Double -> Double -> Bool
forall a. Eq a => a -> a -> Bool
== Double
0
                    then do y <- Trajectory -> Double -> Either MathError [Double]
interpolateTrajectory Trajectory
trajectory Double
left; Right [Event left y]
                    else
                        if Double
fr Double -> Double -> Bool
forall a. Eq a => a -> a -> Bool
== Double
0
                            then do y <- Trajectory -> Double -> Either MathError [Double]
interpolateTrajectory Trajectory
trajectory Double
right; Right [Event right y]
                            else do
                                root <- Int -> Double -> Double -> Double -> Either MathError Double
bisect (Int
0 :: Int) Double
left Double
right Double
fl
                                y <- interpolateTrajectory trajectory root
                                Right [Event root y]
    bisect :: Int -> Double -> Double -> Double -> Either MathError Double
bisect Int
count Double
left Double
right Double
fl
        | Int
count Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
>= Int
100 = MathError -> Either MathError Double
forall a b. a -> Either a b
Left (String -> Int -> MathError
NonConvergence String
"event bisection" Int
count)
        | Double
right Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
left Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
<= Double
tolerance = Double -> Either MathError Double
forall a b. b -> Either a b
Right (Double
left Double -> Double -> Double
forall a. Num a => a -> a -> a
+ (Double
right Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
left) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
2)
        | Bool
otherwise = do
            let middle :: Double
middle = Double
left Double -> Double -> Double
forall a. Num a => a -> a -> a
+ (Double
right Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
left) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
2
            y <- Trajectory -> Double -> Either MathError [Double]
interpolateTrajectory Trajectory
trajectory Double
middle
            fm <- scalar (function middle y)
            if fm == 0
                then Right middle
                else
                    if signum fm == signum fl
                        then bisect (count + 1) middle right fm
                        else bisect (count + 1) left middle fl
    deduplicate :: [Event] -> [Event]
deduplicate [] = []
    deduplicate (Event
event : [Event]
events) =
        Event
event
            Event -> [Event] -> [Event]
forall a. a -> [a] -> [a]
: [Event] -> [Event]
deduplicate ((Event -> Bool) -> [Event] -> [Event]
forall a. (a -> Bool) -> [a] -> [a]
dropWhile (\Event
next -> Double -> Double
forall a. Num a => a -> a
abs (Event -> Double
eventTime Event
next Double -> Double -> Double
forall a. Num a => a -> a -> a
- Event -> Double
eventTime Event
event) Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
<= Double
tolerance) [Event]
events)

{- | Shooting boundary-value solver: all initial-state components are unknown.
The boundary callback takes initial/final states and returns n residuals.
Returns the trajectory matching a locally converged shooting root; convergence
is local and depends on the supplied initial guess.
-}
solveBoundary ::
    Method ->
    ODEOptions ->
    Int ->
    ODE ->
    ([Double] -> [Double] -> [Double]) ->
    Double ->
    Double ->
    [Double] ->
    Either MathError Trajectory
solveBoundary :: Method
-> ODEOptions
-> Int
-> (Double -> [Double] -> [Double])
-> ([Double] -> [Double] -> [Double])
-> Double
-> Double
-> [Double]
-> Either MathError Trajectory
solveBoundary Method
method ODEOptions
options Int
panels Double -> [Double] -> [Double]
f [Double] -> [Double] -> [Double]
boundary Double
start Double
end [Double]
guess = do
    ODEOptions -> Either MathError ()
validateOptions ODEOptions
options
    let residual :: [Double] -> [Double]
residual [Double]
initial = case Method
-> ODEOptions
-> Int
-> (Double -> [Double] -> [Double])
-> Double
-> Double
-> [Double]
-> Either MathError Trajectory
solveIVP Method
method ODEOptions
options Int
panels Double -> [Double] -> [Double]
f Double
start Double
end [Double]
initial of
            Left MathError
_ -> Int -> Double -> [Double]
forall a. Int -> a -> [a]
replicate ([Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
initial) (Double
0 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
0)
            Right Trajectory
trajectory -> [Double] -> [Double] -> [Double]
boundary [Double]
initial ([[Double]] -> [Double]
forall a. HasCallStack => [a] -> a
last (Trajectory -> [[Double]]
trajectoryStates Trajectory
trajectory))
    initial <- ODEOptions
-> ([Double] -> [Double]) -> [Double] -> Either MathError [Double]
newton ODEOptions
options [Double] -> [Double]
residual [Double]
guess
    solveIVP method options panels f start end initial

{- | Index-one implicit DAE driver F(t,y,y')=0 with backward differences.
Initial y and y' must be consistent to the configured tolerance. Each
subsequent state is obtained by a full residual Newton solve.
-}
solveDAE ::
    ODEOptions ->
    Int ->
    (Double -> [Double] -> [Double] -> [Double]) ->
    Double ->
    Double ->
    [Double] ->
    [Double] ->
    Either MathError Trajectory
solveDAE :: ODEOptions
-> Int
-> (Double -> [Double] -> [Double] -> [Double])
-> Double
-> Double
-> [Double]
-> [Double]
-> Either MathError Trajectory
solveDAE ODEOptions
options Int
panels Double -> [Double] -> [Double] -> [Double]
residual Double
start Double
end [Double]
initial [Double]
initialDerivative = do
    ODEOptions -> Either MathError ()
validateOptions ODEOptions
options
    Bool -> String -> Either MathError ()
require
        ( Int
panels Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
> Int
0
            Bool -> Bool -> Bool
&& Int
panels Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
<= ODEOptions -> Int
maximumSteps ODEOptions
options
            Bool -> Bool -> Bool
&& Double -> Bool
finite Double
start
            Bool -> Bool -> Bool
&& Double -> Bool
finite Double
end
            Bool -> Bool -> Bool
&& [Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
initial Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== [Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
initialDerivative
            Bool -> Bool -> Bool
&& (Double -> Bool) -> [Double] -> Bool
forall (t :: * -> *) a. Foldable t => (a -> Bool) -> t a -> Bool
all Double -> Bool
finite ([Double]
initial [Double] -> [Double] -> [Double]
forall a. [a] -> [a] -> [a]
++ [Double]
initialDerivative)
        )
        String
"invalid DAE input"
    let r0 :: [Double]
r0 = Double -> [Double] -> [Double] -> [Double]
residual Double
start [Double]
initial [Double]
initialDerivative
    Bool -> String -> Either MathError ()
require ([Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
r0 Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== [Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
initial Bool -> Bool -> Bool
&& (Double -> Bool) -> [Double] -> Bool
forall (t :: * -> *) a. Foldable t => (a -> Bool) -> t a -> Bool
all Double -> Bool
finite [Double]
r0) String
"invalid initial DAE residual"
    Bool -> String -> Either MathError ()
require (ODEOptions -> [Double] -> [Double] -> Bool
converged ODEOptions
options [Double]
initial [Double]
r0) String
"inconsistent initial DAE state/derivative"
    if Double
start Double -> Double -> Bool
forall a. Eq a => a -> a -> Bool
== Double
end
        then Trajectory -> Either MathError Trajectory
forall a b. b -> Either a b
Right ([Double] -> [[Double]] -> [[Double]] -> Trajectory
Trajectory [Double
start] [[Double]
initial] [[Double]
initialDerivative])
        else do
            let h :: Double
h = (Double
end Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
start) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
panels
                go :: Int
-> Double
-> [Double]
-> [Double]
-> [Double]
-> [[Double]]
-> [[Double]]
-> Either MathError Trajectory
go Int
i Double
t [Double]
y [Double]
d [Double]
times [[Double]]
states [[Double]]
derivatives
                    | Int
i Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
panels = Trajectory -> Either MathError Trajectory
forall a b. b -> Either a b
Right ([Double] -> [[Double]] -> [[Double]] -> Trajectory
Trajectory ([Double] -> [Double]
forall a. [a] -> [a]
reverse [Double]
times) ([[Double]] -> [[Double]]
forall a. [a] -> [a]
reverse [[Double]]
states) ([[Double]] -> [[Double]]
forall a. [a] -> [a]
reverse [[Double]]
derivatives))
                    | Bool
otherwise = do
                        let nextTime :: Double
nextTime = if Int
i Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1 Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
panels then Double
end else Double
start Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral (Int
i Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1) Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
h; dt :: Double
dt = Double
nextTime Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
t
                        Bool -> String -> Either MathError ()
require (Double
dt Double -> Double -> Bool
forall a. Eq a => a -> a -> Bool
/= Double
0) String
"DAE step rounded away"
                        let equation :: [Double] -> [Double]
equation [Double]
z = Double -> [Double] -> [Double] -> [Double]
residual Double
nextTime [Double]
z ((Double -> Double) -> [Double] -> [Double]
forall a b. (a -> b) -> [a] -> [b]
map (Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
dt) ((Double -> Double -> Double) -> [Double] -> [Double] -> [Double]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith (-) [Double]
z [Double]
y))
                        next <- ODEOptions
-> ([Double] -> [Double]) -> [Double] -> Either MathError [Double]
newton ODEOptions
options [Double] -> [Double]
equation ([Double] -> [Double] -> [Double]
add [Double]
y (Double -> [Double] -> [Double]
scale Double
dt [Double]
d))
                        dn <- checked (map (/ dt) (zipWith (-) next y))
                        go (i + 1) nextTime next dn (nextTime : times) (next : states) (dn : derivatives)
            Int
-> Double
-> [Double]
-> [Double]
-> [Double]
-> [[Double]]
-> [[Double]]
-> Either MathError Trajectory
go Int
0 Double
start [Double]
initial [Double]
initialDerivative [Double
start] [[Double]
initial] [[Double]
initialDerivative]

-- | Separable Hamiltonian methods for q'=velocity(p), p'=force(q).
data HamiltonianMethod = SymplecticEuler | Verlet | Leapfrog deriving (HamiltonianMethod -> HamiltonianMethod -> Bool
(HamiltonianMethod -> HamiltonianMethod -> Bool)
-> (HamiltonianMethod -> HamiltonianMethod -> Bool)
-> Eq HamiltonianMethod
forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
$c== :: HamiltonianMethod -> HamiltonianMethod -> Bool
== :: HamiltonianMethod -> HamiltonianMethod -> Bool
$c/= :: HamiltonianMethod -> HamiltonianMethod -> Bool
/= :: HamiltonianMethod -> HamiltonianMethod -> Bool
Eq, Int -> HamiltonianMethod -> ShowS
[HamiltonianMethod] -> ShowS
HamiltonianMethod -> String
(Int -> HamiltonianMethod -> ShowS)
-> (HamiltonianMethod -> String)
-> ([HamiltonianMethod] -> ShowS)
-> Show HamiltonianMethod
forall a.
(Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
$cshowsPrec :: Int -> HamiltonianMethod -> ShowS
showsPrec :: Int -> HamiltonianMethod -> ShowS
$cshow :: HamiltonianMethod -> String
show :: HamiltonianMethod -> String
$cshowList :: [HamiltonianMethod] -> ShowS
showList :: [HamiltonianMethod] -> ShowS
Show)

{- | Fixed-step symplectic integration; output states concatenate q then p.
Verlet is kick-drift-kick; Leapfrog is drift-kick-drift with integer-time
output momenta. These equivalent second-order splittings have distinct order.
-}
solveHamiltonian ::
    HamiltonianMethod ->
    Int ->
    ([Double] -> [Double]) ->
    ([Double] -> [Double]) ->
    Double ->
    Double ->
    [Double] ->
    [Double] ->
    Either MathError Trajectory
solveHamiltonian :: HamiltonianMethod
-> Int
-> ([Double] -> [Double])
-> ([Double] -> [Double])
-> Double
-> Double
-> [Double]
-> [Double]
-> Either MathError Trajectory
solveHamiltonian HamiltonianMethod
method Int
panels [Double] -> [Double]
velocity [Double] -> [Double]
force Double
start Double
end [Double]
q0 [Double]
p0 = do
    Bool -> String -> Either MathError ()
require
        (Int
panels Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
> Int
0 Bool -> Bool -> Bool
&& [Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
q0 Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== [Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
p0 Bool -> Bool -> Bool
&& (Double -> Bool) -> [Double] -> Bool
forall (t :: * -> *) a. Foldable t => (a -> Bool) -> t a -> Bool
all Double -> Bool
finite ([Double]
q0 [Double] -> [Double] -> [Double]
forall a. [a] -> [a] -> [a]
++ [Double]
p0) Bool -> Bool -> Bool
&& Double -> Bool
finite Double
start Bool -> Bool -> Bool
&& Double -> Bool
finite Double
end)
        String
"invalid Hamiltonian initial state or grid"
    let n :: Int
n = [Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
q0
        evaluateVector :: (p -> [Double]) -> p -> Either MathError [Double]
evaluateVector p -> [Double]
f p
x =
            let y :: [Double]
y = p -> [Double]
f p
x
             in if [Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
y Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
/= Int
n
                    then MathError -> Either MathError [Double]
forall a b. a -> Either a b
Left (String -> MathError
InvalidArgument String
"Hamiltonian callback dimension mismatch")
                    else [Double] -> Either MathError [Double]
checked [Double]
y
        deriv :: [Double] -> [Double] -> Either MathError [Double]
deriv [Double]
q [Double]
p = do v <- ([Double] -> [Double]) -> [Double] -> Either MathError [Double]
forall {p}. (p -> [Double]) -> p -> Either MathError [Double]
evaluateVector [Double] -> [Double]
velocity [Double]
p; g <- evaluateVector force q; Right (v ++ g)
        step :: Double
-> [Double] -> [Double] -> Either MathError ([Double], [Double])
step Double
h [Double]
q [Double]
p = case HamiltonianMethod
method of
            HamiltonianMethod
SymplecticEuler -> do
                g <- ([Double] -> [Double]) -> [Double] -> Either MathError [Double]
forall {p}. (p -> [Double]) -> p -> Either MathError [Double]
evaluateVector [Double] -> [Double]
force [Double]
q
                p1 <- checked (add p (scale h g))
                v <- evaluateVector velocity p1
                q1 <- checked (add q (scale h v))
                Right (q1, p1)
            HamiltonianMethod
Verlet -> do
                g <- ([Double] -> [Double]) -> [Double] -> Either MathError [Double]
forall {p}. (p -> [Double]) -> p -> Either MathError [Double]
evaluateVector [Double] -> [Double]
force [Double]
q
                pm <- checked (add p (scale (h / 2) g))
                v <- evaluateVector velocity pm
                q1 <- checked (add q (scale h v))
                g1 <- evaluateVector force q1
                p1 <- checked (add pm (scale (h / 2) g1))
                Right (q1, p1)
            HamiltonianMethod
Leapfrog -> do
                v <- ([Double] -> [Double]) -> [Double] -> Either MathError [Double]
forall {p}. (p -> [Double]) -> p -> Either MathError [Double]
evaluateVector [Double] -> [Double]
velocity [Double]
p
                qm <- checked (add q (scale (h / 2) v))
                g <- evaluateVector force qm
                p1 <- checked (add p (scale h g))
                v1 <- evaluateVector velocity p1
                q1 <- checked (add qm (scale (h / 2) v1))
                Right (q1, p1)
        go :: Int
-> Double
-> [Double]
-> [Double]
-> [Double]
-> [[Double]]
-> [[Double]]
-> Either MathError Trajectory
go Int
i Double
t [Double]
q [Double]
p [Double]
ts [[Double]]
ys [[Double]]
ds
            | Int
i Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
panels = Trajectory -> Either MathError Trajectory
forall a b. b -> Either a b
Right ([Double] -> [[Double]] -> [[Double]] -> Trajectory
Trajectory ([Double] -> [Double]
forall a. [a] -> [a]
reverse [Double]
ts) ([[Double]] -> [[Double]]
forall a. [a] -> [a]
reverse [[Double]]
ys) ([[Double]] -> [[Double]]
forall a. [a] -> [a]
reverse [[Double]]
ds))
            | Bool
otherwise = do
                let nextTime :: Double
nextTime =
                        if Int
i Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1 Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
panels then Double
end else Double
start Double -> Double -> Double
forall a. Num a => a -> a -> a
+ (Double
end Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
start) Double -> Double -> Double
forall a. Num a => a -> a -> a
* Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral (Int
i Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
panels
                Bool -> String -> Either MathError ()
require (Double
nextTime Double -> Double -> Bool
forall a. Eq a => a -> a -> Bool
/= Double
t) String
"Hamiltonian step rounded away"
                (q1, p1) <- Double
-> [Double] -> [Double] -> Either MathError ([Double], [Double])
step (Double
nextTime Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
t) [Double]
q [Double]
p
                d1 <- deriv q1 p1
                go (i + 1) nextTime q1 p1 (nextTime : ts) ((q1 ++ p1) : ys) (d1 : ds)
    d0 <- [Double] -> [Double] -> Either MathError [Double]
deriv [Double]
q0 [Double]
p0
    if start == end
        then Right (Trajectory [start] [q0 ++ p0] [d0])
        else go 0 start q0 p0 [start] [q0 ++ p0] [d0]

{- | Fourth-order explicit Runge-Kutta-Nystrom method for q''=a(t,q), with
acceleration independent of velocity. State output concatenates q and q'.
-}
solveNystrom ::
    Int ->
    (Double -> [Double] -> [Double]) ->
    Double ->
    Double ->
    [Double] ->
    [Double] ->
    Either MathError Trajectory
solveNystrom :: Int
-> (Double -> [Double] -> [Double])
-> Double
-> Double
-> [Double]
-> [Double]
-> Either MathError Trajectory
solveNystrom Int
panels Double -> [Double] -> [Double]
acceleration Double
start Double
end [Double]
q0 [Double]
v0 = do
    Bool -> String -> Either MathError ()
require
        (Int
panels Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
> Int
0 Bool -> Bool -> Bool
&& Double -> Bool
finite Double
start Bool -> Bool -> Bool
&& Double -> Bool
finite Double
end Bool -> Bool -> Bool
&& [Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
q0 Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== [Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
v0 Bool -> Bool -> Bool
&& (Double -> Bool) -> [Double] -> Bool
forall (t :: * -> *) a. Foldable t => (a -> Bool) -> t a -> Bool
all Double -> Bool
finite ([Double]
q0 [Double] -> [Double] -> [Double]
forall a. [a] -> [a] -> [a]
++ [Double]
v0))
        String
"invalid Nystrom initial state or grid"
    let n :: Int
n = [Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
q0
        acc :: Double -> [Double] -> Either MathError [Double]
acc Double
t [Double]
q =
            let a :: [Double]
a = Double -> [Double] -> [Double]
acceleration Double
t [Double]
q
             in if [Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
a Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
/= Int
n then MathError -> Either MathError [Double]
forall a b. a -> Either a b
Left (String -> MathError
InvalidArgument String
"Nystrom callback dimension mismatch") else [Double] -> Either MathError [Double]
checked [Double]
a
        go :: Int
-> Double
-> [Double]
-> [Double]
-> [Double]
-> [[Double]]
-> [[Double]]
-> Either MathError Trajectory
go Int
i Double
t [Double]
q [Double]
v [Double]
ts [[Double]]
ys [[Double]]
ds
            | Int
i Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
panels = Trajectory -> Either MathError Trajectory
forall a b. b -> Either a b
Right ([Double] -> [[Double]] -> [[Double]] -> Trajectory
Trajectory ([Double] -> [Double]
forall a. [a] -> [a]
reverse [Double]
ts) ([[Double]] -> [[Double]]
forall a. [a] -> [a]
reverse [[Double]]
ys) ([[Double]] -> [[Double]]
forall a. [a] -> [a]
reverse [[Double]]
ds))
            | Bool
otherwise = do
                let nextTime :: Double
nextTime =
                        if Int
i Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1 Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
panels then Double
end else Double
start Double -> Double -> Double
forall a. Num a => a -> a -> a
+ (Double
end Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
start) Double -> Double -> Double
forall a. Num a => a -> a -> a
* Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral (Int
i Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
panels
                    h :: Double
h = Double
nextTime Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
t
                Bool -> String -> Either MathError ()
require (Double
h Double -> Double -> Bool
forall a. Eq a => a -> a -> Bool
/= Double
0) String
"Nystrom step rounded away"
                k1 <- Double -> [Double] -> Either MathError [Double]
acc Double
t [Double]
q
                qm <- checked (foldl add q [scale (h / 2) v, scale (h * h / 8) k1])
                k2 <- acc (t + h / 2) qm
                qe <- checked (foldl add q [scale h v, scale (h * h / 2) k2])
                k3 <- acc nextTime qe
                q1 <- checked (foldl add q [scale h v, scale (h * h / 6) (add k1 (scale 2 k2))])
                v1 <- checked (add v (scale (h / 6) (foldl add k1 [scale 4 k2, k3])))
                a1 <- acc nextTime q1
                go (i + 1) nextTime q1 v1 (nextTime : ts) ((q1 ++ v1) : ys) ((v1 ++ a1) : ds)
    a0 <- Double -> [Double] -> Either MathError [Double]
acc Double
start [Double]
q0
    if start == end
        then Right (Trajectory [start] [q0 ++ v0] [v0 ++ a0])
        else go 0 start q0 v0 [start] [q0 ++ v0] [v0 ++ a0]

{- | Terminal-state parameter sensitivities by centered perturbation of the
complete IVP solve. Parameters use the same explicit finite-difference step.
-}
sensitivities ::
    Method ->
    ODEOptions ->
    Int ->
    Double ->
    ([Double] -> ODE) ->
    [Double] ->
    Double ->
    Double ->
    [Double] ->
    Either MathError (Matrix Double)
sensitivities :: Method
-> ODEOptions
-> Int
-> Double
-> ([Double] -> Double -> [Double] -> [Double])
-> [Double]
-> Double
-> Double
-> [Double]
-> Either MathError (Matrix Double)
sensitivities Method
method ODEOptions
options Int
panels Double
h [Double] -> Double -> [Double] -> [Double]
family [Double]
parameters Double
start Double
end [Double]
initial = do
    Bool -> String -> Either MathError ()
require (Double -> Bool
finite Double
h Bool -> Bool -> Bool
&& Double
h Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
> Double
0 Bool -> Bool -> Bool
&& (Double -> Bool) -> [Double] -> Bool
forall (t :: * -> *) a. Foldable t => (a -> Bool) -> t a -> Bool
all Double -> Bool
finite [Double]
parameters) String
"invalid sensitivity step or parameters"
    let terminal :: [Double] -> [Double]
terminal [Double]
p = case Method
-> ODEOptions
-> Int
-> (Double -> [Double] -> [Double])
-> Double
-> Double
-> [Double]
-> Either MathError Trajectory
solveIVP Method
method ODEOptions
options Int
panels ([Double] -> Double -> [Double] -> [Double]
family [Double]
p) Double
start Double
end [Double]
initial of
            Left MathError
_ -> Int -> Double -> [Double]
forall a. Int -> a -> [a]
replicate ([Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
initial) (Double
0 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
0)
            Right Trajectory
trajectory -> [[Double]] -> [Double]
forall a. HasCallStack => [a] -> a
last (Trajectory -> [[Double]]
trajectoryStates Trajectory
trajectory)
    matrix <- Double
-> ([Double] -> [Double])
-> [Double]
-> Either MathError (Matrix Double)
jacobian Double
h [Double] -> [Double]
terminal [Double]
parameters
    Right (matrix)

-- | Find an autonomous equilibrium by a locally convergent Newton solve.
equilibrium :: ODEOptions -> ([Double] -> [Double]) -> [Double] -> Either MathError [Double]
equilibrium :: ODEOptions
-> ([Double] -> [Double]) -> [Double] -> Either MathError [Double]
equilibrium ODEOptions
options [Double] -> [Double]
field [Double]
guess = ODEOptions -> Either MathError ()
validateOptions ODEOptions
options Either MathError ()
-> Either MathError [Double] -> Either MathError [Double]
forall a b.
Either MathError a -> Either MathError b -> Either MathError b
forall (m :: * -> *) a b. Monad m => m a -> m b -> m b
>> ODEOptions
-> ([Double] -> [Double]) -> [Double] -> Either MathError [Double]
newton ODEOptions
options [Double] -> [Double]
field [Double]
guess

-- | Compute autonomous phase-space trajectories for independently supplied starts.
phaseTrajectories ::
    Method ->
    ODEOptions ->
    Int ->
    ([Double] -> [Double]) ->
    Double ->
    Double ->
    [[Double]] ->
    Either MathError [Trajectory]
phaseTrajectories :: Method
-> ODEOptions
-> Int
-> ([Double] -> [Double])
-> Double
-> Double
-> [[Double]]
-> Either MathError [Trajectory]
phaseTrajectories Method
method ODEOptions
options Int
panels [Double] -> [Double]
field Double
start Double
end = ([Double] -> Either MathError Trajectory)
-> [[Double]] -> Either MathError [Trajectory]
forall (t :: * -> *) (f :: * -> *) a b.
(Traversable t, Applicative f) =>
(a -> f b) -> t a -> f (t b)
forall (f :: * -> *) a b.
Applicative f =>
(a -> f b) -> [a] -> f [b]
traverse (Method
-> ODEOptions
-> Int
-> (Double -> [Double] -> [Double])
-> Double
-> Double
-> [Double]
-> Either MathError Trajectory
solveIVP Method
method ODEOptions
options Int
panels (\Double
_ -> [Double] -> [Double]
field) Double
start Double
end)

{- | Full Lyapunov spectrum by the Benettin tangent-space method: integrate
variational equations with RK4, then modified Gram-Schmidt at each sample.
Finite-difference Jacobians are used; results are finite-time estimates.
-}
lyapunovExponents ::
    ODEOptions -> Int -> Double -> ODE -> Double -> Double -> [Double] -> Either MathError [Double]
lyapunovExponents :: ODEOptions
-> Int
-> Double
-> (Double -> [Double] -> [Double])
-> Double
-> Double
-> [Double]
-> Either MathError [Double]
lyapunovExponents ODEOptions
options Int
panels Double
differenceStep Double -> [Double] -> [Double]
f Double
start Double
end [Double]
initial = do
    ODEOptions -> Either MathError ()
validateOptions ODEOptions
options
    Bool -> String -> Either MathError ()
require
        ( Int
panels Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
> Int
0
            Bool -> Bool -> Bool
&& Int
panels Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
<= ODEOptions -> Int
maximumSteps ODEOptions
options
            Bool -> Bool -> Bool
&& Double -> Bool
finite Double
differenceStep
            Bool -> Bool -> Bool
&& Double
differenceStep Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
> Double
0
            Bool -> Bool -> Bool
&& Double -> Bool
finite Double
start
            Bool -> Bool -> Bool
&& Double -> Bool
finite Double
end
            Bool -> Bool -> Bool
&& Double
end Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
> Double
start
            Bool -> Bool -> Bool
&& (Double -> Bool) -> [Double] -> Bool
forall (t :: * -> *) a. Foldable t => (a -> Bool) -> t a -> Bool
all Double -> Bool
finite [Double]
initial
        )
        String
"invalid Lyapunov arguments"
    let n :: Int
n = [Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
initial
        unit :: [Double]
unit = [if Int
i Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
j then Double
1 else Double
0 | Int
i <- [Int
0 .. Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1], Int
j <- [Int
0 .. Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1]]
        augmented :: Double -> [Double] -> [Double]
augmented Double
t [Double]
state =
            let ([Double]
y, [Double]
v) = Int -> [Double] -> ([Double], [Double])
forall a. Int -> [a] -> ([a], [a])
splitAt Int
n [Double]
state
             in case Double
-> ([Double] -> [Double])
-> [Double]
-> Either MathError (Matrix Double)
jacobian Double
differenceStep (Double -> [Double] -> [Double]
f Double
t) [Double]
y of
                    Left MathError
_ -> Int -> Double -> [Double]
forall a. Int -> a -> [a]
replicate (Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
* Int
n) (Double
0 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
0)
                    Right Matrix Double
j ->
                        let values :: [Double]
values = Matrix Double -> [Double]
forall a. Matrix a -> [a]
toList Matrix Double
j; dy :: [Double]
dy = Double -> [Double] -> [Double]
f Double
t [Double]
y
                         in if [Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
dy Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
/= Int
n
                                then Int -> Double -> [Double]
forall a. Int -> a -> [a]
replicate (Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
* Int
n) (Double
0 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
0)
                                else
                                    [Double]
dy
                                        [Double] -> [Double] -> [Double]
forall a. [a] -> [a] -> [a]
++ [ [Double] -> Double
forall a. Num a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum [[Double]
values [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! (Int
i Int -> Int -> Int
forall a. Num a => a -> a -> a
* Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
k) Double -> Double -> Double
forall a. Num a => a -> a -> a
* [Double]
v [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! (Int
k Int -> Int -> Int
forall a. Num a => a -> a -> a
* Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
c) | Int
k <- [Int
0 .. Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1]]
                                           | Int
i <- [Int
0 .. Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1]
                                           , Int
c <- [Int
0 .. Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1]
                                           ]
        orthogonalize :: t [Double] -> Either MathError ([[Double]], [Double])
orthogonalize t [Double]
columns = (([[Double]], [Double])
 -> [Double] -> Either MathError ([[Double]], [Double]))
-> ([[Double]], [Double])
-> t [Double]
-> Either MathError ([[Double]], [Double])
forall (t :: * -> *) (m :: * -> *) b a.
(Foldable t, Monad m) =>
(b -> a -> m b) -> b -> t a -> m b
foldM ([[Double]], [Double])
-> [Double] -> Either MathError ([[Double]], [Double])
orth ([], []) t [Double]
columns
        orth :: ([[Double]], [Double])
-> [Double] -> Either MathError ([[Double]], [Double])
orth ([[Double]]
known, [Double]
logs) [Double]
column = do
            let remaining :: [Double]
remaining = ([Double] -> [Double] -> [Double])
-> [Double] -> [[Double]] -> [Double]
forall b a. (b -> a -> b) -> b -> [a] -> b
forall (t :: * -> *) b a.
Foldable t =>
(b -> a -> b) -> b -> t a -> b
foldl (\[Double]
v [Double]
q -> (Double -> Double -> Double) -> [Double] -> [Double] -> [Double]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith (-) [Double]
v (Double -> [Double] -> [Double]
scale ([Double] -> Double
forall a. Num a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum ((Double -> Double -> Double) -> [Double] -> [Double] -> [Double]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith Double -> Double -> Double
forall a. Num a => a -> a -> a
(*) [Double]
q [Double]
v)) [Double]
q)) [Double]
column [[Double]]
known
                size :: Double
size = [Double] -> Double
norm [Double]
remaining
            Bool -> String -> Either MathError ()
require (Double -> Bool
finite Double
size Bool -> Bool -> Bool
&& Double
size Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
> Double
0) String
"Lyapunov tangent basis lost rank"
            ([[Double]], [Double]) -> Either MathError ([[Double]], [Double])
forall a b. b -> Either a b
Right ([[Double]]
known [[Double]] -> [[Double]] -> [[Double]]
forall a. [a] -> [a] -> [a]
++ [Double -> [Double] -> [Double]
scale (Double
1 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
size) [Double]
remaining], [Double]
logs [Double] -> [Double] -> [Double]
forall a. [a] -> [a] -> [a]
++ [Double -> Double
forall a. Floating a => a -> a
log Double
size])
        go :: Int
-> [Double] -> [Double] -> [Double] -> Either MathError [Double]
go Int
i [Double]
y [Double]
basis [Double]
sums
            | Int
i Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
panels = [Double] -> Either MathError [Double]
checked ((Double -> Double) -> [Double] -> [Double]
forall a b. (a -> b) -> [a] -> [b]
map (Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (Double
end Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
start)) [Double]
sums)
            | Bool
otherwise = do
                let t :: Double
t = Double
start Double -> Double -> Double
forall a. Num a => a -> a -> a
+ (Double
end Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
start) Double -> Double -> Double
forall a. Num a => a -> a -> a
* Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
i Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
panels
                    h :: Double
h = (Double
end Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
start) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
panels
                next <- Method
-> ODEOptions
-> (Double -> [Double] -> [Double])
-> Double
-> Double
-> [Double]
-> Either MathError [Double]
fixedStep Method
RK4 ODEOptions
options Double -> [Double] -> [Double]
augmented Double
t Double
h ([Double]
y [Double] -> [Double] -> [Double]
forall a. [a] -> [a] -> [a]
++ [Double]
basis)
                let (y1, v1) = splitAt n next; columns = [[[Double]
v1 [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! (Int
r Int -> Int -> Int
forall a. Num a => a -> a -> a
* Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
c) | Int
r <- [Int
0 .. Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1]] | Int
c <- [Int
0 .. Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1]]
                (orthonormal, logs) <- orthogonalize columns
                go (i + 1) y1 (concat (transpose orthonormal)) (add sums logs)
    if Int
n Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
0 then [Double] -> Either MathError [Double]
forall a b. b -> Either a b
Right [] else Int
-> [Double] -> [Double] -> [Double] -> Either MathError [Double]
go Int
0 [Double]
initial [Double]
unit (Int -> Double -> [Double]
forall a. Int -> a -> [a]
replicate Int
n Double
0)