]> gitweb.michael.orlitzky.com - numerical-analysis.git/commitdiff
Use RealField/RealRing where possible instead of their constituents.
authorMichael Orlitzky <michael@orlitzky.com>
Wed, 20 Feb 2013 23:56:13 +0000 (18:56 -0500)
committerMichael Orlitzky <michael@orlitzky.com>
Wed, 20 Feb 2013 23:56:13 +0000 (18:56 -0500)
Convert src/ODE/IVP.hs to numeric-prelude.
Add RealRing/RealField instanced to BigFloat.

src/BigFloat.hs
src/ODE/IVP.hs
src/Roots/Fast.hs
src/Roots/Simple.hs

index e23ebe6833bd674be0d5ea1f65dbfe540ec72baf..0bebab9659a2fd8ee22358a39c676661e2fba827 100644 (file)
@@ -10,6 +10,8 @@ import NumericPrelude hiding (abs)
 import qualified Algebra.Absolute as Absolute
 import qualified Algebra.Additive as Additive
 import qualified Algebra.Field as Field
 import qualified Algebra.Absolute as Absolute
 import qualified Algebra.Additive as Additive
 import qualified Algebra.Field as Field
+import qualified Algebra.RealField as RealField
+import qualified Algebra.RealRing as RealRing
 import qualified Algebra.Ring as Ring
 import qualified Algebra.ToRational as ToRational
 import qualified Algebra.ZeroTestable as ZeroTestable
 import qualified Algebra.Ring as Ring
 import qualified Algebra.ToRational as ToRational
 import qualified Algebra.ZeroTestable as ZeroTestable
@@ -38,3 +40,8 @@ instance Epsilon e => ZeroTestable.C (BigFloat e) where
 
 instance Epsilon e => ToRational.C (BigFloat e) where
   toRational = fromRational . P.toRational
 
 instance Epsilon e => ToRational.C (BigFloat e) where
   toRational = fromRational . P.toRational
