Convert src/ODE/IVP.hs to numeric-prelude.
Add RealRing/RealField instanced to BigFloat.
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
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)
{-# LANGUAGE ScopedTypeVariables #-}
{-# LANGUAGE ScopedTypeVariables #-}
+{-# LANGUAGE RebindableSyntax #-}
-- | Numerical solution of the initial value problem,
--
-- | Numerical solution of the initial value problem,
--
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$].
-- >>> 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)
eulers_method1 x0 y0 f h =
y0 + h'*y'
where
eulers_method1 x0 y0 f h =
y0 + h'*y'
where
+ h' = fromRational'$ toRational 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
-- >>> 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
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@
-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@
-- 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@.
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
-- >>> 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@
-- >>> 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@
-- 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@.
-- 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@.
--
-- 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@.
-- >>> 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
-- >>> 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
-- >>> 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
-- >>> 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