1 {-# LANGUAGE ScopedTypeVariables #-}
3 -- | Numerical solution of the initial value problem,
14 import Misc (partition)
17 -- | A single iteration of Euler's method over the interval
26 -- >>> eulers_method1 x0 y0 f h
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.
35 eulers_method1 x0 y0 f h =
38 h' = fromRational $ toRational h
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.
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.
55 -- >>> let ys = eulers_method x0 xN y0 f 10000
56 -- >>> let yN = head $ reverse ys
57 -- >>> abs ((exp 1) - yN) < 1/10^3
62 eulers_method :: forall a b c. (RealFrac a, RealFrac b, Integral c)
63 => a -- ^ x0, the initial point
64 -> b -- ^ y0, the initial value at x0
65 -> a -- ^ xN, the terminal point
66 -> (a -> b -> b) -- ^ The function f(x,y)
67 -> c -- ^ n, the number of intervals to use.
69 eulers_method x0 y0 xN f n =
72 xs = partition n x0 xN
74 -- The 'go' function actually does all the work. It takes a list
75 -- of intervals [(v0,v1), (v1, v2)...] and peels off the first
76 -- one. It then runs the single-step Euler's method on that
77 -- interval, and afterwards recurses down the rest of the list.
79 go :: [(a,a)] -> b -> (a -> b -> b) -> [b]
81 go ((v0,v1):rest) w0 g = w1 : (go rest w1 g)
83 w1 = eulers_method1 v0 w0 g (v1 - v0)
86 -- | Perform as many iterations of Euler's method over the interval
87 -- [$x0$, $xN$] as is necessary for the given step size $h$. The
88 -- number of subintervals `n` will be calculated automatically. A
89 -- list of y-values will be returned.
91 -- The implementation simply computes `n` from the length of the
92 -- interval and the given $h$, and then calls 'eulers_method'.
100 -- >>> let ys = eulers_method x0 xN y0 f 10
101 -- >>> let ys' = eulers_methodH x0 xN y0 f 0.1
105 eulers_methodH :: (RealFrac a, RealFrac b)
106 => a -- ^ x0, the initial point
107 -> b -- ^ y0, the initial value at x0
108 -> a -- ^ xN, the terminal point
109 -> (a -> b -> b) -- ^ The function f(x,y)
110 -> a -- ^ h, the step size.
112 eulers_methodH x0 y0 xN f h =
113 eulers_method x0 y0 xN f n
116 n = floor $ (xN - x0) / h