Fix an off-by-one error in Euler's method.

[numerical-analysis.git] / src / ODE / IVP.hs

diff --git a/src/ODE/IVP.hs b/src/ODE/IVP.hs
index 381ff8a9c98099519dae2e08454fa8f05260127c..233e5f872a97f29e4e46dd3a34bc582484ac46ce 100644 (file)
--- a/src/ODE/IVP.hs
+++ b/src/ODE/IVP.hs
@@ -56,6 +56,8 @@ eulers_method1 x0 y0 f h =
 --   >>> 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
@@ -65,7 +67,7 @@ eulers_method :: forall a b c. (RealFrac a, RealFrac b, Integral c)
                -> 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
 
@@ -78,3 +80,35 @@ eulers_method x0 y0 xN f n =
     go ((x0,x1):rest) y0 f = y1 : (go rest y1 f)
       where
         y1 = eulers_method1 x0 y0 f (x1 - x0)
+
+
+-- | 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 = floor $ (xN - x0) / h