Haskellで関数の近似(ラグランジュ補間)

をラグランジュ補間で近似する

n+1個の点 (x₀, y₀), (x₁, y₁) … (x_n, y_n)
が与えられているとき, これらすべての点を通る n次式は次のように表すことができる.

この式を使って, 与えられた点以外の点の値を求める.

import Text.Printf
import Control.Monad

-- データ点の数 - 1
n = 6 :: Int

-- 元の関数
f::Double->Double
f x = x - (x ^ 3) / (fromIntegral (3 * 2)) + (x ^ 5) / (fromIntegral (5 * 4 * 3 * 2))

-- Lagrange (ラグランジュ) 補間
lagrange::Double->[Double]->[Double]->Double
lagrange d x y = 
    let sum_list = map(\i -> do
        let prod_list = map(\j -> do
            (d - x!!(fromIntegral j)) / 
                (x!!(fromIntegral i)  - 
                 x!!(fromIntegral j))) 
                $ filter (\j -> j /= i) [0..n]
        product $ y!!(fromIntegral i) : prod_list
        ) [0..n]
    in
        sum $ sum_list

main = do
    -- 1.5刻みで -4.5〜4.5 まで, ７点だけ値をセット
    let x = map(\i -> (fromIntegral i) * 1.5 - 4.5) [0..n]
    let y = map(\i -> f(i)) x

    -- 0.5刻みで 与えられていない値を補間
    let d1 = map(\i -> (fromIntegral i) * 0.5 - 4.5) [0..18]
    let d2 = map(\i -> (f i)) d1
    let d3 = map(\i -> (lagrange i x y)) d1

    forM_ (zip (zip d1 d2) d3) $ \((d1, d2), d3) -> do
        printf "%5.2f\t%8.5f\t%8.5f\t%8.5f\n" d1 d2 d3 (d2 - d3)

Z:\>runghc Hs0701.hs
-4.50   -4.68984        -4.68984         0.00000
-4.00   -1.86667        -1.86667        -0.00000
-3.50   -0.73099        -0.73099        -0.00000
-3.00   -0.52500        -0.52500         0.00000
-2.50   -0.70964        -0.70964        -0.00000
-2.00   -0.93333        -0.93333        -0.00000
-1.50   -1.00078        -1.00078         0.00000
-1.00   -0.84167        -0.84167        -0.00000
-0.50   -0.47943        -0.47943        -0.00000
 0.00    0.00000         0.00000         0.00000
 0.50    0.47943         0.47943         0.00000
 1.00    0.84167         0.84167         0.00000
 1.50    1.00078         1.00078         0.00000
 2.00    0.93333         0.93333         0.00000
 2.50    0.70964         0.70964         0.00000
 3.00    0.52500         0.52500         0.00000
 3.50    0.73099         0.73099         0.00000
 4.00    1.86667         1.86667         0.00000
 4.50    4.68984         4.68984         0.00000

参考文献