diff --git a/CHANGES.md b/CHANGES.md new file mode 100644 --- /dev/null +++ b/CHANGES.md @@ -0,0 +1,1 @@ +TBD diff --git a/LICENSE b/LICENSE new file mode 100644 --- /dev/null +++ b/LICENSE @@ -0,0 +1,30 @@ +Copyright (c) 2012, Alexander V Vershilov + +All rights reserved. + +Redistribution and use in source and binary forms, with or without +modification, are permitted provided that the following conditions are met: + + * Redistributions of source code must retain the above copyright + notice, this list of conditions and the following disclaimer. + + * Redistributions in binary form must reproduce the above + copyright notice, this list of conditions and the following + disclaimer in the documentation and/or other materials provided + with the distribution. + + * Neither the name of Alexander V Vershilov nor the names of other + contributors may be used to endorse or promote products derived + from this software without specific prior written permission. + +THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS +"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT +LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR +A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT +OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, +SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT +LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, +DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY +THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT +(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE +OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. diff --git a/README.md b/README.md new file mode 100644 --- /dev/null +++ b/README.md @@ -0,0 +1,6 @@ +Small project for different ODE solvers for haskell, in particular +symplectic solvers. + +This is very experimental and will change. The Störmer-Verlet +generates a correct orbit for Jupiter but no guarantees are given for +any of the other methods. diff --git a/Setup.hs b/Setup.hs new file mode 100644 --- /dev/null +++ b/Setup.hs @@ -0,0 +1,2 @@ +import Distribution.Simple +main = defaultMain diff --git a/diagrams/src_Math_Integrators_StormerVerlet_jupiterOrbit.svg b/diagrams/src_Math_Integrators_StormerVerlet_jupiterOrbit.svg new file mode 100644 --- /dev/null +++ b/diagrams/src_Math_Integrators_StormerVerlet_jupiterOrbit.svg @@ -0,0 +1,1 @@ +1.0e125.0e110.0-5.0e11-1.0e121.0e125.0e110.0-5.0e11-1.0e12 diff --git a/diagrams/src_Math_Integrators_StormerVerlet_mySquare.svg b/diagrams/src_Math_Integrators_StormerVerlet_mySquare.svg new file mode 100644 --- /dev/null +++ b/diagrams/src_Math_Integrators_StormerVerlet_mySquare.svg @@ -0,0 +1,1 @@ + diff --git a/numeric-ode.cabal b/numeric-ode.cabal new file mode 100644 --- /dev/null +++ b/numeric-ode.cabal @@ -0,0 +1,133 @@ +name: numeric-ode +version: 0.0.0.0 +synopsis: Ode solvers +description: Some ode solvers, e.g., Störmer-Verlet +homepage: https://github.com/qnikst/numeric-ode +license: BSD3 +license-file: LICENSE +author: Alexander V Vershilov, Dominic Steinitz +maintainer: dominic@steinitz.org +copyright: Alexander V Vershilov, Dominic Steinitz +category: Math +build-type: Simple +cabal-version: >=1.18 +extra-source-files: README.md, CHANGES.md, diagrams/*.svg +extra-doc-files: diagrams/*.svg + +source-repository head + type: git + location: https://github.com/qnikst/numeric-ode + +library + default-language: Haskell2010 + hs-source-dirs: src + exposed-modules: + Math.Integrators + Math.Integrators.ExplicitEuler + Math.Integrators.ImplicitEuler + Math.Integrators.ImplicitMidpointRule + Math.Integrators.SympleticEuler + Math.Integrators.StormerVerlet + Math.Integrators.StormerVerletAlt + Math.Integrators.RK + Math.Integrators.Implicit + Math.Integrators.Internal + Math.Integrators.RK.Internal + Math.Integrators.RK.Parser + Math.Integrators.RK.Template + Math.Integrators.RK.Types + + ghc-options: -Wall + -- other-modules: + build-depends: base>=4 && < 5, + vector>=0.9 && <1.1, + parallel>=3.2 && <3.3, + parsec == 3.1.*, + template-haskell, + linear, + lens, + primitive>=0.4 && <0.7, + text, + protolude, + mwc-random, + mwc-probability, + primitive, + ad, + reflection, + tdigest, + -- chart-unit, + numhask, + foldl + + other-extensions: TypeFamilies + FlexibleContexts + BangPatterns + QuasiQuotes + +executable Kepler + hs-source-dirs: src/Examples + main-is: KeplerProblem.hs + ghc-options: + build-depends: base, + numeric-ode, + vector>=0.9 && <1.0, + vector-space>=0.8 && <0.11, + colour, + linear, + data-default-class, + diagrams-lib, + diagrams-cairo, + Chart, + Chart-cairo, + data-accessor, + diagrams-rasterific, + diagrams-lib, + JuicyPixels, + plots, + mtl + + default-language: Haskell2010 + +-- executable TestChart +-- hs-source-dirs: src/Examples +-- main-is: TestChart.hs +-- ghc-options: +-- build-depends: base >= 4.7 && < 5, +-- chart-unit, +-- protolude, +-- foldl, +-- text, +-- numhask, +-- -- for data examples +-- mwc-random, +-- mwc-probability, +-- primitive, +-- ad, +-- reflection, +-- tdigest, +-- diagrams-cairo, +-- diagrams-lib, +-- JuicyPixels +-- default-language: Haskell2010 + +-- executable TestRasterific +-- hs-source-dirs: src/Examples +-- main-is: TestRasterific.hs +-- ghc-options: +-- build-depends: base >= 4.7 && < 5, +-- chart-unit, +-- protolude, +-- foldl, +-- text, +-- numhask, +-- -- for data examples +-- mwc-random, +-- mwc-probability, +-- primitive, +-- ad, +-- reflection, +-- tdigest, +-- diagrams-rasterific, +-- diagrams-lib, +-- JuicyPixels +-- default-language: Haskell2010 diff --git a/src/Examples/KeplerProblem.hs b/src/Examples/KeplerProblem.hs new file mode 100644 --- /dev/null +++ b/src/Examples/KeplerProblem.hs @@ -0,0 +1,115 @@ +{-# LANGUAGE NegativeLiterals #-} +{-# LANGUAGE TypeFamilies #-} +{-# LANGUAGE FlexibleContexts #-} + +{-# OPTIONS_GHC -Wall #-} + +module Main (main) where + +import qualified Data.Vector as V +import Control.Monad.ST + +import Math.Integrators.StormerVerlet +import Math.Integrators + +import qualified Linear as L +import Linear.V +import Data.Maybe ( fromJust ) + +import Diagrams.Prelude +import Diagrams.Backend.CmdLine +import Diagrams.Backend.Rasterific.CmdLine + +import Control.Monad +import Control.Monad.State.Class + +import Plots + + +gConst :: Double +gConst = 6.67384e-11 + +nStepsTwoPlanets :: Int +nStepsTwoPlanets = 44 + +stepTwoPlanets :: Double +stepTwoPlanets = 24 * 60 * 60 * 100 + +sunMass, jupiterMass :: Double +sunMass = 1.9889e30 +jupiterMass = 1.8986e27 + +jupiterPerihelion :: Double +jupiterPerihelion = 7.405736e11 + +jupiterV :: [Double] +jupiterV = [-1.0965244901087316e02, -1.3710001990210707e04, 0.0] + +jupiterQ :: [Double] +jupiterQ = [negate jupiterPerihelion, 0.0, 0.0] + +sunV :: [Double] +sunV = [0.0, 0.0, 0.0] + +sunQ :: [Double] +sunQ = [0.0, 0.0, 0.0] + +tm :: V.Vector Double +tm = V.enumFromStepN 0 stepTwoPlanets nStepsTwoPlanets + +kepler :: L.V2 (L.V3 Double) -> L.V2 (L.V3 Double) +kepler (L.V2 q1 q2) = + let r = q2 L.^-^ q1 + ri = r `L.dot` r + rr = ri * (sqrt ri) + q1' = pure gConst * r / pure rr + q2' = negate q1' + q1'' = q1' * pure sunMass + q2'' = q2' * pure jupiterMass + in L.V2 q1'' q2'' + +listToV3 :: [a] -> L.V3 a +listToV3 [x, y, z] = fromV . fromJust . fromVector . V.fromList $ [x, y, z] +listToV3 xs = error $ "Only supply 3 elements not: " ++ show (length xs) + +initPQs :: L.V2 (L.V2 (L.V3 Double)) +initPQs = L.V2 (L.V2 (listToV3 jupiterV) (listToV3 sunV)) + (L.V2 (listToV3 jupiterQ) (listToV3 sunQ)) + +result1 :: V.Vector (L.V2 (L.V2 (L.V3 Double))) +result1 = runST $ integrateV (\h -> stormerVerlet2 kepler (pure h)) initPQs tm + +preMorePts :: [(Double, Double)] +preMorePts = map (\(L.V2 _ (L.V2 (L.V3 x y _z) _)) -> (x,y)) (V.toList result1) + +morePts :: [P2 Double] +morePts = map p2 $ preMorePts + +addPoint :: (Plotable (Diagram B) b, MonadState (Axis b V2 Double) m) => + Double -> (Double, Double) -> m () +addPoint o (x, y) = addPlotable' + ((circle 1e11 :: Diagram B) # + fc brown # + opacity o # + translate (r2 (x, y))) + +jSaxis :: Axis B V2 Double +jSaxis = r2Axis &~ do + addPlotable' ((circle 1e11 :: Diagram B) # fc yellow) + let l = length preMorePts + let os = [0.05,0.1..] + let ps = take (l `div` 4) [0,4..] + zipWithM_ addPoint os (map (preMorePts!!) ps) + linePlot' $ map unp2 $ take 200 morePts + +displayHeader :: FilePath -> Diagram B -> IO () +displayHeader fn = + mainRender ( DiagramOpts (Just 900) (Just 700) fn + , DiagramLoopOpts False Nothing 0 + ) + +main :: IO () +main = do + displayHeader "other/jupiter-sun-line.png" (renderAxis jSaxis # bg ivory) + putStrLn "Finished" + diff --git a/src/Math/Integrators.hs b/src/Math/Integrators.hs new file mode 100644 --- /dev/null +++ b/src/Math/Integrators.hs @@ -0,0 +1,35 @@ +-- | Math integrators if a high level module for different ODE integrators +-- This module provides high-level wrappers over different integration methods +-- +module Math.Integrators + where + +import Data.Vector (Vector,(!)) +import Data.Vector.Mutable +import Control.Monad.Primitive +import qualified Data.Vector as V + +import Math.Integrators.Internal + +{-| + Integrate ODE equation using fixed steps set by a vector, and returns a vector + of solutions corrensdonded to times that was requested. + It takes Vector of time points as a parameter and returns a vector of results + -} +integrateV :: PrimMonad m => Integrator a -- ^ Internal integrator + -> a -- ^ initial value + -> Vector Double -- ^ vector of time points + -> m (Vector a) -- ^ vector of solution +integrateV integrator initial times = do + out <- new (V.length times) + write out 0 initial + compute initial 1 out + V.unsafeFreeze out + where + compute y i out | i == V.length times = return () + | otherwise = do + let h = (times ! i) - (times ! (i-1)) + y' = integrator h y + write out i y' + compute y' (i+1) out + diff --git a/src/Math/Integrators/ExplicitEuler.hs b/src/Math/Integrators/ExplicitEuler.hs new file mode 100644 --- /dev/null +++ b/src/Math/Integrators/ExplicitEuler.hs @@ -0,0 +1,18 @@ +{-# LANGUAGE FlexibleContexts #-} +-- | +-- Module: Math.Integrators.ExplicitEuler. +-- +-- Basic integrator using Euler method. It has following properies: +-- * allows to solve systems of the first order +-- * this method is not symplectic and tends to loose energy +-- +module Math.Integrators.ExplicitEuler + where + +import Linear + +-- | Integrator of form +-- +-- \[ \Phi[h] : y_{n+1} = y_n + h f (y_n) \] +explicitEuler :: (Num (f a), Floating a, Additive f) => (f a -> f a) -> a -> f a -> f a +explicitEuler f = \h y -> y ^+^ h *^ (f y) diff --git a/src/Math/Integrators/Implicit.hs b/src/Math/Integrators/Implicit.hs new file mode 100644 --- /dev/null +++ b/src/Math/Integrators/Implicit.hs @@ -0,0 +1,64 @@ +{-# LANGUAGE FlexibleContexts #-} +-- | Helpers for implicit integration methods +-- +-- TODO: add possibility to make function to create initial value +-- TODO: add possibility to break on step +-- TODO: add possibility to add different initial value based +-- on y0, f +-- TODO: add seq-pseq to make this stuff strict +-- TODO: add Newton iterations +module Math.Integrators.Implicit + ( -- * types + ImplicitSolver + -- * solvers + , fixedPointSolver + , fixedPoint + -- * helpers + , breakNormR + , breakNormIR + ) + where + +import Linear +import Control.Lens + +-- | Implicit solver type +type ImplicitSolver a = (a -> a) -- ^ implicit method + -> (Int -> a -> a -> Bool) -- ^ breakRule + -> a -- ^ initial value + -> a -- ^ final value + +-- | Fixed point method it iterates function f until it will break "" will +-- be reached then it returns one but last iteration +-- +fixedPointSolver :: ImplicitSolver a +fixedPointSolver f break' y0 = inner 0 y0 + where + inner i y = let y' = f y + i' = i+1 + in if break' i y y' + then y' + else inner i' y' + +fixedPoint :: (a -> a) -- ^ function + -> (a -> a -> Bool) -- ^ break rule + -> a -- ^ initial value + -> a -- ^ result +fixedPoint f break' y0 = + let y1 = f y0 + in if break' y0 y1 + then y0 + else fixedPoint f break' y1 + +-- | simple break rule that will break evaluatioin when value less then Eps +breakNormR :: Double -> Double -> Bool +breakNormR eps y = abs y < eps + +-- | same as @breakNormR@ but assume that inner type is an +-- instance of InnerField, so it's possible to use innerproduct to find norm +-- N.B function uses $||v||^2 < eps$, so epsilon should be pre evaluated +breakNormIR :: (Metric f, Floating a, Ord a, Num (f a)) => f a -> a -> Bool +breakNormIR v eps = quadrance v < eps + + + diff --git a/src/Math/Integrators/ImplicitEuler.hs b/src/Math/Integrators/ImplicitEuler.hs new file mode 100644 --- /dev/null +++ b/src/Math/Integrators/ImplicitEuler.hs @@ -0,0 +1,16 @@ +module Math.Integrators.ImplicitEuler + ( implicitEuler + ) where + +import Linear + +import Math.Integrators.Implicit +import Math.Integrators.Internal + +eps :: Floating a => a +eps = 1e-14 + +implicitEuler :: (Metric f, Ord a, Additive f, Num (f a), Floating a) + => (f a -> f a) -> a -> f a -> f a +implicitEuler f = \h y -> + fixedPoint (\x -> y ^+^ (h *^ (f x))) (\x1 x2 -> breakNormIR (x1^-^x2) eps) y diff --git a/src/Math/Integrators/ImplicitMidpointRule.hs b/src/Math/Integrators/ImplicitMidpointRule.hs new file mode 100644 --- /dev/null +++ b/src/Math/Integrators/ImplicitMidpointRule.hs @@ -0,0 +1,19 @@ +{-# LANGUAGE FlexibleContexts #-} +module Math.Integrators.ImplicitMidpointRule + ( imr + ) where + +import Linear + +import Math.Integrators.Implicit +import Math.Integrators.Internal + +eps :: Floating a => a +eps = 1e-14 + +imr :: (Metric f, Num (f a), Floating a, Ord a) + => (f a -> f a) -> a -> f a -> f a +imr f = \h y -> + fixedPoint (\x -> y ^+^ h *^ ( f ( (y^+^x)^/2) )) + (\x1 x2 -> breakNormIR (x1 ^-^ x2) eps) + y diff --git a/src/Math/Integrators/Internal.hs b/src/Math/Integrators/Internal.hs new file mode 100644 --- /dev/null +++ b/src/Math/Integrators/Internal.hs @@ -0,0 +1,10 @@ +module Math.Integrators.Internal + where + +{- | Integrator function + - \Phi [h] |-> y_0 -> y_1 + -} +type Integrator a = Double -- ^ Step + -> a -- ^ Initial value + -> a -- ^ Next value + diff --git a/src/Math/Integrators/RK.hs b/src/Math/Integrators/RK.hs new file mode 100644 --- /dev/null +++ b/src/Math/Integrators/RK.hs @@ -0,0 +1,92 @@ +{-# LANGUAGE QuasiQuotes, FlexibleContexts #-} +{-# OPTIONS_GHC -Wwarn #-} +-- | Runge-Kutta module +-- TODO: add description and history notes +-- add informations about methods properties +module Math.Integrators.RK + ( -- * explicit methods +-- rk45 +-- , rk46 +-- -- * implicit methods +-- , gauss4 +-- , gauss6 +-- , lobattoIIIA4 +-- , lobattoIIIA6 +-- , lobattoIIIB4 + ) where + +-- import Linear + +-- import Math.Integrators.RK.Template +-- import Math.Integrators.RK.Types +-- import Math.Integrators.Internal +-- import Math.Integrators.Implicit + + +{- +rk45 :: (VectorSpace a, Floating (Scalar a)) => (Double -> a -> a) -> Integrator (Double,a) +rk45 = [qrk| +0 | +0.5 | 0.5 +0.5 | 0 & 0.5 +1 | 0 & 0 & 1 +- - + - - - - - - - - + | 1/6 & 2/6 & 2/6 & 1/6 +|] + +rk46 :: (VectorSpace a, Floating (Scalar a)) => (Double -> a -> a) -> Integrator (Double,a) +rk46 = [qrk| +0 | +1/3 | 1/3 +2/3 | -1/3 & 1 +1 | 1 & -1 & 1 +- - + - - - - - - - - + | 1/8 & 3/8 & 3/8 & 1/8 +|] + + +gauss4 :: (VectorSpace a, Floating (Scalar a)) => (ImplicitRkType (a,a)) -> (Double -> a -> a) -> Integrator (Double,a) +gauss4 = [qrk| +0.5 - sqrt(3)/6 | 0.25 & 0.25 - sqrt(3)/6 +0.5 + sqrt(3)/6 | 0.25 + sqrt(3)/6 & 1/4 +- - - - - - - - + - - - - - - - - - - - + | 0.5 & 0.5 +|] + +gauss6 :: (VectorSpace a, Floating (Scalar a)) => (ImplicitRkType (a,a,a)) -> (Double -> a -> a) -> Integrator (Double,a) +gauss6 = [qrk| +0.5 - sqrt(15)/10 | 5/36 & 2/9 - sqrt(15)/15 & 5/36 - sqrt(15)/30 +0.5 | 5/36 + sqrt(15)/24 & 2/9 & 5/36 - sqrt(15)/24 +0.5 + sqrt(15)/10 | 5/36 + sqrt(15)/30 & 2/9+sqrt(15)/15 & 5/36 +- - - - - - - - - + - - - - - - - - - - - - - - - - - - - - - - - + | 5/18 & 4/9 & 5/18 +|] + +lobattoIIIA4 :: (VectorSpace a, Floating (Scalar a)) => (ImplicitRkType (a,a,a)) -> (Double -> a -> a) -> Integrator (Double,a) +lobattoIIIA4 = [qrk| +0 | 0 & 0 & 0 +0.5 | 5/24 & 1/3 & -1/24 +1 | 1/6 & 2/3 & 1/6 +- - + - - - - - - - - - - + | 1/6 & 2/3 & 1/6 +|] + +lobattoIIIA6 :: (VectorSpace a, Floating (Scalar a)) => (ImplicitRkType (a,a,a,a)) -> (Double -> a -> a) -> Integrator (Double,a) +lobattoIIIA6 = [qrk| +0 | 0 & 0 & 0 & 0 +(5-sqrt(5))/10 | (11+sqrt(5))/120 & (25-sqrt(5))/120 & (25 - 13 *sqrt(5)/120) & (-1+sqrt(5))/120 +(5+sqrt(5))/10 | (11-sqrt(5))/120 & (25+13*sqrt(5))/120 & (25+sqrt(5))/120 & (-1-sqrt(5))/120 +1 | 1/12 & 5/12 & 5/12 & 1/12 +- - - - - - - - + - - - - + | 1/12 & 5/12 & 5/12 & 1/12 +|] + +lobattoIIIB4 :: (VectorSpace a, Floating (Scalar a)) => (ImplicitRkType (a,a,a)) -> (Double -> a -> a) -> Integrator (Double,a) +lobattoIIIB4 = [qrk| +0 | 1/6 & -1/6 & 0 +0.5 | 1/6 & 1/3 & 0 +1 | 1/6 & 5/6 & 0 +- - + - - - - - - - - - + | 1/6 & 2/3 & 1/6 +|] +-} diff --git a/src/Math/Integrators/RK/Internal.hs b/src/Math/Integrators/RK/Internal.hs new file mode 100644 --- /dev/null +++ b/src/Math/Integrators/RK/Internal.hs @@ -0,0 +1,14 @@ +module Math.Integrators.RK.Internal + ( MExp(..) + , isExplicit + ) + where + + +-- | Internal type that users by +data MExp = Delimeter | Row (Maybe Double,[Double]) deriving (Show,Eq) + +isExplicit :: [MExp] -> Bool +isExplicit = (all (\(i,(Row (_,x))) -> i> length x)) . (zip [1..]) . top + where + top = takeWhile (/= Delimeter) diff --git a/src/Math/Integrators/RK/Parser.hs b/src/Math/Integrators/RK/Parser.hs new file mode 100644 --- /dev/null +++ b/src/Math/Integrators/RK/Parser.hs @@ -0,0 +1,71 @@ +{-# OPTIONS_GHC -Wwarn #-} -- We need this option, because we want to remove this module in future +module Math.Integrators.RK.Parser + ( readMatrixTable + ) + where + +import Data.Maybe + +-- Parsec stuff +import Text.Parsec +import qualified Text.Parsec.Token as P +import Text.Parsec.Language (haskellDef) +import Text.Parsec.Expr +import Text.Parsec.String + +import Math.Integrators.RK.Internal + +readMatrixTable :: String -> [MExp] +readMatrixTable = + mapMaybe go . lines + where + go ('-':s) = Just Delimeter + go ('#':s) = Nothing + go (ls) | null.filter (==' ')$ ls = Nothing + | otherwise = + let (lhs,_:rhs) = span (/='|') ls + l = case filter (/= ' ') lhs of + "" -> Nothing + ls -> Just $! erun expr ls + in Just $ Row (l, map (erun expr) $! grp '&' rhs) + grp c s = case dropWhile (==c) s of + "" -> [] + s' -> if any (/=' ') w then w : grp c c'' else grp c c'' + where (w,c'') = break (==c) s' + +lexer = P.makeTokenParser haskellDef { P.reservedOpNames = ["*","/","+","-","sqrt","sin","cos"] } + +whiteSpace= P.whiteSpace lexer +lexeme = P.lexeme lexer +symbol = P.symbol lexer +float = P.float lexer +parens = P.parens lexer +natural = P.natural lexer +identifier= P.identifier lexer +reserved = P.reserved lexer +reservedOp= P.reservedOp lexer + +expr :: Parser Double +expr = buildExpressionParser table factor + "expression" +factor = parens expr + <|> try float + <|> fmap realToFrac natural + "simple expression" +table = [ [prefix "-" negate] + , [prefix "sqrt" sqrt,prefix "sin" sin,prefix "cos" cos] + , [op "*" (*) AssocLeft, op "/" (/) AssocLeft] + , [op "+" (+) AssocLeft, op "-" (-) AssocLeft] + ] + where + op s f assoc = Infix (do{ reservedOp s; return f} "operator") assoc + prefix s f = Prefix (do { reservedOp s; return f} "prefix") + + +erun :: Parser Double -> String -> Double +erun p input = erun' (do { whiteSpace ; x <- p ; eof; return x}) + where + erun' p' = case (parse p' "" input) of + Left err -> error $ "Parse error at "++(show err) + Right x -> x + diff --git a/src/Math/Integrators/RK/Template.hs b/src/Math/Integrators/RK/Template.hs new file mode 100644 --- /dev/null +++ b/src/Math/Integrators/RK/Template.hs @@ -0,0 +1,169 @@ +{-# LANGUAGE TemplateHaskell #-} +{-# OPTIONS_GHC -Wwarn #-} -- this module will be removed in future versions +module Math.Integrators.RK.Template + where + +import Data.Maybe + +import Language.Haskell.TH.Quote +import Language.Haskell.TH + +import Control.Monad + +import Math.Integrators.RK.Types +import Math.Integrators.RK.Internal +import Math.Integrators.RK.Parser + + +qrk :: QuasiQuoter +qrk = QuasiQuoter {quoteExp = x} + where + x s = rk $! readMatrixTable s + +----------------------------------------------------------- +-- List of helpers +jv = Just . VarE +jv' = Just . varE +jld = Just . LitE. {-DoublePrimL-} RationalL . toRational +ld = LitE . RationalL . toRational +plus = VarE $! mkName "+" +vplus = VarE $! mkName "^+^" +vmult = VarE $! mkName "*^" +mult = VarE $! mkName "*" +ld' = litE . RationalL . toRational +plus' = varE $! mkName "+" +vplus'= varE $! mkName "^+^" +vmult'= varE $! mkName "*^" +mult' = varE $! mkName "*" +vN s = varE (mkName s) + +foldOp op = foldl1 (\x y -> infixE (Just x) op (Just y)) + +realToFracN = varE (mkName "realToFrac") +zeroVN = varE (mkName "zeroV") +f = mkName "f" +t = mkName "t" +h = mkName "h" +y = mkName "y" +tpy = mkName "tpy" + +rk :: [MExp] -> Q Exp +rk mExp = do + let (ab,_:c:[]) = break (==Delimeter) mExp + lenA = length ab + kn <- forM [1..lenA] (\_ -> newName "k") + let kvv = zip kn ab + ks' <- forM kvv $ \(k,r) -> do + t <- rkt1 r kn + return $ ValD (VarP k) (NormalB t) [] + y' <- rkt2 c kn + if isExplicit mExp + then return $ LamE [VarP f, VarP h, TupP [VarP t,VarP y]] $ LetE ks' y' + else irk mExp + where + rkt1 (Row (Just c,ls)) ks = do + let ft = infixE (jv' t) plus' (Just $ infixE (Just $ ld' c) mult' (jv' h)) + st = if null ls + then (varE y) + else + infixE (jv' y) vplus' + (Just $ infixE (Just $ appE realToFracN (varE h)) vmult' + (Just $ foldOp vplus' $ + zipWith (\k l -> infixE (Just $ appE realToFracN (ld' l)) vmult' (jv' k)) + ks + ls + ) + ) + appE (appE (varE f) ft) st + rkt2 (Row (_,ls)) ks = do + tupE [ infixE (jv' t) plus' (jv' h) + , infixE (jv' y) vplus' + (Just $ infixE (Just $ appE realToFracN (varE h)) vmult' + (Just $ foldOp vplus' $ + zipWith (\k l -> infixE (Just $ appE realToFracN (ld' l)) vmult' (jv' k)) ks ls + ) + ) + ] + + +test = [Row (Just 1,[2,3]),Row (Just 4,[5,6]),Delimeter, Row (Nothing, [7,8])] + +irk mExp = do + fpoint' <- fpoint mExp + lamE [varP tpy, varP f, varP h, tupP [varP t,varP y]] $ + caseE (varE tpy) [match (conP (mkName "Math.Integrators.RK.Types.FixedPoint") [varP $ mkName "breakRule"]) + (normalB (fpointRun mExp)) + fpoint' + ,match (conP (mkName "Math.Integrators.RK.Types.NewtonIteration") []) (normalB (varE $ mkName "undefined")) [] + ] + +fpointRun mExp = do + let (ab,_:(Row (_,ls)):[]) = break (==Delimeter) mExp + lenA = length ab + zs <- forM [1..lenA] (\_ -> newName "z") + letE [valD (tupP $ map varP zs) + (normalB $ + appE + (appE + (appE + (varE $! mkName "Math.Integrators.Implicit.fixedPointSolver") + (varE $! mkName "method") + ) + (varE $ mkName "breakRule") + ) + (tupE $ replicate lenA zeroVN) {- TODO: give avaliability to user -} + ) + []] + (appE (varE (mkName "solution")) (tupE $ map varE zs)) + +fpoint mExp = do + let (ab,_:(Row (_,ls)):[]) = break (==Delimeter) mExp + lenA = length ab + zs <- forM [1..lenA] (\_ -> newName "z") + zs' <- forM [1..lenA] (const $ newName "z'") + + return $ + [ funD (mkName "method") [clause [tupP $ map varP zs] (normalB $ letE (map (topRow zs) $! zip zs' ab) (tupE $ map varE zs')) [] ] + , funD (mkName "solution") [clause [tupP $ map varP zs] (normalB $ solutionRow zs ls) []] + ] + where + topRow zs (x,(Row (Just c,ls))) = + valD (varP x) + (normalB $ infixE (Just $ appE realToFracN (varE h)) + vmult' + (Just $ foldOp vplus' $ + zipWith (\z l -> infixE (Just $ appE realToFracN (ld' l)) + vmult' + (Just $ appE (appE (varE f) + (infixE (jv' t) + plus' + (Just $ infixE (Just $ appE realToFracN $ varE h) mult' (Just $ ld' c)) + ) + ) + (infixE (jv' y) vplus' (jv' z)) + ) + ) + zs + ls + ) + ) + [] + topRow _ _ = error "not a row" + solutionRow zs ls = + (tupE [infixE (jv' t) plus' (jv' h) + ,infixE (jv' y) + vplus' + (Just $ infixE (Just $ appE realToFracN $ varE h) + vmult' + (Just $ foldOp vplus' + (zipWith (\z b -> infixE (Just $ appE realToFracN + (ld' b)) + vmult' + (jv' z)) + zs + ls + ) + ) + ) + ] + ) diff --git a/src/Math/Integrators/RK/Types.hs b/src/Math/Integrators/RK/Types.hs new file mode 100644 --- /dev/null +++ b/src/Math/Integrators/RK/Types.hs @@ -0,0 +1,6 @@ +module Math.Integrators.RK.Types + where + + +-- | type implicit solver +data ImplicitRkType a = FixedPoint (Int -> a -> a -> Bool) | NewtonIteration diff --git a/src/Math/Integrators/StormerVerlet.hs b/src/Math/Integrators/StormerVerlet.hs new file mode 100644 --- /dev/null +++ b/src/Math/Integrators/StormerVerlet.hs @@ -0,0 +1,203 @@ +{-# LANGUAGE FlexibleContexts #-} + +-- | +-- Module: Math.Integrators.StormerVerlet +-- +-- +-- Störmer-Verlet is an order 2 symplectic method. This means it will +-- preserve the Hamiltonian for the system the differential equations +-- describe, for example, important for modelling planetary motion; +-- the application of something like the much-loved Runge-Kutta 4th +-- order method would either model the planet spiralling toward or +-- away from the Sun! +-- +-- Here's a diagram showing the orbit of Jupiter around the Sun. +-- +-- <> +-- +-- To create this, consider the \(n\)-body problem. The Hamiltonian is +-- +-- \[ +-- {\mathbb H} = \frac{1}{2}\sum_{i=0}^n \frac{p_i^\top p_i}{m_i} - \frac{G}{2}\sum_{i=0}^n\sum_{j \neq i} \frac{m_i m_j}{\|q_i - q_j\|} +-- \] +-- +-- Apply Hamilton's equations will gives \(2n\) first order +-- equations. To use 'stormerVerlet2' this needs to be \(n\) second order +-- equations. In this case, the Lagrangian is easy +-- +-- \[ +-- {\mathcal{L}} = \frac{1}{2}\sum_{i=0}^n \frac{p_i^\top p_i}{m_i} + \frac{G}{2}\sum_{i=0}^n\sum_{j \neq i} \frac{m_i m_j}{\|q_i - q_j\|} +-- \] +-- +-- Applying Lagrange's equation +-- +-- \[ +-- \frac{\mathrm{d}}{\mathrm{d}t}\bigg(\frac{\partial{\mathcal{L}}}{\partial\dot{q}_j}\bigg) = \frac{\partial{\mathcal{L}}}{\partial{q}_j} +-- \] +-- +-- gives +-- +-- \[ +-- m_j\ddot{q}_j = G\sum_{k \neq j}m_k m_j \frac{q_k - q_j}{\|q_k - q_j\|^3} +-- \] +-- +-- For \(n = 2\) this gives +-- +-- \[ +-- \begin{aligned} +-- \ddot{q}_1 &= m_2G\frac{q_1 - q_2}{\|q_1 - q_2\|^3} \\ +-- \ddot{q}_2 &= m_1G\frac{q_2 - q_1}{\|q_2 - q_1\|^3} +-- \end{aligned} +-- \] +-- +-- > {-# LANGUAGE NegativeLiterals #-} +-- > {-# LANGUAGE TypeFamilies #-} +-- > {-# LANGUAGE FlexibleContexts #-} +-- > {-# LANGUAGE MultiParamTypeClasses #-} +-- > +-- > import qualified Data.Vector as V +-- > import Control.Monad.ST +-- > +-- > import Math.Integrators.StormerVerlet +-- > import Math.Integrators +-- > +-- > import qualified Linear as L +-- > import Linear.V +-- > import Data.Maybe ( fromJust ) +-- > +-- > import Diagrams.Prelude +-- > +-- > import Control.Monad +-- > import Control.Monad.State.Class +-- > +-- > import Plots +-- > +-- > -- First some constants describing the system +-- > +-- > gConst :: Double +-- > gConst = 6.67384e-11 +-- > +-- > nStepsTwoPlanets :: Int +-- > nStepsTwoPlanets = 44 +-- > +-- > -- A step size of 100 days! +-- > +-- > stepTwoPlanets :: Double +-- > stepTwoPlanets = 24 * 60 * 60 * 100 +-- > +-- > sunMass, jupiterMass :: Double +-- > sunMass = 1.9889e30 +-- > jupiterMass = 1.8986e27 +-- > +-- > jupiterPerihelion :: Double +-- > jupiterPerihelion = 7.405736e11 +-- > +-- > jupiterV :: L.V3 Double +-- > jupiterV = L.V3 (-1.0965244901087316e02) (-1.3710001990210707e04) 0.0 +-- > +-- > jupiterQ :: L.V3 Double +-- > jupiterQ = L.V3 (-jupiterPerihelion) 0.0 0.0 +-- > +-- > sunV :: L.V3 Double +-- > sunV = L.V3 0.0 0.0 0.0 +-- > +-- > sunQ :: L.V3 Double +-- > sunQ = L.V3 0.0 0.0 0.0 +-- > +-- > -- The right hand side of the second order differential equation system. +-- > +-- > kepler :: L.V2 (L.V3 Double) -> L.V2 (L.V3 Double) +-- > kepler (L.V2 q1 q2) = +-- > let r = q2 L.^-^ q1 +-- > ri = r `L.dot` r +-- > rr = ri * (sqrt ri) +-- > q1' = pure gConst * r / pure rr +-- > q2' = negate q1' +-- > q1'' = q1' * pure sunMass +-- > q2'' = q2' * pure jupiterMass +-- > in L.V2 q1'' q2'' +-- > +-- > -- Initial values +-- > +-- > initPQs :: L.V2 (L.V2 (L.V3 Double)) +-- > initPQs = L.V2 (L.V2 jupiterV sunV) (L.V2 jupiterQ sunQ) +-- > +-- > -- Steps at which to evolve the system +-- > +-- > tm :: V.Vector Double +-- > tm = V.enumFromStepN 0 stepTwoPlanets nStepsTwoPlanets +-- > +-- > -- The results +-- > +-- > result1 :: V.Vector (L.V2 (L.V2 (L.V3 Double))) +-- > result1 = runST $ integrateV (\h -> stormerVerlet2 kepler (pure h)) initPQs tm +-- > +-- > preMorePts :: [(Double, Double)] +-- > preMorePts = map (\(L.V2 _ (L.V2 (L.V3 x y _z) _)) -> (x,y)) (V.toList result1) +-- > +-- > morePts :: [P2 Double] +-- > morePts = map p2 $ preMorePts +-- > +-- > -- Finally plot the results +-- > +-- > addPoint :: (Plotable (Diagram B) b, MonadState (Axis b V2 Double) m) => +-- > Double -> (Double, Double) -> m () +-- > addPoint o (x, y) = addPlotable' +-- > ((circle 1e11 :: Diagram B) # +-- > fc brown # +-- > opacity o # +-- > translate (r2 (x, y))) +-- > +-- > jSaxis :: Axis B V2 Double +-- > jSaxis = r2Axis &~ do +-- > addPlotable' ((circle 1e11 :: Diagram B) # fc yellow) +-- > let l = length preMorePts +-- > let os = [0.05,0.1..] +-- > let ps = take (l `div` 4) [0,4..] +-- > zipWithM_ addPoint os (map (preMorePts!!) ps) +-- > linePlot' $ map unp2 $ take 200 morePts +-- > +-- > jupiterOrbit = renderAxis jSaxis # bg ivory +-- + +module Math.Integrators.StormerVerlet + ( stormerVerlet2H + , stormerVerlet2 + ) where + +import Linear +import Control.Lens + +-- | Störmer-Verlet integration scheme for systems of the form +-- \(\mathbb{H}(p,q) = T(p,q) + V(p,q)\) +stormerVerlet2H :: (Applicative f, Num (f a), Fractional a) => + a -- ^ Step size + -> (f a -> f a) -- ^ \(\frac{\partial H}{\partial q}\) + -> (f a -> f a) -- ^ \(\frac{\partial H}{\partial p}\) + -> V2 (f a) -- ^ Current \((p, q)\) as a 2-dimensional vector + -> V2 (f a) -- ^ New \((p, q)\) as a 2-dimensional vector +stormerVerlet2H hh nablaQ nablaP prev = V2 qNew pNew + where + h2 = hh / 2 + hhs = pure hh + hh2s = pure h2 + qsPrev = prev ^. _x + psPrev = prev ^. _y + pp2 = psPrev - hh2s * nablaQ qsPrev + qNew = qsPrev + hhs * nablaP pp2 + pNew = pp2 - hh2s * nablaQ qNew + +-- | Störmer-Verlet integration scheme for system: \(\ddot{\mathbf{q}} = f(\mathbf{q})\) +stormerVerlet2 :: (Applicative f, Num (f a), Fractional a) + => (f a -> f a) -- ^ \(f\) + -> a -- ^ Step size + -> V2 (f a) -- ^ Current \((p, q)\) as a 2-dimensional vector + -> V2 (f a) -- ^ New \((p, q)\) as a 2-dimensional vector +stormerVerlet2 f h prev = + let h' = h + h2' = 0.5 * h + p1 = prev ^. _x + pure h2' * (f (prev ^. _y)) + q' = prev ^. _y + pure h' * p1 + p' = p1 + pure h2' * (f q') + in V2 p' q' + diff --git a/src/Math/Integrators/StormerVerletAlt.hs b/src/Math/Integrators/StormerVerletAlt.hs new file mode 100644 --- /dev/null +++ b/src/Math/Integrators/StormerVerletAlt.hs @@ -0,0 +1,47 @@ +{-# OPTIONS_GHC -Wall #-} +{-# OPTIONS_GHC -fno-warn-type-defaults #-} + +module Math.Integrators.StormerVerletAlt where + +import Linear +import Control.Lens + +oneStepH98 :: (Applicative f, Num (f a), Fractional a) => + a -- ^ Step size + -> (f a -> f a) -- ^ \(\frac{\partial H}{\partial q}\) + -> (f a -> f a) -- ^ \(\frac{\partial H}{\partial p}\) + -> V2 (f a) -- ^ Current \((p, q)\) as a 2-dimensional vector + -> V2 (f a) -- ^ New \((p, q)\) as a 2-dimensional vector +oneStepH98 hh nablaQ nablaP prev = V2 qNew pNew + where + h2 = hh / 2 + hhs = pure hh + hh2s = pure h2 + qsPrev = prev ^. _x + psPrev = prev ^. _y + pp2 = psPrev - hh2s * nablaQ qsPrev + qNew = qsPrev + hhs * nablaP pp2 + pNew = pp2 - hh2s * nablaQ qNew + +-- And now can apply this to the two body problem with the following +-- derivatives of the Hamiltonian. + +nablaQ' :: V2 Double -> V2 Double +nablaQ' qs = V2 (qq1 / r) (qq2 / r) + where + qq1 = qs ^. _x + qq2 = qs ^. _y + r = (qq1 ^ 2 + qq2 ^ 2) ** (3/2) + +nablaP' :: V2 Double -> V2 Double +nablaP' ps = ps + +e, q10, q20, p10, p20 :: Double +e = 0.6 +q10 = 1 - e +q20 = 0.0 +p10 = 0.0 +p20 = sqrt ((1 + e) / (1 - e)) + +inits :: V2 (V2 Double) +inits = V2 (V2 q10 q20) (V2 p10 p20) diff --git a/src/Math/Integrators/SympleticEuler.hs b/src/Math/Integrators/SympleticEuler.hs new file mode 100644 --- /dev/null +++ b/src/Math/Integrators/SympleticEuler.hs @@ -0,0 +1,33 @@ +{-# LANGUAGE FlexibleContexts #-} +module Math.Integrators.SympleticEuler + where + +import Linear +import Control.Lens + +import Math.Integrators.Implicit + +eps :: Floating a => a +eps = 1e-10 + +sympleticEuler1 :: (Metric f, Num (f a), Floating a, Ord a) + => (f a -> f a -> f a) + -> (f a -> f a -> f a) + -> a -- ^ Step size + -> V2 (f a) -- ^ Current \((p,q)\) as a 2-dimentional vector + -> V2 (f a) -- ^ New \((p, q)\) as a 2-dimetional vector +sympleticEuler1 f g = \h prev -> + -- explicit coordinate + let u' = (prev^._x) ^+^ h *^ (f (prev^._x) v') + -- implicit coordinate + v' = fixedPoint (\x -> (prev^._y) ^+^ h *^ (g (prev^._x) x)) + (\x1 x2 -> breakNormIR (x1^-^x2) eps) (prev^._x) + in V2 u' v' + +{- +sEuler2 :: ((a->a->a),(a->a->a)) -> Double -> (a,a) -> (a,a) +sEuler2 (a,b) h (u,v) = + let u' = u + ( h * (a u' v) ) + v' = v + ( h * (b u' v) ) + in (u',v') +-}