Numerical interpolation

Interpolation.hs

module Main where

import Plotting

skipNthElement :: Int -> [a] -> [a]
skipNthElement d ys = take d ys <> drop (d + 1) ys

lagrangePolynome :: [Double] -> Int -> Double -> Double
lagrangePolynome xs i x = product $ skipNthElement i $ map pn [0 .. n - 1]
  where
    pn :: Int -> Double
    pn j = (x - (xs !! j)) / ((xs !! i) - (xs !! j))
    n = length xs

interpolate :: [(Double, Double)] -> Double -> Double
interpolate points x =
  sum $ zipWith (\y i -> y * lagrangePolynome xs i x) ys [0 ..]
  where
    xs = map fst points
    ys = map snd points

genDatapointsEq :: Int -> (Double, Double) -> (Double -> Double) -> [(Double, Double)]
genDatapointsEq n interval f =
  map (\x -> (x, f x)) [lerp n' (0.0, nD) interval | n' <- [0.0 .. nD]]
  where
    nD = fromIntegral n

genDatapointsTchebychev ::
     Int -> (Double, Double) -> (Double -> Double) -> [(Double, Double)]
genDatapointsTchebychev n (a, b) f =
  map (\x -> (x, f x)) [xit x | x <- [0.0 .. fromIntegral n]]
  where
    xit :: Double -> Double
    xit i =
      (((b - a) * 0.5) * cos (((2.0 * (n' - i) + 1.0) * pi) / (2.0 * n' + 2.0))) +
      ((a + b) * 0.5)
    n' :: Double
    n' = fromIntegral n

runge :: Double -> Double
runge x = recip $ 1 + x**2

main :: IO ()
main = do
  let knotsCount = 20
  let interval = (-5.0, 5.0) :: (Double, Double)
  let valsEq = genDatapointsEq knotsCount interval runge
  let valsT = genDatapointsTchebychev knotsCount interval runge
  let x = 4.68
  let params =
        PlottingParameters
          { minX = -5.0
          , maxX = 5.0
          , minY = -5.0
          , maxY = 5.0
          , renderGrid = True
          , windowWidth = 800
          , windowHeight = 800
          , unboundVar = "x"
          , ticSize = 10
          }
  putStrLn
    "Comparing Runge-function interpolated equidistant and with Tchebychev points"
  putStrLn $ "Real        value at x = " <> show x <> ": " <> show (runge x)
  putStrLn $
    "Equidistant value at x = " <> show x <> ": " <> show (interpolate valsEq x)
  putStrLn $
    "Tchebychev  value at x = " <> show x <> ": " <> show (interpolate valsT x)
  putStrLn "Plotting: red - real, green - equal, yellow - Tchebychev"
  plotFunctions [runge, interpolate valsEq, interpolate valsT] params

Plotting.hs

module Plotting where

import Graphics.Gloss
import Text.Printf

data PlottingParameters =
    PlottingParameters
        { minX :: Double
        , maxX :: Double
        , minY :: Double
        , maxY :: Double
        , renderGrid :: Bool
        , windowWidth :: Int
        , windowHeight :: Int
        , unboundVar :: String
        , ticSize :: Int
        }

plotFunctions :: [(Double -> Double)] -> PlottingParameters -> IO ()
plotFunctions fs params = do
  let bgCol = greyN 0.2
  let paths = map (flip functionToPath params) fs
  let psOfPath = zipWith (\col -> color col . line) colors paths
  let axes = axis params
  display
    (InWindow
       ("Plot in " <> unboundVar params)
       (windowWidth params, windowHeight params)
       (100, 100))
    bgCol
    (translate
       ((* (-0.5)) $ fromIntegral $ windowWidth params)
       ((* (-0.5)) $ fromIntegral $ windowHeight params) $
     pictures $
     (if renderGrid params
        then grid params
        else blank) :
     axes : psOfPath)

colors :: [Color]
colors = cycle [red, green, yellow, blue, violet, azure, orange, magenta, white]

grid :: PlottingParameters -> Picture
grid params =
    color (withAlpha 0.7 $ greyN 0.4) $
    pictures $
    [ line
        [ (fromIntegral x * xScale, 0)
        , (fromIntegral x * xScale, fromIntegral $ windowHeight params)
        ]
    | x <- [0 .. ticSize params]
    ] <>
    [ line
        [ (0, fromIntegral y * yScale)
        , (fromIntegral (windowWidth params), fromIntegral y * yScale)
        ]
    | y <- [0 .. ticSize params]
    ]
  where
    xScale = fromIntegral $ windowWidth params `div` ticSize params
    yScale = fromIntegral $ windowHeight params `div` ticSize params

functionToPath :: (Double -> Double) -> PlottingParameters -> Path
functionToPath f params =
  map (convertXCoordToPoint f params) [0 .. windowWidth params]

convertXCoordToPoint ::
       (Double -> Double) -> PlottingParameters -> Int -> Point
convertXCoordToPoint f params xcoord =
  ( fromIntegral xcoord
  , realToFrac $
    lerp
      (f $
       lerp
         (toEnum xcoord)
         (0.0, toEnum $ windowWidth params)
         (minX params, maxX params))
      (minY params, maxY params)
      (0, (fromIntegral $ windowHeight params)))

lerp :: Double -> (Double, Double) -> (Double, Double) -> Double
lerp x (x0, x1) (y0, y1) = y0 + ((y1 - y0) / (x1 - x0)) * (x - x0)

axis :: PlottingParameters -> Picture
axis params = color white $ pictures [xAxis params, yAxis params]

yAxis :: PlottingParameters -> Picture
yAxis params =
    pictures $
    line
        [ (fromIntegral $ windowWidth params `div` 2, 0)
        , ( fromIntegral $ windowWidth params `div` 2
          , fromIntegral $ windowHeight params)
        ] :
    map (makeYLabel params . (* (windowHeight params `div` ticSize params)))
        [0 .. ticSize params]

makeYLabel :: PlottingParameters -> Int -> Picture
makeYLabel params n =
    translate (fromIntegral $ windowWidth params `div` 2) (fromIntegral n) $
    pictures [dash, translate 2 (-5) label]
  where
    dash = line [(-3, 0), (3, 0)]
    label = scale 0.1 0.1 $ text $ printf "%.2f" yCoord
    yCoord =
        lerp
            (fromIntegral n)
            (0, fromIntegral $ windowHeight params)
            (minY params, maxY params)

xAxis :: PlottingParameters -> Picture
xAxis params =
    pictures $
    line
        [ (0, fromIntegral $ windowHeight params `div` 2)
        , ( (fromIntegral $ windowWidth params)
          , fromIntegral $ windowHeight params `div` 2)
        ] :
    map (makeXLabel params . (* (windowWidth params `div` ticSize params)))
        [0 .. ticSize params]

makeXLabel :: PlottingParameters -> Int -> Picture
makeXLabel params n =
    translate (fromIntegral n) (fromIntegral $ windowHeight params `div` 2) $
    pictures [dash, translate 0 (-13) $ label]
  where
    dash = line [(0, 3), (0, -3)]
    label = scale 0.1 0.1 $ text $ printf "%.2f" xCoord
    xCoord =
        lerp
            (fromIntegral n)
            (0, fromIntegral $ windowWidth params)
            (minX params, maxX params)