差商
对于一系列互不相等的 x 0 , x 1 , ⋯ , x n x_0,x_1,\cdots,x_n x0,x1,⋯,xn,有差商定义如下:
f [ x i , x j ] = f ( x i ) − f ( x j ) x i − x j ( i ≠ j , i , j = 0 , 1 , ⋯ , n ) f[x_i,x_j]=\frac{f(x_i)-f(x_j)}{x_i-x_j}(i\neq j,i,j=0,1,\cdots,n) f[xi,xj]=xi−xjf(xi)−f(xj)(i=j,i,j=0,1,⋯,n)
称此为 f ( x ) f(x) f(x)在 x i , x j x_i,x_j xi,xj处的一阶差商,又称
f ( x i , x j , x k ) = f [ x i , x j ] − f [ x j , x k ] x i − x k f(x_i,x_j,x_k)=\frac{f[x_i,x_j]-f[x_j,x_k]}{x_i-x_k} f(xi,xj,xk)=xi−xkf[xi,xj]−f[xj,xk]
为 f ( x ) f(x) f(x)在 x i , x j , x k x_i,x_j,x_k xi,xj,xk处的二阶差商。
因此可以得到n阶差商的定义:
f [ x 0 , x 1 , ⋯ , x n ] = f [ x 0 , x 1 , ⋯ , x n − 1 ] − f [ x 1 , x 2 , ⋯ , x n ] x 0 − x n f[x_0,x_1,\cdots,x_n]=\frac{f[x_0,x_1,\cdots,x_{n-1}]-f[x_1,x_2,\cdots,x_n]}{x_0-x_n} f[x0,x1,⋯,xn]=x0−xnf[x0,x1,⋯,xn−1]−f[x1,x2,⋯,xn]
差商的计算方法如下:
牛顿插值法
与拉格朗日插值法不同,如果需要增加插值点,拉格朗日插值法需要重新计算所有插值点的值,重新计算基函数,而牛顿插值法只需要重附上一项即可,不需要重新计算。
牛顿插值法需要利用到差商,其插值公式如下:
N n ( x ) = f ( x 0 ) + f [ x 0 , x 1 ] ∗ ( x − x 0 ) + f [ x 0 , x 1 , x 2 ] ∗ ( x − x 0 ) ( x − x 1 ) + f [ x 0 , x 1 , x 2 , x 3 ] ∗ ( x − x 0 ) ( x − x 1 ) ( x − x 2 ) + ⋯ + f [ x 0 , x 1 , ⋯ , x n ] ∗ ( x − x 0 ) ( x − x 1 ) ⋯ ( x − x n − 1 ) \begin{aligned} N_n(x)=&f(x_0)+f[x_0,x_1]*(x-x_0)+f[x_0,x_1,x_2]*(x-x_0)(x-x_1)+\\&f[x_0,x_1,x_2,x_3]*(x-x_0)(x-x_1)(x-x_2)+\cdots +\\&f[x_0,x_1,\cdots,x_n]*(x-x_0)(x-x_1)\cdots(x-x_{n-1}) \end{aligned} Nn(x)=f(x0)+f[x0,x1]∗(x−x0)+f[x0,x1,x2]∗(x−x0)(x−x1)+f[x0,x1,x2,x3]∗(x−x0)(x−x1)(x−x2)+⋯+f[x0,x1,⋯,xn]∗(x−x0)(x−x1)⋯(x−xn−1)
代码
//牛顿插值法
#include<iostream>
#include<cmath>
//求差商
void poor_quotient(double *xi,double *yi,double *Ni,int n)
{
double **a_temp=new double*[n];
for (int i=0; i < n;i++)
{
a_temp[i]=new double[n];
a_temp[i][0]=yi[i];
}
for (int j = 1; j < n;j++)
{
for (int i=j; i < n;i++)
{
a_temp[i][j]=(a_temp[i][j-1]-a_temp[i-1][j-1])/(xi[i]-xi[i-j]);
}
}
for (int i = 0; i < n; i++)
{
Ni[i] = a_temp[i][i];
}
for (int i = 0; i < n;i++)
{
delete[] a_temp[i];
}
delete[] a_temp;
}
double Newton(double *Ni,double *xi,double x,int n)
{
double result =Ni[0];
double temp = 1.0f;
for (int i = 1;i<n;i++)
{
temp *= x - xi[i-1];
result += Ni[i] * temp;
}
return result;
}
double fx(double x)
{
//这里填入被插值函数
//如果插值区间包括分母为0的情况需要手动调整
return sqrt(x + 3);
}
int main()
{
int n;//插值点数
std::cout<<"请输入插值点数:";
std::cin>>n;
double error[3]={0.0f,0.0f,0.0f};//误差评价
double *x = new double [n];
double *y = new double [n];
double *Ni = new double [n];
double a=3, b=10;//定义插值区间
//生成插值数据
for (int i = 0;i<n;i++)
{
x[i] = a + (b-a)*i/(n-1);
y[i] = fx(x[i]);
}
//计算差商
poor_quotient(x, y, Ni, n);
//std::cout << std::endl;
// 验证误差,抽取区间内的100个点进行误差检测
for (int i = 0; i < 100; i++)
{
double x_temp = a + (b - a) * i / (100 - 1);
double y_temp = fx(x_temp);
double y_temp2 = Newton(Ni, x, x_temp, n);
if (fabs(y_temp - y_temp2) > error[0])
{
error[0] = fabs(y_temp - y_temp2);
// std::cout << y_temp2 << std::endl;
}
error[1] += fabs((y_temp - y_temp2) / y_temp);
error[2] += ((y_temp - y_temp2)) * ((y_temp - y_temp2));
}
//err[0]得到的是在区间内最大的误差
//err[1]得到的是平均最大相对误差
//err[2]是均方根误差
error[1] = error[1]/100;
error[2] = sqrt(error[2]/100);
for(int i = 0;i < 3;i++)
{
std::cout<<"误差"<<i+1<<":"<<error[i]<<std::endl;
}
delete[] x;
delete[] y;
delete[] Ni;
return 0;
}
结果分析
与多项式插值、拉格朗日插值相比,对于相同的测试函数,其具有更高的精度,并且其龙格现象也更加不明显。