From: Michael Orlitzky <michael@orlitzky.com>
Date: Thu, 6 Dec 2012 02:14:15 +0000 (-0500)
Subject: Add the ODE.IVP module with Euler's method.
X-Git-Url: https://gitweb.michael.orlitzky.com/?a=commitdiff_plain;h=ddad82e5f42a7d38d542fbb77e981dc7309089f4;p=numerical-analysis.git

Add the ODE.IVP module with Euler's method.
Clean up .ghci a little bit and add the new module.
---

diff --git a/.ghci b/.ghci
index a04730c..c9819c7 100644
--- a/.ghci
+++ b/.ghci
@@ -2,15 +2,26 @@
 :set -isrc
 
 -- Load everything.
-:l src/Roots/Simple.hs src/Aliases.hs src/TwoTuple.hs src/Integration/Simpson.hs src/Integration/Trapezoid.hs
+:{
+:load src/Aliases.hs
+  src/Integration/Simpson.hs
+  src/Integration/Trapezoid.hs
+  src/Misc.hs
+  src/ODE/IVP.hs
+  src/Roots/Simple.hs
+  src/TwoTuple.hs
+:}
 
 -- Just for convenience.
 import Data.Number.BigFloat
 
 import Aliases
-import TwoTuple
 import Integration.Simpson
 import Integration.Trapezoid
+import Misc
+import ODE.IVP
+import Roots.Simple
+import TwoTuple
 
 -- Use a calmer prompt.
 :set prompt "numerical-analysis> "
diff --git a/src/ODE/IVP.hs b/src/ODE/IVP.hs
new file mode 100644
index 0000000..381ff8a
--- /dev/null
+++ b/src/ODE/IVP.hs
@@ -0,0 +1,80 @@
+{-# LANGUAGE ScopedTypeVariables #-}
+
+-- | Numerical solution of the initial value problem,
+--
+--     y'(x)  = f(x, y(x))
+--     y'(x0) = y0
+--
+--   for x in [x0, xN].
+--
+
+module ODE.IVP
+where
+
+import Misc (partition)
+
+
+-- | A single iteration of Euler's method over the interval
+--   [$x0$, $x0$+$h$].
+--
+--   Examples:
+--
+--   >>> let x0 = 0.0
+--   >>> let y0 = 1.0
+--   >>> let h  = 1.0
+--   >>> let f x y = y
+--   >>> eulers_method1 x0 y0 f h
+--   2.0
+--
+eulers_method1 :: (RealFrac a, RealFrac b)
+               => a -- ^ x0, the initial point
+               -> b -- ^ y0, the initial value at x0
+               -> (a -> b -> b) -- ^ The function f(x,y)
+               -> a -- ^ The step size h.
+               -> b
+eulers_method1 x0 y0 f h =
+  y0 +  h'*y'
+  where
+    h' = fromRational $ toRational h
+    y' = (f x0 y0)
+
+
+-- | Perform $n$ iterations of Euler's method over the interval [$x0$,
+--   $xN$]. The step size `h` will be calculated automatically. A list
+--   of y-values will be returned.
+--
+--   The explicit 'forall' in the type signature allows us to refer
+--   back to the type variables 'a' and 'b' in the 'where' clause.
+--
+--   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 10000
+--   >>> let yN = head $ reverse ys
+--   >>> abs ((exp 1) - yN) < 1/10^3
+--   True
+--
+eulers_method :: forall a b c. (RealFrac a, RealFrac b, Integral c)
+               => 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)
+               -> c -- ^ n, the number of intervals to use.
+               -> [b]
+eulers_method x0 y0 xN f n =
+  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
+    -- 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)
+      where
+        y1 = eulers_method1 x0 y0 f (x1 - x0)