Permalink
Browse files

moved interation & picard iteration to seperate files, initial implem…

…entation of taylor models
  • Loading branch information...
1 parent 601a431 commit 777a26ca32c61f014b7bf4d258692c9b29eed2c5 @wires committed Feb 28, 2012
Showing with 149 additions and 87 deletions.
  1. +18 −0 Data/Multinomial.hs
  2. +1 −87 Data/Real/DReal.hs
  3. +83 −0 Data/Real/Integration.hs
  4. +47 −0 Data/Real/Picard.hs
View
18 Data/Multinomial.hs
@@ -7,6 +7,7 @@ module Data.Multinomial
(Polynomial, xP, yP, zP, constP, o,
evalP, degree,
BoundPolynomial, boundPolynomial,
+ lagrange, integrate,
derive, dx, dy, dz,
Zero, One, Two, zero, one, two,
VectorQ, (<+>), (<*>)) where
@@ -188,3 +189,20 @@ instance (VectorQ a, Num a) => VectorQ (Polynomial a) where
(<+>) = (+)
q <*> (Poly p) = Poly $ map (q<*>) p
zeroVector = fromInteger 0
+
+-- lagrange polynomial helper
+littlel :: Fractional a => [(a,a)] -> Int -> Polynomial a
+littlel pts j = product . map base $ filter (/= j) [0..length pts - 1]
+ where base m = (constP . recip $ x j - x m) * (xP - constP (x m))
+ x i = fst $ pts !! i
+
+-- | given a list of points [(x0,y0)...] create lagrange polynomial p
+-- such that p(x0) = y0, etc. This assumes all x values are distinct
+lagrange :: Fractional a => [(a,a)] -> Polynomial a
+lagrange pts = sum . map f $ [0..length pts - 1]
+ where f j = constP (snd $ pts !! j) * littlel pts j
+
+-- Using \int a*x^n = (a/n+1)*x^(n+1)
+integrate :: Fractional a => Polynomial a -> Polynomial a
+integrate (Poly p) = Poly $ 0:xs
+ where xs = map (\(n,a) -> a/fromInteger n) $ zip [1..] p
View
88 Data/Real/DReal.hs
@@ -425,90 +425,4 @@ danswer :: Int -> DReal -> String
danswer n x = show (round $ toQ $ drealScale (10^n) x (1 / 2)) ++ "x10^-" ++ (show n)
instance Show DReal where
- show = danswer 50
-
--- | half of (ceiling of 1 / 4-th root)
-g :: Rational -> Int
-g x | x == 0 = 0
- | True = (ceiling . (1/) . (2*) . toRational . sqrt . sqrt . fromRational) x
-
--- 'g' should be replaced by an exact function; the idea is to
--- approximate g to one digit and then check if this is indeed
--- the right value. If not it must be either g+1 or g-1
-
-
-{-| In order to approximate an integral up to precision eps
- using Simpsons rule one must modify the interval at which
- the function is evaluated. This function computes the
- interval to the forth power: h^4 = 180 * epsilon / (M * (b-a))
- This is done using rationals and gives an exact answer,
- assuming b > a, m > 0, eps > 0.
--}
-h4 :: Rational -> Rational -> Rational -> Gauge -> Rational
-h4 a b m eps | (b-a) > 0 = 180 * eps / (m * (b-a) ^^ 5)
- | otherwise = 0
-
-fromInt :: Num a => Int -> a
-fromInt = fromInteger . toEnum
-
--- | To apply an uniform continuous function to a rational and
--- get a result with precision eps, we need to approximate the
--- input at precision mu eps
-applyToRational :: (Dyadic :=> DReal) -> Rational -> DReal
-applyToRational f = bind f . fromRational
-
-dsimpson :: Rational -> Rational -> Rational -> (Dyadic :=> DReal) -> DReal
-dsimpson a b m f eps = h3 * (sum $ map f' ps)
- where
- h3 :: Dyadic
- h3 = fromRational (h/3) eps
- ps = [(1,a)] ++ pts ++ [(1,b)]
-
- n :: Int
- n = 1 + 2 * (g $ h4 a b m eps)
- -- ^ number of points (a b excluded)
- h :: Rational
- h = (b - a) / (toInteger (n + 1) % 1)
- -- ^ distance between two points
-
- pts :: [(Int,Rational)]
- pts = take n $ zip (cycle [4,2]) ahh
- -- ^ actual points 4*x1, 2*x2, 4*x3, ...
-
- ahh = scanl (+) (a+h) (repeat h)
-
- f' (c, x) = (fromInt c) * (applyToRational f x e)
- e = eps / (toInteger n % 1)
- -- ^ we are applying f n times, so need precision eps/n
-
--- example: dsimpson 1 2 1 (dLnUniformCts $ fromInteger 1 ) (1/10000000000000)
--- example: dsimpson 0 (31415729/5000000) 1 (dSinCts) (1/100)
--- example: dsimpson 0 1 1 idCts (1/1000) ≈ 1/2
-
--- | compose a function n times
-ncomp :: Int -> (a -> a) -> (a -> a)
-ncomp 1 = id
-ncomp n = \f -> f . ncomp (n - 1) f
-
--- contraction operator on the space of uniform cts real fns
-type Operator = (Dyadic :=> DReal) :=>> (Dyadic :=> DReal)
-
--- | computes for a given contraction and precision the required
--- picard steps, ie, nr of compositions in (F o F o ... o F)
-picard_steps :: Operator -> Gauge -> Int
-picard_steps f eps = ceiling $ logBase c e -- TODO use precise log!
- where c = fromRational . lipschitzConst $ f
- e = fromRational eps
-
-picard :: Operator -> Dyadic -> DReal
-picard op x eps = (forgetUniformCts opn) x eps
- where n = picard_steps op eps
- opn = ncomp n (forgetContraction op) idCts
-
--- | The initial value problem y'=y for y(0)=1 as intergral eqn
--- The exp function is the obvious solution to this IVP.
-expOp :: Operator
-expOp = mkContraction (1/2) fCts
- where fCts u = mkUniformCts (/2) (\t -> 1 + dsimpson 0 (toQ t) 3 u)
-
--- example: compute (exp .5): picard expOp (Dyadic 1 (- 1)) 1/10
+ show = danswer 50
View
83 Data/Real/Integration.hs
@@ -0,0 +1,83 @@
+{-# OPTIONS -XTypeOperators #-}
+{-# OPTIONS -XTypeSynonymInstances #-}
+
+module Data.Real.Integration
+where
+
+import Data.Real.Dyadic
+import Data.Real.DReal
+import Data.Real.Complete
+import Data.Ratio
+import Data.Multinomial
+
+-- | half of (ceiling of 1 / 4-th root)
+g :: Rational -> Int
+g x | x == 0 = 0
+ | True = (ceiling . (1/) . (2*) . toRational . sqrt . sqrt . fromRational) x
+-- 'g' should be replaced by an exact function; the idea is to
+-- approximate g to one digit and then check if this is indeed
+-- the right value. If not it must be either g+1 or g-1
+-- Alternatively compute as real real number and just add 2*eps
+-- and take ceiling of that. This gives either g or g+1.
+
+{-| In order to approximate an integral up to precision eps
+ using Simpsons rule one must modify the interval at which
+ the function is evaluated. This function computes the
+ interval to the forth power: h^4 = 180 * epsilon / (M * (b-a))
+ This is done using rationals and gives an exact answer,
+ assuming b > a, m > 0, eps > 0.
+-}
+h4 :: Rational -> Rational -> Rational -> Gauge -> Rational
+h4 a b m eps | (b-a) > 0 = 180 * eps / (m * (b-a) ^^ 5)
+ | otherwise = 0
+
+fromInt :: Num a => Int -> a
+fromInt = fromInteger . toEnum
+
+-- | To apply an uniform continuous function to a rational and
+-- get a result with precision eps, we need to approximate the
+-- input at precision mu eps
+applyToRational :: (Dyadic :=> DReal) -> Rational -> DReal
+applyToRational f = bind f . fromRational
+
+simpson_steps :: Rational -> Rational -> Rational -> Gauge -> Int
+simpson_steps a b m eps = 1 + 2 * (g $ h4 a b m eps)
+
+dsimpson :: Rational -> Rational -> Rational -> (Dyadic :=> DReal) -> DReal
+dsimpson a b m f eps = h3 * (sum $ map f' ps)
+ where
+ h3 :: Dyadic
+ h3 = fromRational (h/3) eps
+ ps = [(1,a)] ++ pts ++ [(1,b)]
+
+ n :: Int
+ n = 1 + 2 * (g $ h4 a b m eps)
+ -- ^ number of points (a b excluded)
+ h :: Rational
+ h = (b - a) / (toInteger (n + 1) % 1)
+ -- ^ distance between two points
+
+ pts :: [(Int,Rational)]
+ pts = take n $ zip (cycle [4,2]) ahh
+ -- ^ actual points 4*x1, 2*x2, 4*x3, ...
+
+ ahh = scanl (+) (a+h) (repeat h)
+
+ f' (c, x) = (fromInt c) * (applyToRational f x e)
+ e = eps / (toInteger n % 1)
+ -- ^ we are applying f n times, so need precision eps/n
+
+-- example: dsimpson 1 2 1 (dLnUniformCts $ fromInteger 1 ) (1/10000000000000)
+-- example: dsimpson 0 (31415729/5000000) 1 (dSinCts) (1/100)
+-- example: dsimpson 0 1 1 idCts (1/1000) ≈ 1/2
+
+subdivisions :: (Fractional a) => a -> a -> Int -> [a]
+subdivisions lb ub n = (take n pts) ++ [ub]
+ where pts = scanl (+) lb (repeat h)
+ h = (ub - lb) / (fromInteger $ toEnum n)
+
+polyapprox a b m f eps = lagrange (zip pts $ map f' pts)
+ where pts = subdivisions a b (simpson_steps a b m eps)
+ f' = toQ . ($ eps) . (bind f) . fromRational
+
+symbolicSimpson a b m f eps = integrate $ polyapprox a b m f eps
View
47 Data/Real/Picard.hs
@@ -0,0 +1,47 @@
+{-# OPTIONS -XTypeOperators #-}
+{-# OPTIONS -XTypeSynonymInstances #-}
+
+module Data.Real.Picard
+where
+
+import Data.Real.Dyadic
+import Data.Real.DReal
+import Data.Real.Integration
+import Data.Real.Complete
+import Data.Ratio
+
+-- | compose a function n times
+ncomp :: Int -> (a -> a) -> (a -> a)
+ncomp 1 = id
+ncomp n = \f -> f . ncomp (n - 1) f
+
+-- contraction operator on the space of uniform cts real fns
+type Operator a b = (a :=> b) :=>> (a :=> b)
+
+-- | computes for a given contraction and precision the required
+-- picard steps, ie, nr of compositions in (F o F o ... o F)
+picard_steps :: Operator a b -> Gauge -> Int
+picard_steps f eps = ceiling $ logBase c e -- TODO use precise log!
+ where c = fromRational . lipschitzConst $ f
+ e = fromRational eps
+
+picard :: Operator a (Complete a) -> a -> Complete a
+picard op x eps = (forgetUniformCts opn) x eps
+ where n = picard_steps op eps
+ opn = ncomp n (forgetContraction op) idCts
+
+-- | The initial value problem y'=y for y(0)=1 as intergral eqn
+-- The exp function is the obvious solution to this IVP.
+expOp :: Operator Dyadic DReal
+expOp = mkContraction (1/2) fCts
+ where fCts u = mkUniformCts (/2) (\t -> 1 + dsimpson 0 (toQ t) 3 u)
+
+-- example: compute (exp .5): picard expOp (Dyadic 1 (- 1)) 1/10
+
+{- bla
+symbolicPicard :: Fractional a => Operator a (Complete a) -> a -> Polynomial a
+symbolicPicard op x eps = foldl (flip (.)) toP (replicate n integrate)
+ where toP :: a -> Polynomial a
+ toP x = symbolicSimpson 0 x 1 (asUniformCts op) eps
+ n = max 0 $ picard_steps op eps - 1
+-}

0 comments on commit 777a26c

Please sign in to comment.