ONLY DO WHAT ONLY YOU CAN DO

こけたら立ちなはれ 立ったら歩きなはれ

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



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

n+1個の点 (x0, y0), (x1, y1) … (xn, yn)
が与えられているとき, これらすべての点を通る n次式は次のように表すことができる.

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

module Fs0701

// データ点の数 - 1
let N = 6 

// 元の関数
let f (x:double):double =
    x - System.Math.Pow(x,3.0) / (float (3 * 2)) + System.Math.Pow(x,5.0) / (float (5 * 4 * 3 * 2))

// Lagrange (ラグランジュ) 補間
let lagrange(d:double) (x:double list) (y:double list) =
    let mutable sum_list = [0.0]
    for i in [0..N] do
        let mutable prod_list = [y.[i]]
        for j in [0..N] do
            if i <> j then
                prod_list <- ((d - x.[j]) / (x.[i] - x.[j]))::prod_list
        sum_list <- (prod_list |> List.reduce(*))::sum_list
    sum_list |> List.sum

// 1.5刻みで -4.5〜4.5 まで, 7点だけ値をセット
let x = [0..N] |> List.map(fun i -> (float i) * 1.5 - 4.5)
let y = x |> List.map(f)

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

(List.zip (List.zip d1 d2) d3)
|> List.iter (fun ((d1, d2), d3) -> printfn "%5.2f\t%8.5f\t%8.5f\t%8.5f" d1 d2 d3 (d2 - d3))

exit 0
Z:\>fsi Fs0701.fs --nologo --quiet
-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
参考文献