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