C++で関数の近似(ネヴィル補間)
与えられた5個の関数値 f(x0), f(x1), f(x2), f(x3), f(x4) と、 未知の点を通る4次式を求める場合, 次のような順序で計算する.
(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 まで, 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