Skip to content
This repository

HTTPS clone URL

Subversion checkout URL

You can clone with HTTPS or Subversion.

Download ZIP
Browse code

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

…entation of taylor models
  • Loading branch information...
commit 777a26ca32c61f014b7bf4d258692c9b29eed2c5 1 parent 601a431
wires authored
18 Data/Multinomial.hs
@@ -7,6 +7,7 @@ module Data.Multinomial
7 7 (Polynomial, xP, yP, zP, constP, o,
8 8 evalP, degree,
9 9 BoundPolynomial, boundPolynomial,
  10 + lagrange, integrate,
10 11 derive, dx, dy, dz,
11 12 Zero, One, Two, zero, one, two,
12 13 VectorQ, (<+>), (<*>)) where
@@ -188,3 +189,20 @@ instance (VectorQ a, Num a) => VectorQ (Polynomial a) where
188 189 (<+>) = (+)
189 190 q <*> (Poly p) = Poly $ map (q<*>) p
190 191 zeroVector = fromInteger 0
  192 +
  193 +-- lagrange polynomial helper
  194 +littlel :: Fractional a => [(a,a)] -> Int -> Polynomial a
  195 +littlel pts j = product . map base $ filter (/= j) [0..length pts - 1]
  196 + where base m = (constP . recip $ x j - x m) * (xP - constP (x m))
  197 + x i = fst $ pts !! i
  198 +
  199 +-- | given a list of points [(x0,y0)...] create lagrange polynomial p
  200 +-- such that p(x0) = y0, etc. This assumes all x values are distinct
  201 +lagrange :: Fractional a => [(a,a)] -> Polynomial a
  202 +lagrange pts = sum . map f $ [0..length pts - 1]
  203 + where f j = constP (snd $ pts !! j) * littlel pts j
  204 +
  205 +-- Using \int a*x^n = (a/n+1)*x^(n+1)
  206 +integrate :: Fractional a => Polynomial a -> Polynomial a
  207 +integrate (Poly p) = Poly $ 0:xs
  208 + where xs = map (\(n,a) -> a/fromInteger n) $ zip [1..] p
