C++で関数の近似(ネヴィル補間)

をネヴィル補間で近似する

与えられた５個の関数値 f(x₀), f(x₁), f(x₂), f(x₃), f(x₄) と、未知の点を通る４次式を求める場合, 次のような順序で計算する.
(1)

(2)

(3)

(4)

(5)

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

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

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

// 元の関数
double f(double x);

// Neville (ネヴィル) 補間
double neville(double d, double x[], double y[]);

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 = neville(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;
}

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

// Neville (ネヴィル) 補間
double neville(double d, double x[], double y[])
{
    static double w[N][N];

    //w[0][0] = y[0];
    //w[1][1] = y[1];
    //w[2][2] = y[2];
    //w[3][3] = y[3];
    //w[4][4] = y[4];

    for (int i = 0; i < N; i++)
        w[i][i] = y[i];

    //w[1][0] = w[1][1] + (w[1][1] - w[0][0]) * (t - x[1]) / (x[1] - x[0]);
    //w[2][1] = w[2][2] + (w[2][2] - w[1][1]) * (t - x[2]) / (x[2] - x[1]);
    //w[3][2] = w[3][3] + (w[3][3] - w[2][2]) * (t - x[3]) / (x[3] - x[2]);
    //w[4][3] = w[4][4] + (w[4][4] - w[3][3]) * (t - x[4]) / (x[4] - x[3]);

    //w[2][0] = w[2][1] + (w[2][1] - w[1][0]) * (t - x[2]) / (x[2] - x[0]);
    //w[3][1] = w[3][2] + (w[3][2] - w[2][1]) * (t - x[3]) / (x[3] - x[1]);
    //w[4][2] = w[4][3] + (w[4][3] - w[3][2]) * (t - x[4]) / (x[4] - x[2]);

    //w[3][0] = w[3][1] + (w[3][1] - w[2][0]) * (t - x[3]) / (x[3] - x[0]);
    //w[4][1] = w[4][2] + (w[4][2] - w[3][1]) * (t - x[4]) / (x[4] - x[1]);

    //w[4][0] = w[4][1] + (w[4][1] - w[3][0]) * (t - x[4]) / (x[4] - x[0]);

    for (int j = 1; j < N; j++)
    {
        for (int i = j; i < N; i++)
            w[i][i-j] = w[i][i-j+1] + (w[i][i-j+1] - w[i-1][i-j]) * (d - x[i]) / (x[i] - x[i-j]);
    }
    
    return w[N-1][0];
}

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

Z:\>CP0703
-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

ONLY DO WHAT ONLY YOU CAN DO

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

参考文献