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

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

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

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

#include <iostream>
#include <iomanip>
#include <math.h>
using namespace std;

const int  N = 5; // データ点の数

// 元の関数
double f(double x)
{
    return 4 * pow(x,4) + 3 * pow(x,3) - 2 * pow(x,2) - x + 1;
}

// Lagrange (ラグランジュ) 補間
double lagrange(double d, double x[], double y[])
{
    //double nume[N];
    //nume[0] =                 (d    - x[1]) * (d    - x[2]) * (d    - x[3]) * (d    - x[4]);
    //nume[1] = (d    - x[0]) *                 (d    - x[2]) * (d    - x[3]) * (d    - x[4]);
    //nume[2] = (d    - x[0]) * (d    - x[1]) *                 (d    - x[3]) * (d    - x[4]);
    //nume[3] = (d    - x[0]) * (d    - x[1]) * (d    - x[2]) *                 (d    - x[4]);
    //nume[4] = (d    - x[0]) * (d    - x[1]) * (d    - x[2]) * (d    - x[3])                ;

    //double deno[N];
    //deno[0] =                 (x[0] - x[1]) * (x[0] - x[2]) * (x[0] - x[3]) * (x[0] - x[4]);
    //deno[1] = (x[1] - x[0]) *                 (x[1] - x[2]) * (x[1] - x[3]) * (x[1] - x[4]);
    //deno[2] = (x[2] - x[0]) * (x[2] - x[1])                 * (x[2] - x[3]) * (x[2] - x[4]);
    //deno[3] = (x[3] - x[0]) * (x[3] - x[1]) * (x[3] - x[2]) *                 (x[3] - x[4]);
    //deno[4] = (x[4] - x[0]) * (x[4] - x[1]) * (x[4] - x[2]) * (x[4] - x[3])                ;

    double sum = 0;
    //sum += y[0] * nume[0] / deno[0];
    //sum += y[1] * nume[1] / deno[1];
    //sum += y[2] * nume[2] / deno[2];
    //sum += y[3] * nume[3] / deno[3];
    //sum += y[4] * nume[4] / deno[4];

    for (int i = 0; i < N; i++)
    {
        double prod = y[i];
        for (int j = 0; j < N; j++)
        {
            if (j != i)
                prod *= (d - x[j]) / (x[i] - x[j]);
        }
        sum += prod;
    }
    return sum;
}

int main()
{
    double x[N];
    double y[N];

    // 0.25刻みで -0.5〜0.5 まで, ５点だけ値をセット
    for (int i = 0; i < N; i++)
    {
        double d = i * 0.25 - 0.5;
        x[i] = d;
        y[i] = f(d);
    }

    // 0.05刻みで 与えられていない値を補間 
    for (int i = 0; i <= 20; i++)
    {
        double d  = i * 0.05 - 0.5;
        double d1 = lagrange(d, x, y);
        double d2 = f(d);

        // 元の関数と比較
        cout << setw(5) << fixed << setprecision(2) << d       << ": "; 
        cout << setw(8) << fixed << setprecision(5) << d1      << " - "; 
        cout << setw(8) << fixed << setprecision(5) << d2      << " = "; 
        cout << setw(8) << fixed << setprecision(5) << d1 - d2 << endl; 
    }
    
   return 0;
}

Z:\>bcc32 CP0701.cpp
Borland C++ 5.5.1 for Win32 Copyright (c) 1993, 2000 Borland
cp0701.cpp:
Turbo Incremental Link 5.00 Copyright (c) 1997, 2000 Borland

Z:\>CP0701
-0.50:  0.87500 -  0.87500 =  0.00000
-0.45:  0.93565 -  0.93565 =  0.00000
-0.40:  0.99040 -  0.99040 =  0.00000
-0.35:  1.03640 -  1.03640 =  0.00000
-0.30:  1.07140 -  1.07140 =  0.00000
-0.25:  1.09375 -  1.09375 =  0.00000
-0.20:  1.10240 -  1.10240 = -0.00000
-0.15:  1.09690 -  1.09690 =  0.00000
-0.10:  1.07740 -  1.07740 =  0.00000
-0.05:  1.04465 -  1.04465 = -0.00000
 0.00:  1.00000 -  1.00000 = -0.00000
 0.05:  0.94540 -  0.94540 =  0.00000
 0.10:  0.88340 -  0.88340 =  0.00000
 0.15:  0.81715 -  0.81715 = -0.00000
 0.20:  0.75040 -  0.75040 =  0.00000
 0.25:  0.68750 -  0.68750 =  0.00000
 0.30:  0.63340 -  0.63340 = -0.00000
 0.35:  0.59365 -  0.59365 = -0.00000
 0.40:  0.57440 -  0.57440 =  0.00000
 0.45:  0.58240 -  0.58240 = -0.00000
 0.50:  0.62500 -  0.62500 =  0.00000

参考文献