88 Data/Real/DReal.hs
@@ -425,90 +425,4 @@ danswer :: Int -> DReal -> String
425 425 danswer n x = show (round $ toQ $ drealScale (10^n) x (1 / 2)) ++ "x10^-" ++ (show n)
426 426
427 427 instance Show DReal where
428   - show = danswer 50
429   -
430   --- | half of (ceiling of 1 / 4-th root)
431   -g :: Rational -> Int
432   -g x | x == 0 = 0
433   - | True = (ceiling . (1/) . (2*) . toRational . sqrt . sqrt . fromRational) x
434   -
435   --- 'g' should be replaced by an exact function; the idea is to
436   --- approximate g to one digit and then check if this is indeed
437   --- the right value. If not it must be either g+1 or g-1
438   -
439   -
440   -{-| In order to approximate an integral up to precision eps
441   - using Simpsons rule one must modify the interval at which
442   - the function is evaluated. This function computes the
443   - interval to the forth power: h^4 = 180 * epsilon / (M * (b-a))
444   - This is done using rationals and gives an exact answer,
445   - assuming b > a, m > 0, eps > 0.
446   --}
447   -h4 :: Rational -> Rational -> Rational -> Gauge -> Rational
448   -h4 a b m eps | (b-a) > 0 = 180 * eps / (m * (b-a) ^^ 5)
449   - | otherwise = 0
450   -
451   -fromInt :: Num a => Int -> a
452   -fromInt = fromInteger . toEnum
453   -
454   --- | To apply an uniform continuous function to a rational and
455   --- get a result with precision eps, we need to approximate the
456   --- input at precision mu eps
457   -applyToRational :: (Dyadic :=> DReal) -> Rational -> DReal
458   -applyToRational f = bind f . fromRational
459   -
460   -dsimpson :: Rational -> Rational -> Rational -> (Dyadic :=> DReal) -> DReal
461   -dsimpson a b m f eps = h3 * (sum $ map f' ps)
462   - where
463   - h3 :: Dyadic
464   - h3 = fromRational (h/3) eps
465   - ps = [(1,a)] ++ pts ++ [(1,b)]
466   -
467   - n :: Int
468   - n = 1 + 2 * (g $ h4 a b m eps)
469   - -- ^ number of points (a b excluded)
470   - h :: Rational
471   - h = (b - a) / (toInteger (n + 1) % 1)
472   - -- ^ distance between two points
473   -
474   - pts :: [(Int,Rational)]
475   - pts = take n $ zip (cycle [4,2]) ahh
476   - -- ^ actual points 4*x1, 2*x2, 4*x3, ...
477   -
478   - ahh = scanl (+) (a+h) (repeat h)
479   -
480   - f' (c, x) = (fromInt c) * (applyToRational f x e)
481   - e = eps / (toInteger n % 1)
482   - -- ^ we are applying f n times, so need precision eps/n
483   -
484   --- example: dsimpson 1 2 1 (dLnUniformCts $ fromInteger 1 ) (1/10000000000000)
485   --- example: dsimpson 0 (31415729/5000000) 1 (dSinCts) (1/100)
486   --- example: dsimpson 0 1 1 idCts (1/1000) ≈ 1/2
487   -
488   --- | compose a function n times
489   -ncomp :: Int -> (a -> a) -> (a -> a)
490   -ncomp 1 = id
491   -ncomp n = \f -> f . ncomp (n - 1) f
492   -
493   --- contraction operator on the space of uniform cts real fns
494   -type Operator = (Dyadic :=> DReal) :=>> (Dyadic :=> DReal)
495   -
496   --- | computes for a given contraction and precision the required
497   --- picard steps, ie, nr of compositions in (F o F o ... o F)
498   -picard_steps :: Operator -> Gauge -> Int
499   -picard_steps f eps = ceiling $ logBase c e -- TODO use precise log!
500   - where c = fromRational . lipschitzConst $ f
501   - e = fromRational eps
502   -
503   -picard :: Operator -> Dyadic -> DReal
504   -picard op x eps = (forgetUniformCts opn) x eps
505   - where n = picard_steps op eps
506   - opn = ncomp n (forgetContraction op) idCts
507   -
508   --- | The initial value problem y'=y for y(0)=1 as intergral eqn
509   --- The exp function is the obvious solution to this IVP.
510   -expOp :: Operator
511   -expOp = mkContraction (1/2) fCts
512   - where fCts u = mkUniformCts (/2) (\t -> 1 + dsimpson 0 (toQ t) 3 u)
513   -
514   --- example: compute (exp .5): picard expOp (Dyadic 1 (- 1)) 1/10
  428 + show = danswer 50
83 Data/Real/Integration.hs
... ... @@ -0,0 +1,83 @@
  1 +{-# OPTIONS -XTypeOperators #-}
  2 +{-# OPTIONS -XTypeSynonymInstances #-}
  3 +
  4 +module Data.Real.Integration
  5 +where
  6 +
  7 +import Data.Real.Dyadic
  8 +import Data.Real.DReal
  9 +import Data.Real.Complete
  10 +import Data.Ratio
  11 +import Data.Multinomial
  12 +
  13 +-- | half of (ceiling of 1 / 4-th root)
  14 +g :: Rational -> Int
  15 +g x | x == 0 = 0
  16 + | True = (ceiling . (1/) . (2*) . toRational . sqrt . sqrt . fromRational) x
  17 +-- 'g' should be replaced by an exact function; the idea is to
  18 +-- approximate g to one digit and then check if this is indeed
  19 +-- the right value. If not it must be either g+1 or g-1
  20 +-- Alternatively compute as real real number and just add 2*eps
  21 +-- and take ceiling of that. This gives either g or g+1.
  22 +
  23 +{-| In order to approximate an integral up to precision eps
  24 + using Simpsons rule one must modify the interval at which
  25 + the function is evaluated. This function computes the
  26 + interval to the forth power: h^4 = 180 * epsilon / (M * (b-a))
  27 + This is done using rationals and gives an exact answer,
  28 + assuming b > a, m > 0, eps > 0.
  29 +-}
  30 +h4 :: Rational -> Rational -> Rational -> Gauge -> Rational
  31 +h4 a b m eps | (b-a) > 0 = 180 * eps / (m * (b-a) ^^ 5)
  32 + | otherwise = 0
  33 +
  34 +fromInt :: Num a => Int -> a
  35 +fromInt = fromInteger . toEnum
  36 +
  37 +-- | To apply an uniform continuous function to a rational and
  38 +-- get a result with precision eps, we need to approximate the
  39 +-- input at precision mu eps
  40 +applyToRational :: (Dyadic :=> DReal) -> Rational -> DReal
  41 +applyToRational f = bind f . fromRational
  42 +
  43 +simpson_steps :: Rational -> Rational -> Rational -> Gauge -> Int
  44 +simpson_steps a b m eps = 1 + 2 * (g $ h4 a b m eps)
  45 +
  46 +dsimpson :: Rational -> Rational -> Rational -> (Dyadic :=> DReal) -> DReal
  47 +dsimpson a b m f eps = h3 * (sum $ map f' ps)
  48 + where
  49 + h3 :: Dyadic
  50 + h3 = fromRational (h/3) eps
  51 + ps = [(1,a)] ++ pts ++ [(1,b)]
  52 +
  53 + n :: Int
  54 + n = 1 + 2 * (g $ h4 a b m eps)
  55 + -- ^ number of points (a b excluded)
  56 + h :: Rational
  57 + h = (b - a) / (toInteger (n + 1) % 1)
  58 + -- ^ distance between two points
  59 +
  60 + pts :: [(Int,Rational)]
  61 + pts = take n $ zip (cycle [4,2]) ahh
  62 + -- ^ actual points 4*x1, 2*x2, 4*x3, ...
  63 +
  64 + ahh = scanl (+) (a+h) (repeat h)
  65 +
  66 + f' (c, x) = (fromInt c) * (applyToRational f x e)
  67 + e = eps / (toInteger n % 1)
  68 + -- ^ we are applying f n times, so need precision eps/n
  69 +
  70 +-- example: dsimpson 1 2 1 (dLnUniformCts $ fromInteger 1 ) (1/10000000000000)
  71 +-- example: dsimpson 0 (31415729/5000000) 1 (dSinCts) (1/100)
  72 +-- example: dsimpson 0 1 1 idCts (1/1000) ≈ 1/2
  73 +
  74 +subdivisions :: (Fractional a) => a -> a -> Int -> [a]
  75 +subdivisions lb ub n = (take n pts) ++ [ub]
  76 + where pts = scanl (+) lb (repeat h)
  77 + h = (ub - lb) / (fromInteger $ toEnum n)
  78 +
  79 +polyapprox a b m f eps = lagrange (zip pts $ map f' pts)
  80 + where pts = subdivisions a b (simpson_steps a b m eps)
  81 + f' = toQ . ($ eps) . (bind f) . fromRational
  82 +
  83 +symbolicSimpson a b m f eps = integrate $ polyapprox a b m f eps
47 Data/Real/Picard.hs
... ... @@ -0,0 +1,47 @@
  1 +{-# OPTIONS -XTypeOperators #-}
  2 +{-# OPTIONS -XTypeSynonymInstances #-}
  3 +
  4 +module Data.Real.Picard
  5 +where
  6 +
  7 +import Data.Real.Dyadic
  8 +import Data.Real.DReal
  9 +import Data.Real.Integration
  10 +import Data.Real.Complete
  11 +import Data.Ratio
  12 +
  13 +-- | compose a function n times
  14 +ncomp :: Int -> (a -> a) -> (a -> a)
  15 +ncomp 1 = id
  16 +ncomp n = \f -> f . ncomp (n - 1) f
  17 +
  18 +-- contraction operator on the space of uniform cts real fns
  19 +type Operator a b = (a :=> b) :=>> (a :=> b)
  20 +
  21 +-- | computes for a given contraction and precision the required
  22 +-- picard steps, ie, nr of compositions in (F o F o ... o F)
  23 +picard_steps :: Operator a b -> Gauge -> Int
  24 +picard_steps f eps = ceiling $ logBase c e -- TODO use precise log!
  25 + where c = fromRational . lipschitzConst $ f
  26 + e = fromRational eps
  27 +
  28 +picard :: Operator a (Complete a) -> a -> Complete a
  29 +picard op x eps = (forgetUniformCts opn) x eps
  30 + where n = picard_steps op eps
  31 + opn = ncomp n (forgetContraction op) idCts
  32 +
  33 +-- | The initial value problem y'=y for y(0)=1 as intergral eqn
  34 +-- The exp function is the obvious solution to this IVP.
  35 +expOp :: Operator Dyadic DReal
  36 +expOp = mkContraction (1/2) fCts
  37 + where fCts u = mkUniformCts (/2) (\t -> 1 + dsimpson 0 (toQ t) 3 u)
  38 +
  39 +-- example: compute (exp .5): picard expOp (Dyadic 1 (- 1)) 1/10
  40 +
  41 +{- bla
  42 +symbolicPicard :: Fractional a => Operator a (Complete a) -> a -> Polynomial a
  43 +symbolicPicard op x eps = foldl (flip (.)) toP (replicate n integrate)
  44 + where toP :: a -> Polynomial a
  45 + toP x = symbolicSimpson 0 x 1 (asUniformCts op) eps
  46 + n = max 0 $ picard_steps op eps - 1
  47 +-}

0 comments on commit 777a26c

Please sign in to comment.
Something went wrong with that request. Please try again.