eulers_method1 x0 y0 f h =
y0 + h'*y'
where
- h' = fromRational $ toRational h
+ h' = realToFrac h
y' = (f x0 y0)
-- >>> let yN = head $ reverse ys
-- >>> abs ((exp 1) - yN) < 1/10^3
-- True
+-- >>> head ys == y0
+-- True
--
eulers_method :: forall a b c. (RealFrac a, RealFrac b, Integral c)
=> a -- ^ x0, the initial point
-> c -- ^ n, the number of intervals to use.
-> [b]
eulers_method x0 y0 xN f n =
- go xs y0 f
+ y0 : go xs y0 f
where
xs = partition n x0 xN
-- The 'go' function actually does all the work. It takes a list
- -- of intervals [(x0,x1), (x1, x2)...] and peels off the first
+ -- of intervals [(v0,v1), (v1, v2)...] and peels off the first
-- one. It then runs the single-step Euler's method on that
-- interval, and afterwards recurses down the rest of the list.
+ --
go :: [(a,a)] -> b -> (a -> b -> b) -> [b]
go [] _ _ = []
- go ((x0,x1):rest) y0 f = y1 : (go rest y1 f)
+ go ((v0,v1):rest) w0 g = w1 : (go rest w1 g)
where
- y1 = eulers_method1 x0 y0 f (x1 - x0)
+ w1 = eulers_method1 v0 w0 g (v1 - v0)
+
+
+-- | Perform as many iterations of Euler's method over the interval
+-- [$x0$, $xN$] as is necessary for the given step size $h$. The
+-- number of subintervals `n` will be calculated automatically. A
+-- list of y-values will be returned.
+--
+-- The implementation simply computes `n` from the length of the
+-- interval and the given $h$, and then calls 'eulers_method'.
+--
+-- Examples:
+--
+-- >>> let x0 = 0.0
+-- >>> let xN = 1.0
+-- >>> let y0 = 1.0
+-- >>> let f x y = y
+-- >>> let ys = eulers_method x0 xN y0 f 10
+-- >>> let ys' = eulers_methodH x0 xN y0 f 0.1
+-- >>> ys == ys'
+-- True
+--
+eulers_methodH :: (RealFrac a, RealFrac b)
+ => a -- ^ x0, the initial point
+ -> b -- ^ y0, the initial value at x0
+ -> a -- ^ xN, the terminal point
+ -> (a -> b -> b) -- ^ The function f(x,y)
+ -> a -- ^ h, the step size.
+ -> [b]
+eulers_methodH x0 y0 xN f h =
+ eulers_method x0 y0 xN f n
+ where
+ n :: Integer
+ n = floor $ (xN - x0) / h