+
+instance Epsilon e => RealRing.C (BigFloat e) where
+  floor = fromInteger . P.floor
+
+instance Epsilon e => RealField.C (BigFloat e)
index 6bde763a8e264549b8dae24ecbe945b57c0b4e99..b4fd3dc5a37d32b29ef020a06bba47e0ef4de3ba 100644 (file)
@@ -1,4 +1,5 @@
 {-# LANGUAGE ScopedTypeVariables #-}
 {-# LANGUAGE ScopedTypeVariables #-}
+{-# LANGUAGE RebindableSyntax #-}
 
 -- | Numerical solution of the initial value problem,
 --
 
 -- | Numerical solution of the initial value problem,
 --
@@ -12,7 +13,12 @@ module ODE.IVP
 where
 
 import Misc (partition)
 where
 
 import Misc (partition)
-
+import NumericPrelude hiding (abs)
+import Algebra.Absolute (abs)
+import qualified Algebra.Field as Field
+import qualified Algebra.ToInteger as ToInteger
+import qualified Algebra.ToRational as ToRational
+import qualified Algebra.RealField as RealField
 
 -- | A single iteration of Euler's method over the interval
 --   [$x0$, $x0$+$h$].
 
 -- | A single iteration of Euler's method over the interval
 --   [$x0$, $x0$+$h$].
@@ -26,7 +32,7 @@ import Misc (partition)
 --   >>> eulers_method1 x0 y0 f h
 --   2.0
 --
 --   >>> eulers_method1 x0 y0 f h
 --   2.0
 --
-eulers_method1 :: (RealFrac a, RealFrac b)
+eulers_method1 :: (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)
                => a -- ^ x0, the initial point
                -> b -- ^ y0, the initial value at x0
                -> (a -> b -> b) -- ^ The function f(x,y)
@@ -35,7 +41,7 @@ eulers_method1 :: (RealFrac a, RealFrac b)
 eulers_method1 x0 y0 f h =
   y0 +  h'*y'
   where
 eulers_method1 x0 y0 f h =
   y0 +  h'*y'
   where
-    h' = realToFrac h
+    h' = fromRational'$ toRational h
     y' = (f x0 y0)
 
 
     y' = (f x0 y0)
 
 
@@ -59,7 +65,11 @@ eulers_method1 x0 y0 f h =
 --   >>> head ys == y0
 --   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
                => a -- ^ x0, the initial point
                -> b -- ^ y0, the initial value at x0
                -> a -- ^ xN, the terminal point
@@ -102,7 +112,7 @@ eulers_method x0 y0 xN f n =
 --   >>> ys == ys'
 --   True
 --
 --   >>> ys == ys'
 --   True
 --
-eulers_methodH :: (RealFrac a, RealFrac b)
+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 -- ^ x0, the initial point
                -> b -- ^ y0, the initial value at x0
                -> a -- ^ xN, the terminal point
index 47fa512a75476310403d1faed03b9c54bdabe171..78c299ad6a2cc39629b409bff5eed86787a3ca20 100644 (file)
@@ -13,15 +13,14 @@ import Data.List (find)
 import Normed
 
 import NumericPrelude hiding (abs)
 import Normed
 
 import NumericPrelude hiding (abs)
-import Algebra.Absolute
-import Algebra.Field
-import Algebra.Ring
-
-has_root :: (Algebra.Field.C a,
-             Ord a,
-             Algebra.Ring.C b,
-             Ord b,
-             Algebra.Absolute.C b)
+import qualified Algebra.Absolute as Absolute
+import qualified Algebra.Field as Field
+import qualified Algebra.RealRing as RealRing
+import qualified Algebra.RealField as RealField
+
+has_root :: (RealField.C a,
+             RealRing.C b,
+             Absolute.C b)
          => (a -> b) -- ^ The function @f@
          -> a       -- ^ The \"left\" endpoint, @a@
          -> a       -- ^ The \"right\" endpoint, @b@
          => (a -> b) -- ^ The function @f@
          -> a       -- ^ The \"left\" endpoint, @a@
          -> a       -- ^ The \"right\" endpoint, @b@
@@ -61,11 +60,9 @@ has_root f a b epsilon f_of_a f_of_b =
     c = (a + b)/2
 
 
     c = (a + b)/2
 
 
-bisect :: (Algebra.Field.C a,
-           Ord a,
-           Algebra.Ring.C b,
-           Ord b,
-           Algebra.Absolute.C b)
+bisect :: (RealField.C a,
+           RealRing.C b,
+           Absolute.C b)
        => (a -> b) -- ^ The function @f@ whose root we seek
        -> a       -- ^ The \"left\" endpoint of the interval, @a@
        -> a       -- ^ The \"right\" endpoint of the interval, @b@
        => (a -> b) -- ^ The function @f@ whose root we seek
        -> a       -- ^ The \"left\" endpoint of the interval, @a@
        -> a       -- ^ The \"right\" endpoint of the interval, @b@
@@ -119,8 +116,8 @@ fixed_point_iterations f x0 =
 --   We also return the number of iterations required.
 --
 fixed_point_with_iterations :: (Normed a,
 --   We also return the number of iterations required.
 --
 fixed_point_with_iterations :: (Normed a,
-                                Algebra.Field.C b,
-                                Algebra.Absolute.C b,
+                                Field.C b,
+                                Absolute.C b,
                                 Ord b)
                             => (a -> a)  -- ^ The function @f@ to iterate.
                             -> b        -- ^ The tolerance, @epsilon@.
                                 Ord b)
                             => (a -> a)  -- ^ The function @f@ to iterate.
                             -> b        -- ^ The tolerance, @epsilon@.
index 44d3d62112d2c4f339fa16b65c143e65ed5bb83a..0a1debff0a0b1db653efabdf2be00ce9a4b91ed4 100644 (file)
@@ -18,9 +18,11 @@ import Normed
 import qualified Roots.Fast as F
 
 import NumericPrelude hiding (abs)
 import qualified Roots.Fast as F
 
 import NumericPrelude hiding (abs)
-import Algebra.Absolute
-import Algebra.Field
-import Algebra.Ring
+import qualified Algebra.Absolute as Absolute
+import Algebra.Absolute (abs)
+import qualified Algebra.Field as Field
+import qualified Algebra.RealField as RealField
+import qualified Algebra.RealRing as RealRing
 
 -- | Does the (continuous) function @f@ have a root on the interval
 --   [a,b]? If f(a) <] 0 and f(b) ]> 0, we know that there's a root in
 
 -- | Does the (continuous) function @f@ have a root on the interval
 --   [a,b]? If f(a) <] 0 and f(b) ]> 0, we know that there's a root in
