X-Git-Url: http://gitweb.michael.orlitzky.com/?a=blobdiff_plain;f=src%2FODE%2FIVP.hs;h=5bffe08cfff387210c3bdd7cc3f5483dba81c13e;hb=e1bf90ef57aaa37620e001d04cabec0ea2e8ddbf;hp=381ff8a9c98099519dae2e08454fa8f05260127c;hpb=ddad82e5f42a7d38d542fbb77e981dc7309089f4;p=numerical-analysis.git diff --git a/src/ODE/IVP.hs b/src/ODE/IVP.hs index 381ff8a..5bffe08 100644 --- a/src/ODE/IVP.hs +++ b/src/ODE/IVP.hs @@ -1,4 +1,5 @@ {-# LANGUAGE ScopedTypeVariables #-} +{-# LANGUAGE RebindableSyntax #-} -- | Numerical solution of the initial value problem, -- @@ -8,11 +9,18 @@ -- for x in [x0, xN]. -- -module ODE.IVP +module ODE.IVP ( + eulers_method, + eulers_method1, + eulers_methodH ) where -import Misc (partition) - +import Misc ( partition ) +import NumericPrelude hiding ( abs ) +import qualified Algebra.Field as Field ( C ) +import qualified Algebra.ToInteger as ToInteger ( C ) +import qualified Algebra.ToRational as ToRational ( C ) +import qualified Algebra.RealField as RealField ( C ) -- | A single iteration of Euler's method over the interval -- [$x0$, $x0$+$h$]. @@ -26,7 +34,7 @@ import Misc (partition) -- >>> eulers_method1 x0 y0 f h -- 2.0 -- -eulers_method1 :: (RealFrac a, RealFrac b) +eulers_method1 :: forall a b. (Field.C a, ToRational.C a, Field.C b) => a -- ^ x0, the initial point -> b -- ^ y0, the initial value at x0 -> (a -> b -> b) -- ^ The function f(x,y) @@ -35,7 +43,7 @@ eulers_method1 :: (RealFrac a, RealFrac b) eulers_method1 x0 y0 f h = y0 + h'*y' where - h' = fromRational $ toRational h + h' = fromRational'$ toRational h :: b y' = (f x0 y0) @@ -48,6 +56,7 @@ eulers_method1 x0 y0 f h = -- -- Examples: -- +-- >>> import Algebra.Absolute (abs) -- >>> let x0 = 0.0 -- >>> let xN = 1.0 -- >>> let y0 = 1.0 @@ -56,8 +65,14 @@ 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) +eulers_method :: forall a b c. (Field.C a, + ToRational.C a, + Field.C b, + ToInteger.C c, + Enum c) => a -- ^ x0, the initial point -> b -- ^ y0, the initial value at x0 -> a -- ^ xN, the terminal point @@ -65,16 +80,50 @@ 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 -- 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 :: (RealField.C a, ToRational.C a, Field.C 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