diff options
author | idontgetoutmuch <dominic@steinitz.org> | 2018-04-22 03:17:40 -0700 |
---|---|---|
committer | GitHub <noreply@github.com> | 2018-04-22 03:17:40 -0700 |
commit | df9980a1bc95cbe96b6766dc232c4abc8ca28c38 (patch) | |
tree | 72eae1364abfd4ef4b022f94cb126c189f5cde85 /examples/sundials.hs | |
parent | 94db29b5d25f28708c1d430554f0c337c26172df (diff) | |
parent | 492b7e13c47bc7e42583c6a1754bdb2bb445b94f (diff) |
Merge pull request #264 from idontgetoutmuch/sundials-clean
Add Sundials as an ODE Engine
Diffstat (limited to 'examples/sundials.hs')
-rw-r--r-- | examples/sundials.hs | 56 |
1 files changed, 56 insertions, 0 deletions
diff --git a/examples/sundials.hs b/examples/sundials.hs new file mode 100644 index 0000000..99f662d --- /dev/null +++ b/examples/sundials.hs | |||
@@ -0,0 +1,56 @@ | |||
1 | {-# OPTIONS_GHC -fno-warn-missing-signatures #-} | ||
2 | |||
3 | {-# LANGUAGE ViewPatterns #-} | ||
4 | |||
5 | import Numeric.Sundials.ARKode.ODE | ||
6 | import Numeric.LinearAlgebra | ||
7 | import Graphics.Plot | ||
8 | |||
9 | vanderpol mu = do | ||
10 | let xdot nu _t [x,v] = [v, -x + nu * v * (1-x*x)] | ||
11 | xdot _ _ _ = error "vanderpol RHS not defined" | ||
12 | ts = linspace 1000 (0,50) | ||
13 | sol = toColumns $ odeSolve (xdot mu) [1,0] ts | ||
14 | mplot (ts : sol) | ||
15 | mplot sol | ||
16 | |||
17 | |||
18 | harmonic w d = do | ||
19 | let xdot u dd _t [x,v] = [v, a*x + b*v] where a = -u*u; b = -2*dd*u | ||
20 | xdot _ _ _ _ = error "harmonic RHS not defined" | ||
21 | ts = linspace 100 (0,20) | ||
22 | sol = odeSolve (xdot w d) [1,0] ts | ||
23 | mplot (ts : toColumns sol) | ||
24 | |||
25 | |||
26 | kepler v a = mplot (take 2 $ toColumns sol) where | ||
27 | xdot _t [x,y,vx,vy] = [vx,vy,x*k,y*k] | ||
28 | where g=1 | ||
29 | k=(-g)*(x*x+y*y)**(-1.5) | ||
30 | xdot _ _ = error "kepler RHS not defined" | ||
31 | ts = linspace 100 (0,30) | ||
32 | sol = odeSolve xdot [4, 0, v * cos (a*degree), v * sin (a*degree)] ts | ||
33 | degree = pi/180 | ||
34 | |||
35 | |||
36 | main = do | ||
37 | vanderpol 2 | ||
38 | harmonic 1 0 | ||
39 | harmonic 1 0.1 | ||
40 | kepler 0.3 60 | ||
41 | kepler 0.4 70 | ||
42 | vanderpol' 2 | ||
43 | |||
44 | -- example of odeSolveV with jacobian | ||
45 | vanderpol' mu = do | ||
46 | let xdot nu _t (toList->[x,v]) = fromList [v, -x + nu * v * (1-x*x)] | ||
47 | xdot _ _ _ = error "vanderpol' RHS not defined" | ||
48 | jac _ (toList->[x,v]) = (2><2) [ 0 , 1 | ||
49 | , -1-2*x*v*mu, mu*(1-x**2) ] | ||
50 | jac _ _ = error "vanderpol' Jacobian not defined" | ||
51 | ts = linspace 1000 (0,50) | ||
52 | hi = pure $ (ts!1 - ts!0) / 100.0 | ||
53 | sol = toColumns $ odeSolveV (SDIRK_5_3_4 jac) hi 1E-8 1E-8 (xdot mu) (fromList [1,0]) ts | ||
54 | mplot sol | ||
55 | |||
56 | |||