{-# LANGUAGE ParallelListComp #-} module Math.Polynomial.Interpolation where import Math.Polynomial import Math.Polynomial.Lagrange import Data.List -- |Evaluate a polynomial passing through the specified set of points. The -- order of the interpolating polynomial will (at most) be one less than -- the number of points given. polyInterp :: Fractional a => [(a,a)] -> a -> a polyInterp xys = head . last . neville xys -- |Computes the tableau generated by Neville's algorithm. Each successive -- row of the table is a list of interpolants one order higher than the previous, -- using a range of input points starting at the same position in the input -- list as the interpolant's position in the output list. neville :: Fractional a => [(a,a)] -> a -> [[a]] neville xys x = table where (xs,ys) = unzip xys table = ys : [ [ ((x - x_j) * p1 + (x_i - x) * p0) / (x_i - x_j) | p0:p1:_ <- tails row | x_j <- xs | x_i <- x_is ] | row <- table | x_is <- tail (tails xs) , not (null x_is) ] -- |Computes the tableau generated by a modified form of Neville's algorithm -- described in Numerical Recipes, Ch. 3, Sec. 2, which records the differences -- between interpolants at each level. Each pair (c,d) is the amount to add -- to the previous level's interpolant at either the same or the subsequent -- position (respectively) in order to obtain the new level's interpolant. -- Mathematically, either sum yields the same value, but due to numerical -- errors they may differ slightly, and some \"paths\" through the table -- may yield more accurate final results than others. nevilleDiffs :: Fractional a => [(a,a)] -> a -> [[(a,a)]] nevilleDiffs xys x = table where (xs,ys) = unzip xys table = zip ys ys : [ [ ( {-c-} (x_j - x) * (c1 - d0) / (x_j - x_i) , {-d-} (x_i - x) * (c1 - d0) / (x_j - x_i) ) | (_c0,d0):(c1,_d1):_ <- tails row | x_j <- xs | x_i <- x_is ] | row <- table | x_is <- tail (tails xs) , not (null x_is) ] -- |Fit a polynomial to a set of points by iteratively evaluating the -- interpolated polynomial (using 'polyInterp') at 0 to establish the -- constant coefficient and reducing the polynomial by subtracting that -- coefficient from all y's and dividing by their corresponding x's. -- -- Slower than 'lagrangePolyFit' but stable under different sets of -- conditions. -- -- Note that computing the coefficients of a fitting polynomial is an -- inherently ill-conditioned problem. In most cases it is both faster and -- more accurate to use 'polyInterp' or 'nevilleDiffs' instead of evaluating -- a fitted polynomial. iterativePolyFit :: Fractional a => [(a,a)] -> Poly a iterativePolyFit = poly LE . loop where loop [] = [] loop xys = c0 : loop (drop 1 xys') where c0 = polyInterp xys 0 xys' = [ (x,(y - c0) / x) | (x,y) <- xys ] -- |Fit a polynomial to a set of points using barycentric Lagrange polynomials. -- -- Note that computing the coefficients of a fitting polynomial is an -- inherently ill-conditioned problem. In most cases it is both faster and -- more accurate to use 'polyInterp' or 'nevilleDiffs' instead of evaluating -- a fitted polynomial. lagrangePolyFit :: Fractional a => [(a,a)] -> Poly a lagrangePolyFit xys = sumPolys [ scalePoly f (fst (contractPoly p x)) | f <- zipWith (/) ys phis | x <- xs ] where (xs,ys) = unzip xys p = lagrange xs phis = map (snd . evalPolyDeriv p) xs