\documentclass{article}
%include polycode.fmt
\usepackage{fontspec}
\usepackage{amsmath}
\usepackage{unicode-math}
\usepackage{lualatex-math}
\setmainfont{latinmodern-math.otf}
\setmathfont{latinmodern-math.otf}
\usepackage{verbatim}
\author{Sophie Taylor}
\title{haskell-clifford: A Haskell Clifford algebra dynamics library}
\begin{document}
This is the numeric integration portion of the library.
\begin{code}
\end{code}
%if False
\begin{code}
\end{code}
%endif
\begin{code}
module Numeric.Clifford.NumericIntegration where
import Numeric.Clifford.Multivector as MV
import NumericPrelude hiding (iterate, head, map, tail, reverse, scanl, zipWith, drop, (++), filter, null, length, foldr, foldl1, zip, foldl, concat, (!!), concatMap,any, repeat, replicate, elem, replicate)
import Algebra.Absolute
import Algebra.Algebraic
import Algebra.Additive hiding (elementAdd, elementSub)
import Algebra.Ring
import Data.Monoid
import Algebra.ToInteger
import Algebra.Module
import Algebra.Field
import Data.List.Stream
import Numeric.Natural
import qualified Data.Vector as V
import NumericPrelude.Numeric (sum)
import qualified NumericPrelude.Numeric as NPN
import Test.QuickCheck
import Math.Sequence.Converge (convergeBy)
import Number.Ratio hiding (scale)
import Algebra.ToRational
import Control.Lens hiding (indices)
import Control.Exception (assert)
import Data.Maybe
import GHC.TypeLits
import Data.DeriveTH
import Debug.Trace
elementAdd = zipWith (+)
elementScale = zipWith (*>)
a `elementSub` b = zipWith () a b
a `elementMul` b = zipWith (*) a b
rk4ClassicalTableau = ButcherTableau [[0,0,0,0],[0.5,0,0,0],[0,0.5,0,0],[0,0,1,0]] [1.0 NPN./6,1.0 NPN./3, 1.0 NPN./3, 1.0 NPN./6] [0, 0.5, 0.5, 1]
implicitEulerTableau = ButcherTableau [[1.0::NPN.Double]] [1] [1]
sumVector = sumList . V.toList
systemBroydensMethod f x0 x1 = map fst $ update (x1,ident) x0 where
update (xm1,jm1) xm2 | zero `elem` dx = [(xm1,undefined)]
| otherwise = if x == xm1 then [(x,undefined)] else (x,j) : update (x,j) xm1 where
x = xm1 `elementSub` ( (fm1 `elementMul` dx) `elementMul` ody)
j = undefined
fm1 = f xm1
fm2 = f xm2
dx = elementSub xm1 xm2
dy = elementSub fm1 fm2
ody = map inverse dy
ident = undefined
convergeList ::(Show f, Ord f) => [[f]] -> [f]
convergeList = converge
showOutput name x = trace ("output of " ++ name ++" is " ++ show x) x
convergeTolLists :: (Ord f, Algebra.Absolute.C f, Algebra.Algebraic.C f, Show f, SingI d)
=> f -> [[Multivector d f]] -> [Multivector d f]
convergeTolLists t [] = error "converge: empty list"
convergeTolLists t xs = fromMaybe empty (convergeBy check Just xs)
where
empty = error "converge: error in impl"
check (a:b:c:_)
| (trace ("Converging at " ++ show a) a) == b = Just b
| a == c = Just c
| ((showOutput ("convergence check with tolerance " ++ show t) $
magnitude (sumList $ (zipWith (\x y -> NPN.abs (xy)) b c))) <= t) = showOutput ("convergence with tolerance "++ show t )$ Just c
check _ = Nothing
type RKStepper (d::Nat) t stateType =
(Ord t, Show t, Algebra.Module.C t (Multivector d t), Algebra.Field.C t) =>
t -> (t -> stateType -> stateType) ->
([Multivector d t] -> stateType) ->
(stateType ->[Multivector d t]) ->
(t,stateType) -> (t,stateType)
data ButcherTableau f = ButcherTableau {_tableauA :: [[f]], _tableauB :: [f], _tableauC :: [f]}
makeLenses ''ButcherTableau
type ConvergerFunction f = forall (d::Nat) f . [[Multivector d f]] -> [Multivector d f]
type AdaptiveStepSizeFunction f state = f -> state -> f
data RKAttribute f state = Explicit
| HamiltonianFunction
| AdaptiveStepSize {sigma :: AdaptiveStepSizeFunction f state}
| ConvergenceTolerance {epsilon :: f}
| ConvergenceFunction {converger :: ConvergerFunction f }
| RootSolver
| UseAutomaticDifferentiationForRootSolver
| StartingGuessMethod
$( derive makeIs ''RKAttribute)
genericRKMethod :: forall (d::Nat) t stateType .
( Ord t, Show t, Algebra.Module.C t (Multivector d t),Algebra.Absolute.C t, Algebra.Algebraic.C t, SingI d)
=> ButcherTableau t -> [RKAttribute t stateType] -> RKStepper d t stateType
genericRKMethod tableau attributes = rkMethodImplicitFixedPoint where
s :: Int
s = length (_tableauC tableau)
c :: Int -> t
c n = l !! (n1) where
l = _tableauC tableau
a :: Int -> [t]
a n = (l !! (n1)) & filter (/= zero) where
l = _tableauA tableau
b :: Int -> t
b i = l !! (i 1) where
l = _tableauB tableau
sumListOfLists :: [[Multivector d t]] -> [Multivector d t]
sumListOfLists = map sumList . transpose
converger :: [[Multivector d t]] -> [Multivector d t]
converger = case find (\x -> isConvergenceTolerance x || isConvergenceFunction x) attributes of
Just (ConvergenceFunction conv) -> conv
Just (ConvergenceTolerance tol) -> convergeTolLists (trace ("Convergence tolerance set to " ++ show tol)tol)
Nothing -> trace "No convergence tolerance specified, defaulting to equality" convergeList
stepSizeAdapter :: AdaptiveStepSizeFunction t stateType
stepSizeAdapter = case find isAdaptiveStepSize attributes of
Just (AdaptiveStepSize sigma) -> sigma
Nothing -> (\_ _ -> one)
rkMethodImplicitFixedPoint :: RKStepper d t stateType
rkMethodImplicitFixedPoint h f project unproject (time, state) =
(time + (stepSizeAdapter time state)*h*(c s), newState) where
zi :: Int -> [Multivector d t]
zi i = (\out -> trace ("initialGuess is " ++ show initialGuess++" whereas the final one is " ++ show out) out) $
assert (i <= s && i>= 1) $ converger $ iterate (zkp1 i) initialGuess where
initialGuess :: [Multivector d t]
initialGuess = if i == 1 || null (zi (i1)) then map (h'*>) $ unproject $ f guessTime state else zi (i1)
adaptiveStepSizeFraction :: t
adaptiveStepSizeFraction = stepSizeAdapter time state
h' :: t
h' = adaptiveStepSizeFraction * h * (c i)
guessTime :: t
guessTime = time + h'
zkp1 :: NPN.Int -> [Multivector d t] -> [Multivector d t]
zkp1 i zk = map (h*>) (sumOfJs i zk) where
sumOfJs :: Int -> [Multivector d t] -> [Multivector d t]
sumOfJs i zk = sumListOfLists $ map (scaledByAij zk) (a i) where
scaledByAij :: [Multivector d t] -> t -> [Multivector d t]
scaledByAij guess a = map (a*>) $ evalDerivatives guessTime $ elementAdd state' guess
state' = unproject state :: [Multivector d t]
newState :: stateType
newState = project $ elementAdd state' (assert (not $ null dy) dy)
dy = sumListOfLists [map ((b i) *>) (zi i) | i <- [1..s]] :: [Multivector d t]
evalDerivatives :: t -> [Multivector d t] -> [Multivector d t]
evalDerivatives time stateAtTime= unproject $ (f time) $ project stateAtTime
\end{code}
Look at creating an exponential integrator: https://en.wikipedia.org/wiki/Exponential_integrators
\bibliographystyle{IEEEtran}
\bibliography{biblio.bib}
\end{document}