C++で関数の近似(エルミート補間)
をエルミート補間で近似する
２個の関数値 f(x0), f(x1) と,それぞれの微分係数f'(x0), f'(x1) とを
与えられたとき, 与えられた２個の関数値を通る５３次式を求めるには, まず次のような表を作る.
このとき,
.
あとは, ニュートン補間と同じ
この式を使って, 与えられた点以外の点の値を求める.
#include <iostream> #include <iomanip> #include <math.h> using namespace std; // データ点の数 const int N = 3; const int N2 = 6; // 元の関数 double f(double x); // 導関数 double fd(double x); // Hermite (エルミート) 補間 double hermite(double d, double z[], double a[]); int main() { double x[N]; double y[N]; double yy[N]; // 1.5刻みで -1.5〜1.5 まで, ３点だけ値をセット for (int i = 0; i < N; i++) { double d = i * 1.5 - 1.5; x[i] = d; y[i] = f(d); yy[i] = fd(d); } // 差分商の表を作る double z[N2]; //z[0] = x[0]; //z[1] = x[0]; //z[2] = x[1]; //z[3] = x[1]; //z[4] = x[2]; //z[5] = x[2]; double d[N2][N2]; //d[0][0] = y[0]; //d[0][1] = y[0]; //d[0][2] = y[1]; //d[0][3] = y[1]; //d[0][4] = y[2]; //d[0][5] = y[2]; for (int i = 0; i < N2; i++) { int j = i / 2; z[i] = x[j]; d[0][i] = y[j]; } //d[1][0] = yy[0]; //d[1][1] = (d[0][2] - d[0][1]) / (z[2] - z[1]); //d[1][2] = yy[1]; //d[1][3] = (d[0][4] - d[0][3]) / (z[4] - z[3]); //d[1][5] = yy[2]; //d[2][0] = (d[1][1] - d[1][0]) / (z[2] - z[0]); //d[2][1] = (d[1][2] - d[1][1]) / (z[3] - z[1]); //d[2][2] = (d[1][3] - d[1][2]) / (z[4] - z[2]); //d[2][3] = (d[1][4] - d[1][3]) / (z[5] - z[3]); //d[3][0] = (d[2][1] - d[2][0]) / (z[3] - z[0]); //d[3][1] = (d[2][2] - d[2][1]) / (z[4] - z[1]); //d[3][2] = (d[2][3] - d[2][2]) / (z[5] - z[2]); //d[4][0] = (d[3][1] - d[3][0]) / (z[4] - z[0]); //d[4][1] = (d[3][2] - d[3][1]) / (z[5] - z[1]); //d[5][0] = (d[4][1] - d[4][0]) / (z[5] - z[0]); for (int i = 1; i < N2; i++) { for (int j = 0; j < N2 - i; j++) { if (i == 1 && j % 2 == 0) d[i][j] = yy[j / 2]; else d[i][j] = (d[i-1][j+1] - d[i-1][j]) / (z[j+i] - z[j]); } } // ｎ階差分商 double a[N2]; //a[0] = d[0][0]; //a[1] = d[1][0]; //a[2] = d[2][0]; //a[3] = d[3][0]; //a[4] = d[4][0]; //a[5] = d[5][0]; for (int j = 0; j < N2; j++) a[j] = d[j][0]; // 0.1刻みで 与えられていない値を補間 for (int i = 0; i <= 80; i++) { double d = i * 0.1 - 4.0; double d1 = hermite(d, z, a); 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(9) << fixed << setprecision(5) << d1 - d2 << endl; } return 0; } // 元の関数 double f(double x) { return x - pow(x,3)/6 + pow(x,5)/120; } // 導関数 double fd(double x) { return (24 - 12 * pow(x,2) + pow(x,4)) / 24; } // Hermite (エルミート) 補間 double hermite(double d, double z[], double a[]) { double sum = a[0]; //sum += a[1] * (d - z[0]); //sum += a[2] * (d - z[0]) * (d - z[1]); //sum += a[3] * (d - z[0]) * (d - z[1]) * (d - z[2]); //sum += a[4] * (d - z[0]) * (d - z[1]) * (d - z[2]) * (d - z[3]); //sum += a[5] * (d - z[0]) * (d - z[1]) * (d - z[2]) * (d - z[3]) * (d - z[4]); for (int i = 1; i < N2; i++) { double prod = 1; for (int j = 0; j < i; j++) prod *= (d - z[j]); sum += a[i] * prod; } return sum; }
Z:\>bcc32 CP0704.cpp Borland C++ 5.5.1 for Win32 Copyright (c) 1993, 2000 Borland cp0704.cpp: Turbo Incremental Link 5.00 Copyright (c) 1997, 2000 Borland Z:\>CP0704 -4.00: -1.86667 - -1.86667 = -0.00000 -3.90: -1.53218 - -1.53218 = -0.00000 -3.80: -1.25760 - -1.25760 = -0.00000 -3.70: -1.03650 - -1.03650 = -0.00000 -3.60: -0.86285 - -0.86285 = -0.00000 -3.50: -0.73099 - -0.73099 = -0.00000 -3.40: -0.63562 - -0.63562 = -0.00000 -3.30: -0.57178 - -0.57178 = -0.00000 -3.20: -0.53487 - -0.53487 = -0.00000 -3.10: -0.52060 - -0.52060 = -0.00000 -3.00: -0.52500 - -0.52500 = -0.00000 -2.90: -0.54443 - -0.54443 = -0.00000 -2.80: -0.57553 - -0.57553 = -0.00000 -2.70: -0.61524 - -0.61524 = -0.00000 -2.60: -0.66078 - -0.66078 = -0.00000 -2.50: -0.70964 - -0.70964 = -0.00000 -2.40: -0.75955 - -0.75955 = -0.00000 -2.30: -0.80853 - -0.80853 = -0.00000 -2.20: -0.85480 - -0.85480 = -0.00000 -2.10: -0.89684 - -0.89684 = -0.00000 -2.00: -0.93333 - -0.93333 = -0.00000 以下略
エルミート補間に３個の関数値と,それぞれの微分係数を与えた結果
ちなみに、ニュートン補間に６点を与えた場合...
