1 {-# LANGUAGE ScopedTypeVariables #-}
2 {-# LANGUAGE RebindableSyntax #-}
4 -- | Numerical solution of the initial value problem,
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 )
25 -- | A single iteration of Euler's method over the interval
34 -- >>> eulers_method1 x0 y0 f h
37 eulers_method1 :: (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.
43 eulers_method1 x0 y0 f h =
46 h' = fromRational'$ toRational h
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.
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.
59 -- >>> import Algebra.Absolute (abs)
64 -- >>> let ys = eulers_method x0 xN y0 f 10000
65 -- >>> let yN = head $ reverse ys
66 -- >>> abs ((exp 1) - yN) < 1/10^3
71 eulers_method :: forall a b c. (Field.C a,
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.
82 eulers_method x0 y0 xN f n =
85 xs = partition n x0 xN
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.
92 go :: [(a,a)] -> b -> (a -> b -> b) -> [b]
94 go ((v0,v1):rest) w0 g = w1 : (go rest w1 g)
96 w1 = eulers_method1 v0 w0 g (v1 - v0)
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.
104 -- The implementation simply computes `n` from the length of the
105 -- interval and the given $h$, and then calls 'eulers_method'.
113 -- >>> let ys = eulers_method x0 xN y0 f 10
114 -- >>> let ys' = eulers_methodH x0 xN y0 f 0.1
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.
125 eulers_methodH x0 y0 xN f h =
126 eulers_method x0 y0 xN f n
129 n = floor $ (xN - x0) / h