From cb09d0e99ae4f10cd2b3f3ac667df83946a9744d Mon Sep 17 00:00:00 2001 From: Dominic Steinitz Date: Sun, 30 Jun 2019 08:58:16 +0100 Subject: Remove sundials as it has its own repo now. Fix * https://github.com/haskell-numerics/hmatrix/issues/304 * https://github.com/haskell-numerics/hmatrix/issues/302 * https://github.com/haskell-numerics/hmatrix/issues/290 --- packages/sundials/src/Main.hs | 186 ------------------------------------------ 1 file changed, 186 deletions(-) delete mode 100644 packages/sundials/src/Main.hs (limited to 'packages/sundials/src/Main.hs') diff --git a/packages/sundials/src/Main.hs b/packages/sundials/src/Main.hs deleted file mode 100644 index 16c21c5..0000000 --- a/packages/sundials/src/Main.hs +++ /dev/null @@ -1,186 +0,0 @@ -{-# OPTIONS_GHC -Wall #-} - -import qualified Numeric.Sundials.ARKode.ODE as ARK -import qualified Numeric.Sundials.CVode.ODE as CV -import Numeric.LinearAlgebra - -import Plots as P -import qualified Diagrams.Prelude as D -import Diagrams.Backend.Rasterific - -import Control.Lens - -import Test.Hspec - - -lorenz :: Double -> [Double] -> [Double] -lorenz _t u = [ sigma * (y - x) - , x * (rho - z) - y - , x * y - beta * z - ] - where - rho = 28.0 - sigma = 10.0 - beta = 8.0 / 3.0 - x = u !! 0 - y = u !! 1 - z = u !! 2 - -_lorenzJac :: Double -> Vector Double -> Matrix Double -_lorenzJac _t u = (3><3) [ (-sigma), rho - z, y - , sigma , -1.0 , x - , 0.0 , (-x) , (-beta) - ] - where - rho = 28.0 - sigma = 10.0 - beta = 8.0 / 3.0 - x = u ! 0 - y = u ! 1 - z = u ! 2 - -brusselator :: Double -> [Double] -> [Double] -brusselator _t x = [ a - (w + 1) * u + v * u * u - , w * u - v * u * u - , (b - w) / eps - w * u - ] - where - a = 1.0 - b = 3.5 - eps = 5.0e-6 - u = x !! 0 - v = x !! 1 - w = x !! 2 - -_brussJac :: Double -> Vector Double -> Matrix Double -_brussJac _t x = (3><3) [ (-(w + 1.0)) + 2.0 * u * v, w - 2.0 * u * v, (-w) - , u * u , (-(u * u)) , 0.0 - , (-u) , u , (-1.0) / eps - u - ] - where - y = toList x - u = y !! 0 - v = y !! 1 - w = y !! 2 - eps = 5.0e-6 - -stiffish :: Double -> [Double] -> [Double] -stiffish t v = [ lamda * u + 1.0 / (1.0 + t * t) - lamda * atan t ] - where - lamda = -100.0 - u = v !! 0 - -stiffishV :: Double -> Vector Double -> Vector Double -stiffishV t v = fromList [ lamda * u + 1.0 / (1.0 + t * t) - lamda * atan t ] - where - lamda = -100.0 - u = v ! 0 - -_stiffJac :: Double -> Vector Double -> Matrix Double -_stiffJac _t _v = (1><1) [ lamda ] - where - lamda = -100.0 - -predatorPrey :: Double -> [Double] -> [Double] -predatorPrey _t v = [ x * a - b * x * y - , d * x * y - c * y - e * y * z - , (-f) * z + g * y * z - ] - where - x = v!!0 - y = v!!1 - z = v!!2 - a = 1.0 - b = 1.0 - c = 1.0 - d = 1.0 - e = 1.0 - f = 1.0 - g = 1.0 - -lSaxis :: [[Double]] -> P.Axis B D.V2 Double -lSaxis xs = P.r2Axis &~ do - let ts = xs!!0 - us = xs!!1 - vs = xs!!2 - ws = xs!!3 - P.linePlot' $ zip ts us - P.linePlot' $ zip ts vs - P.linePlot' $ zip ts ws - -kSaxis :: [(Double, Double)] -> P.Axis B D.V2 Double -kSaxis xs = P.r2Axis &~ do - P.linePlot' xs - -main :: IO () -main = do - - let res1 = ARK.odeSolve brusselator [1.2, 3.1, 3.0] (fromList [0.0, 0.1 .. 10.0]) - renderRasterific "diagrams/brusselator.png" - (D.dims2D 500.0 500.0) - (renderAxis $ lSaxis $ [0.0, 0.1 .. 10.0]:(toLists $ tr res1)) - - let res1a = ARK.odeSolve brusselator [1.2, 3.1, 3.0] (fromList [0.0, 0.1 .. 10.0]) - renderRasterific "diagrams/brusselatorA.png" - (D.dims2D 500.0 500.0) - (renderAxis $ lSaxis $ [0.0, 0.1 .. 10.0]:(toLists $ tr res1a)) - - let res2 = ARK.odeSolve stiffish [0.0] (fromList [0.0, 0.1 .. 10.0]) - renderRasterific "diagrams/stiffish.png" - (D.dims2D 500.0 500.0) - (renderAxis $ kSaxis $ zip [0.0, 0.1 .. 10.0] (concat $ toLists res2)) - - let res2a = ARK.odeSolveV (ARK.SDIRK_5_3_4') Nothing 1e-6 1e-10 stiffishV (fromList [0.0]) (fromList [0.0, 0.1 .. 10.0]) - - let res2b = ARK.odeSolveV (ARK.TRBDF2_3_3_2') Nothing 1e-6 1e-10 stiffishV (fromList [0.0]) (fromList [0.0, 0.1 .. 10.0]) - - let maxDiffA = maximum $ map abs $ - zipWith (-) ((toLists $ tr res2a)!!0) ((toLists $ tr res2b)!!0) - - let res2c = CV.odeSolveV (CV.BDF) Nothing 1e-6 1e-10 stiffishV (fromList [0.0]) (fromList [0.0, 0.1 .. 10.0]) - - let maxDiffB = maximum $ map abs $ - zipWith (-) ((toLists $ tr res2a)!!0) ((toLists $ tr res2c)!!0) - - let maxDiffC = maximum $ map abs $ - zipWith (-) ((toLists $ tr res2b)!!0) ((toLists $ tr res2c)!!0) - - let res3 = ARK.odeSolve lorenz [-5.0, -5.0, 1.0] (fromList [0.0, 0.01 .. 10.0]) - - renderRasterific "diagrams/lorenz.png" - (D.dims2D 500.0 500.0) - (renderAxis $ kSaxis $ zip ((toLists $ tr res3)!!0) ((toLists $ tr res3)!!1)) - - renderRasterific "diagrams/lorenz1.png" - (D.dims2D 500.0 500.0) - (renderAxis $ kSaxis $ zip ((toLists $ tr res3)!!0) ((toLists $ tr res3)!!2)) - - renderRasterific "diagrams/lorenz2.png" - (D.dims2D 500.0 500.0) - (renderAxis $ kSaxis $ zip ((toLists $ tr res3)!!1) ((toLists $ tr res3)!!2)) - - let res4 = CV.odeSolve predatorPrey [0.5, 1.0, 2.0] (fromList [0.0, 0.01 .. 10.0]) - - renderRasterific "diagrams/predatorPrey.png" - (D.dims2D 500.0 500.0) - (renderAxis $ kSaxis $ zip ((toLists $ tr res4)!!0) ((toLists $ tr res4)!!1)) - - renderRasterific "diagrams/predatorPrey1.png" - (D.dims2D 500.0 500.0) - (renderAxis $ kSaxis $ zip ((toLists $ tr res4)!!0) ((toLists $ tr res4)!!2)) - - renderRasterific "diagrams/predatorPrey2.png" - (D.dims2D 500.0 500.0) - (renderAxis $ kSaxis $ zip ((toLists $ tr res4)!!1) ((toLists $ tr res4)!!2)) - - let res4a = ARK.odeSolve predatorPrey [0.5, 1.0, 2.0] (fromList [0.0, 0.01 .. 10.0]) - - let maxDiffPpA = maximum $ map abs $ - zipWith (-) ((toLists $ tr res4)!!0) ((toLists $ tr res4a)!!0) - - hspec $ describe "Compare results" $ do - it "for SDIRK_5_3_4' and TRBDF2_3_3_2'" $ maxDiffA < 1.0e-6 - it "for SDIRK_5_3_4' and BDF" $ maxDiffB < 1.0e-6 - it "for TRBDF2_3_3_2' and BDF" $ maxDiffC < 1.0e-6 - it "for CV and ARK for the Predator Prey model" $ maxDiffPpA < 1.0e-3 - -- cgit v1.2.3