@@ -41,11 +43,7 @@ import Algebra.Ring
 --   >>> has_root cos (-2) 2 (Just 0.001)
 --   True
 --
 --   >>> has_root cos (-2) 2 (Just 0.001)
 --   True
 --
-has_root :: (Algebra.Field.C a,
-             Ord a,
-             Algebra.Ring.C b,
-             Algebra.Absolute.C b,
-             Ord b)
+has_root :: (RealField.C a, RealRing.C b)
          => (a -> b) -- ^ The function @f@
          -> a       -- ^ The \"left\" endpoint, @a@
          -> a       -- ^ The \"right\" endpoint, @b@
          => (a -> b) -- ^ The function @f@
          -> a       -- ^ The \"left\" endpoint, @a@
          -> a       -- ^ The \"right\" endpoint, @b@
@@ -75,11 +73,7 @@ has_root f a b epsilon =
 --   >>> bisect sin (-1) 1 0.001
 --   Just 0.0
 --
 --   >>> bisect sin (-1) 1 0.001
 --   Just 0.0
 --
-bisect :: (Algebra.Field.C a,
-           Ord a,
-           Algebra.Ring.C b,
-           Algebra.Absolute.C b,
-           Ord b)
+bisect :: (RealField.C a, RealRing.C b)
        => (a -> b) -- ^ The function @f@ whose root we seek
        -> a       -- ^ The \"left\" endpoint of the interval, @a@
        -> a       -- ^ The \"right\" endpoint of the interval, @b@
        => (a -> b) -- ^ The function @f@ whose root we seek
        -> a       -- ^ The \"left\" endpoint of the interval, @a@
        -> a       -- ^ The \"right\" endpoint of the interval, @b@
@@ -93,10 +87,7 @@ bisect f a b epsilon =
 --   at x0. We delegate to the version that returns the number of
 --   iterations and simply discard the number of iterations.
 --
 --   at x0. We delegate to the version that returns the number of
 --   iterations and simply discard the number of iterations.
 --
-fixed_point :: (Normed a,
-                Algebra.Field.C b,
-                Algebra.Absolute.C b,
-                Ord b)
+fixed_point :: (Normed a, RealField.C b)
             => (a -> a) -- ^ The function @f@ to iterate.
             -> b       -- ^ The tolerance, @epsilon@.
             -> a       -- ^ The initial value @x0@.
             => (a -> a) -- ^ The function @f@ to iterate.
             -> b       -- ^ The tolerance, @epsilon@.
             -> a       -- ^ The initial value @x0@.
@@ -109,10 +100,7 @@ fixed_point f epsilon x0 =
 --   the function @f@ with the search starting at x0 and tolerance
 --   @epsilon@. We delegate to the version that returns the number of
 --   iterations and simply discard the fixed point.
 --   the function @f@ with the search starting at x0 and tolerance
 --   @epsilon@. We delegate to the version that returns the number of
 --   iterations and simply discard the fixed point.
-fixed_point_iteration_count :: (Normed a,
-                                Algebra.Field.C b,
-                                Algebra.Absolute.C b,
-                                Ord b)
+fixed_point_iteration_count :: (Normed a, RealField.C b)
                             => (a -> a) -- ^ The function @f@ to iterate.
                             -> b       -- ^ The tolerance, @epsilon@.
                             -> a       -- ^ The initial value @x0@.
                             => (a -> a) -- ^ The function @f@ to iterate.
                             -> b       -- ^ The tolerance, @epsilon@.
                             -> a       -- ^ The initial value @x0@.
