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

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

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

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

Option Explicit

'データ点の数 - 1
Private Const N = 6

Dim x(): ReDim x(N)
Dim y(): ReDim y(N)

'1.5刻みで -4.5〜4.5 まで, ７点だけ値をセット
Dim i
For i = 0 To N
    Dim d: d = i * 1.5 - 4.5
    x(i) = d
    y(i) = f(d)
Next

'0.5刻みで 与えられていない値を補間 
For i = 0 To 18
    d = i * 0.5 - 4.5
    Dim d1: d1 = f(d)
    Dim d2: d2 = lagrange(d, x, y)

    '元の関数と比較
    WScript.StdOut.Write     Right(Space(5) & FormatNumber(d,       2, -1, 0, 0), 5) & vbTab
    WScript.StdOut.Write     Right(Space(8) & FormatNumber(d1,      5, -1, 0, 0), 8) & vbTab
    WScript.StdOut.Write     Right(Space(8) & FormatNumber(d2,      5, -1, 0, 0), 8) & vbTab
    WScript.StdOut.WriteLine Right(Space(8) & FormatNumber(d1 - d2, 5, -1, 0, 0), 8)
Next

'元の関数
Private Function f(ByVal x)
    f = x - (x ^ 3) / (3 * 2) + (x ^ 5) / (5 * 4 * 3 * 2)
End Function

'Lagrange (ラグランジュ) 補間
Private Function lagrange(ByVal d, ByVal x(), ByVal y())
    Dim sum: sum = 0
    Dim i, j
    For i = 0 To N
        Dim prod: prod = y(i)
        For j = 0 To N
            If j <> i Then
                prod = prod * (d - x(j)) / (x(i) - x(j))
            End If
        Next
        sum = sum + prod
    Next
    lagrange = sum
End Function

Z:\>cscript //nologo Z:\0701.vbs
-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

参考文献