Ticket #10948: Interpolation.hs

File Interpolation.hs, 3.3 KB (added by devnull, 4 years ago)
Line 
1{-# LANGUAGE ParallelListComp #-}
2{-# OPTIONS_GHC -O3 #-}
3
4import           Data.List               (zip5)
5import           Graphics.Gnuplot.Simple
6import           Numeric.FastMath
7import           Text.Printf             (printf)
8
9f_sum :: Num a => [(t -> a)] -> t -> a
10f_sum xs = \x -> sum [f x | f <- xs]
11
12f_product :: Num a => [(t -> a)] -> t -> a
13f_product xs = \x -> product [f x | f <- xs]
14
15delete :: Int -> [a] -> [a]
16delete k xs = let (ys, zs) = splitAt k xs in ys ++ tail zs
17
18fact :: (Enum a, Num a) => a -> a
19fact n = product [1..n]
20
21lagrange :: Fractional a => [a] -> [a] -> a -> a
22lagrange xs ys
23    | n /= length ys = error "Incorrect input size"
24    | otherwise   = f_sum [\x -> f x * (y / a) | f <- l1 | a <- l2 | y <- ys]
25    where n  = length xs
26          l1 = map f_product [[\x -> x - a | a <- delete k xs] | k <- [0..n-1]]
27          l2 = map product [let xk = xs !! k in [xk - a | a <- delete k xs] | k <- [0..n-1]]
28
29step :: Fractional a => a -> a -> Int -> a
30step a b n = abs (b - a) / (fromIntegral $ n - 1)
31
32memo2 :: (Enum a, Num a) => (a -> a -> b) -> [[b]]
33memo2 f = map (\x -> map (f x) [0..]) [0..]
34
35newtonK :: Fractional a => Int -> a -> a -> [a] -> a -> a
36newtonK k x0 h deltas x = q k x / (fromIntegral . fact $ k) * (deltas !! k)
37    where q k           = f_product [\x -> (x - x0) / h - (fromIntegral a) | a <- [0..k-1]]
38
39newton :: Fractional a => a -> a -> [a] -> a -> a
40newton a b ys        = f_sum [newtonK k a h deltas | k <- [0..n-1]]
41    where h          = step a b n
42          n          = length ys
43          deltas     = map (\(a, b) -> delta' a b) [(i, 0) | i <- [0..n-1]]
44          delta 0 n  = ys !! n
45          delta 1 n  = ys !! (n + 1) - ys !! n
46          delta m n  = delta' (m - 1) (n + 1) - delta' (m - 1) n
47          deltastore = memo2 delta
48          delta' i j = deltastore !! i !! j
49
50splitting :: (Enum t, Fractional t) => t -> t -> Int -> [t]
51splitting a b n = [a, a + h .. b]
52    where h     = step a b n
53
54optimalSplitting :: Floating t => t -> t -> Int -> [t]
55optimalSplitting a b n = [0.5 * (b - a) * c i + (b + a) | i <- [n,n-1..0]]
56    where c i          = cos $ pi * (fromIntegral $ 2 * i + 1) / (fromIntegral $ 2 * n + 2)
57
58
59printTable :: [Double] -> [Double] -> [Double] -> [Double] -> [Double] -> IO ()
60printTable s1 s2 s3 s4 s5 = sequence_ [printRow s | s <- (zip5 s1 s2 s3 s4 s5)]
61    where printRow (c1, c2, c3, c4, c5) = printf "%13.7f | %13.10f | %13.10f | %13.10f | %13.10f\n" c1 c2 c3 c4 c5
62
63main :: IO ()
64main = sequence_ $
65       plotFuncs options xs [f1, f2, g] :
66       printf "%13s | %13s | %13s | %13s | %13s\n" "x" "y" "lagrange" "lagrange opt" "newton" :
67       printTable xs1 ys1 fs1 fs2 gs :
68       []
69    where a   = 1
70          b   = 11
71          n   = 21
72          h   = \x -> 10 `logBase` x + 7 / (2 * x + 6)
73          xs  = splitting a b (n * 10)
74          xs1 = splitting a b n
75          xs2 = optimalSplitting a b n
76          ys1 = map h xs1
77          ys2 = map h xs2
78          fs1 = map f1 xs1
79          fs2 = map f2 xs1
80          gs  = map g xs1
81          f1  = lagrange xs1 ys1 :: Double -> Double
82          f2  = lagrange xs2 ys2 :: Double -> Double
83          g   = newton a b ys1   :: Double -> Double
84          options = [Title "Interpolation", XRange (a - 1, b + 1)]