@@ -130,10 +118,7 @@ fixed_point_iteration_count f epsilon x0 =
 --
 --   This is used to determine the rate of convergence.
 --
 --
 --   This is used to determine the rate of convergence.
 --
-fixed_point_error_ratios :: (Normed a,
-                             Algebra.Field.C b,
-                             Algebra.Absolute.C b,
-                             Ord b)
+fixed_point_error_ratios :: (Normed a, RealField.C b)
                    => (a -> a) -- ^ The function @f@ to iterate.
                    -> a       -- ^ The initial value @x0@.
                    -> a       -- ^ The true solution, @x_star@.
                    => (a -> a) -- ^ The function @f@ to iterate.
                    -> a       -- ^ The initial value @x0@.
                    -> a       -- ^ The true solution, @x_star@.
@@ -160,7 +145,7 @@ fixed_point_error_ratios f x0 x_star p =
 --   >>> tail $ take 4 $ newton_iterations f f' 2
 --   [1.6806282722513088,1.4307389882390624,1.2549709561094362]
 --
 --   >>> tail $ take 4 $ newton_iterations f f' 2
 --   [1.6806282722513088,1.4307389882390624,1.2549709561094362]
 --
-newton_iterations :: (Algebra.Field.C a)
+newton_iterations :: (Field.C a)
                     => (a -> a) -- ^ The function @f@ whose root we seek
                     -> (a -> a) -- ^ The derivative of @f@
                     -> a       -- ^ Initial guess, x-naught
                     => (a -> a) -- ^ The function @f@ whose root we seek
                     -> (a -> a) -- ^ The derivative of @f@
                     -> a       -- ^ Initial guess, x-naught
@@ -195,7 +180,7 @@ newton_iterations f f' x0 =
 --   >>> abs (f root) < eps
 --   True
 --
 --   >>> abs (f root) < eps
 --   True
 --
-newtons_method :: (Algebra.Field.C a, Algebra.Absolute.C a, Ord a)
+newtons_method :: (RealField.C a)
                  => (a -> a) -- ^ The function @f@ whose root we seek
                  -> (a -> a) -- ^ The derivative of @f@
                  -> a       -- ^ The tolerance epsilon
                  => (a -> a) -- ^ The function @f@ whose root we seek
                  -> (a -> a) -- ^ The derivative of @f@
                  -> a       -- ^ The tolerance epsilon
@@ -245,7 +230,7 @@ iterate2 f x0 x1 =
 --   >>> take 4 $ secant_iterations f 2 1
 --   [2.0,1.0,1.0161290322580645,1.190577768676638]
 --
 --   >>> take 4 $ secant_iterations f 2 1
 --   [2.0,1.0,1.0161290322580645,1.190577768676638]
 --
-secant_iterations :: (Algebra.Field.C a)
+secant_iterations :: (Field.C a)
                     => (a -> a) -- ^ The function @f@ whose root we seek
                     -> a       -- ^ Initial guess, x-naught
                     -> a       -- ^ Second initial guess, x-one
                     => (a -> a) -- ^ The function @f@ whose root we seek
                     -> a       -- ^ Initial guess, x-naught
                     -> a       -- ^ Second initial guess, x-one
@@ -273,7 +258,7 @@ secant_iterations f x0 x1 =
 --   >>> abs (f root) < (1/10^9)
 --   True
 --
 --   >>> abs (f root) < (1/10^9)
 --   True
 --
-secant_method :: (Algebra.Field.C a, Algebra.Absolute.C a, Ord a)
+secant_method :: (RealField.C a)
                  => (a -> a) -- ^ The function @f@ whose root we seek
                  -> a       -- ^ The tolerance epsilon
                  -> a       -- ^ Initial guess, x-naught
                  => (a -> a) -- ^ The function @f@ whose root we seek
                  -> a       -- ^ The tolerance epsilon
                  -> a       -- ^ Initial guess, x-naught