1 {-# LANGUAGE ScopedTypeVariables #-}
2 {-# LANGUAGE RebindableSyntax #-}
4 -- | Numerical solution of the initial value problem,
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
22 -- | A single iteration of Euler's method over the interval
31 -- >>> eulers_method1 x0 y0 f h
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.
40 eulers_method1 x0 y0 f h =
43 h' = fromRational'$ toRational h
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.
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.
56 -- >>> import Algebra.Absolute (abs)
61 -- >>> let ys = eulers_method x0 xN y0 f 10000
62 -- >>> let yN = head $ reverse ys
63 -- >>> abs ((exp 1) - yN) < 1/10^3
68 eulers_method :: forall a b c. (Field.C a,
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.
79 eulers_method x0 y0 xN f n =
82 xs = partition n x0 xN
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.
89 go :: [(a,a)] -> b -> (a -> b -> b) -> [b]
91 go ((v0,v1):rest) w0 g = w1 : (go rest w1 g)
93 w1 = eulers_method1 v0 w0 g (v1 - v0)
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.
101 -- The implementation simply computes `n` from the length of the
102 -- interval and the given $h$, and then calls 'eulers_method'.
110 -- >>> let ys = eulers_method x0 xN y0 f 10
111 -- >>> let ys' = eulers_methodH x0 xN y0 f 0.1
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.
122 eulers_methodH x0 y0 xN f h =
123 eulers_method x0 y0 xN f n
126 n = floor $ (xN - x0) / h