src/ODE/IVP.hs

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