src/ODE/IVP.hs

   1 {-# LANGUAGE ScopedTypeVariables #-}
   2 {-# LANGUAGE RebindableSyntax #-}
   3
   4 -- | Numerical solution of the initial value problem,
   5 --
   6 --     y'(x)  = f(x, y(x))
   7 --     y'(x0) = y0
   8 --
   9 --   for x in [x0, xN].
  10 --
  11
  12 module ODE.IVP
  13 where
  14
  15 import Misc (partition)
  16 import NumericPrelude hiding (abs)
  17 import qualified Algebra.Field as Field
  18 import qualified Algebra.ToInteger as ToInteger
  19 import qualified Algebra.ToRational as ToRational
  20 import qualified Algebra.RealField as RealField
  21
  22 -- | A single iteration of Euler's method over the interval
  23 --   [$x0$, $x0$+$h$].
  24 --
  25 --   Examples:
  26 --
  27 --   >>> let x0 = 0.0
  28 --   >>> let y0 = 1.0
  29 --   >>> let h  = 1.0
  30 --   >>> let f x y = y
  31 --   >>> eulers_method1 x0 y0 f h
  32 --   2.0
  33 --
  34 eulers_method1 :: (Field.C a, ToRational.C a, Field.C b)
  35                => a -- ^ x0, the initial point
  36                -> b -- ^ y0, the initial value at x0
  37                -> (a -> b -> b) -- ^ The function f(x,y)
  38                -> a -- ^ The step size h.
  39                -> b
  40 eulers_method1 x0 y0 f h =
  41   y0 +  h'*y'
  42   where
  43     h' = fromRational'$ toRational h
  44     y' = (f x0 y0)
  45
  46
  47 -- | Perform $n$ iterations of Euler's method over the interval [$x0$,
  48 --   $xN$]. The step size `h` will be calculated automatically. A list
  49 --   of y-values will be returned.
  50 --
  51 --   The explicit 'forall' in the type signature allows us to refer
  52 --   back to the type variables 'a' and 'b' in the 'where' clause.
  53 --
  54 --   Examples:
  55 --
  56 --   >>> import Algebra.Absolute (abs)
  57 --   >>> let x0 = 0.0
  58 --   >>> let xN = 1.0
  59 --   >>> let y0 = 1.0
  60 --   >>> let f x y = y
  61 --   >>> let ys = eulers_method x0 xN y0 f 10000
  62 --   >>> let yN = head $ reverse ys
  63 --   >>> abs ((exp 1) - yN) < 1/10^3
  64 --   True
  65 --   >>> head ys == y0
  66 --   True
  67 --
  68 eulers_method :: forall a b c. (Field.C a,
  69                                 ToRational.C a,
  70                                 Field.C b,
  71                                 ToInteger.C c,
  72                                 Enum c)
  73                => a -- ^ x0, the initial point
  74                -> b -- ^ y0, the initial value at x0
  75                -> a -- ^ xN, the terminal point
  76                -> (a -> b -> b) -- ^ The function f(x,y)
  77                -> c -- ^ n, the number of intervals to use.
  78                -> [b]
  79 eulers_method x0 y0 xN f n =
  80   y0 : go xs y0 f
  81   where
  82     xs = partition n x0 xN
  83
  84     -- The 'go' function actually does all the work. It takes a list
  85     -- of intervals [(v0,v1), (v1, v2)...] and peels off the first
  86     -- one. It then runs the single-step Euler's method on that
  87     -- interval, and afterwards recurses down the rest of the list.
  88     --
  89     go :: [(a,a)] -> b -> (a -> b -> b) -> [b]
  90     go [] _ _ = []
  91     go ((v0,v1):rest) w0 g = w1 : (go rest w1 g)
  92       where
  93         w1 = eulers_method1 v0 w0 g (v1 - v0)
  94
  95
  96 -- | Perform as many iterations of Euler's method over the interval
  97 --   [$x0$, $xN$] as is necessary for the given step size $h$. The
  98 --   number of subintervals `n` will be calculated automatically. A
  99 --   list of y-values will be returned.
 100 --
 101 --   The implementation simply computes `n` from the length of the
 102 --   interval and the given $h$, and then calls 'eulers_method'.
 103 --
 104 --   Examples:
 105 --
 106 --   >>> let x0 = 0.0
 107 --   >>> let xN = 1.0
 108 --   >>> let y0 = 1.0
 109 --   >>> let f x y = y
 110 --   >>> let ys = eulers_method x0 xN y0 f 10
 111 --   >>> let ys' = eulers_methodH x0 xN y0 f 0.1
 112 --   >>> ys == ys'
 113 --   True
 114 --
 115 eulers_methodH :: (RealField.C a, ToRational.C a, Field.C b)
 116                => a -- ^ x0, the initial point
 117                -> b -- ^ y0, the initial value at x0
 118                -> a -- ^ xN, the terminal point
 119                -> (a -> b -> b) -- ^ The function f(x,y)
 120                -> a -- ^ h, the step size.
 121                -> [b]
 122 eulers_methodH x0 y0 xN f h =
 123   eulers_method x0 y0 xN f n
 124   where
 125     n :: Integer
 126     n = floor $ (xN - x0) / h