二次插值

已知（xi，yi），i＝0，1，2是三个互异点，二次插值就是求过这三个已知点的抛物线方程，即构造一个二次代数多项式L2（x）：

地球物理数据处理基础

使其满足：

地球物理数据处理基础

显然，满足条件（6－10）的插值多项式函数（6－9）可由下面的方程组确定：

地球物理数据处理基础

此方程中a0，a1，a2为未知数，写成矩阵形式为

地球物理数据处理基础

于是，只要从式（6－12）中求解出a0，a1，a2，则插值多项式L2（x）便唯一确定。

但是，通过解方程组的办法来求取插值多项式通常是比较麻烦的，尤其是当插值节点增多时，计算量大。因此，我们有必要寻找构造插值函数的其他途径，下面我们介绍构造二次插值函数的Lagrange差值法。

根据插值问题的唯一性，参照线性插值函数形式，令

地球物理数据处理基础

li（x）（i＝0，1，2）应该具有相似的形式，由于插值函数L2（x）是二次多项式，那么li（x）（i＝0，1，2）也应该是二次多项式形式。不妨令

地球物理数据处理基础

由插值条件式（6－10）得

地球物理数据处理基础

于是得到二次插值函数：

地球物理数据处理基础

其中，li（x）（i＝0，1，2）表示如下：

地球物理数据处理基础

li（x）称为二次插值的基函数或形函数。li（x）满足

地球物理数据处理基础

同样，对于n＋1个已知观测点，我们仍可以先将观测区间［x0，xn］划分为n/2（n为偶数）段，在相邻的三个观测点构成的子区间［xi－1，xi＋1］（i＝1，3，…，n－1）进行分段二次插值。

不难得到插值函数

地球物理数据处理基础

其中：

地球物理数据处理基础

因此，对于区间［xi－1，xi＋1］内的任意点x，只要求出相应的lj（x），即可根据式（6－18）计算出x处的近似值。

#include<stdio.h>

#include<math.h>

double Lagrange1(double *x, double *y, double xx) //拉格郎日插值

{

int i,j

double *a,yy=0.000

a=new double[6]

for(i=0i<6i++)

{

a[i]=y[i]

for(j=0j<6j++)

if(j!=i)

a[i]*=(xx-x[j])/(x[i]-x[j])

yy+=a[i]

}

delete a

return yy

}

double Lagrange2(double *x, double *y, double input) //分段线性插值

{

double output

int i

for (i=0i<5i++)

{

if (x[i] <= input &&x[i+1] >= input)

{

output=y[i] +(y[i+1]-y[i])*(input-x[i])/(x[i+1]-x[i])

break

}

}

return output

}

double Lagrange3(double *x,double *y,double u) //分段二次插值

{

int i,k=0

double v

for(i=0i<6i++)

{

if(u<x[1])

{

k=0

v=y[k]*(u-x[k+1])*(u-x[k+2])/((x[k]-x[k+1])*(x[k]-x[k+2]))+y[k+1]*(u-x[k])*(u-x[k+2])/((x[k+1]-x[k])*(x[k+1]-x[k+2]))+y[k+2]*(u-x[k])*(u-x[k+1])/((x[k+2]-x[k])*(x[k+2]-x[k+1]))

}

if((x[i]<u&&u<=x[i+1])&&(fabs(u-x[i])<=fabs(u-x[i+1])))

{

k=i-1

v=y[k]*(u-x[k+1])*(u-x[k+2])/((x[k]-x[k+1])*(x[k]-x[k+2]))+y[k+1]*(u-x[k])*(u-x[k+2])/((x[k+1]-x[k])*(x[k+1]-x[k+2]))+y[k+2]*(u-x[k])*(u-x[k+1])/((x[k+2]-x[k])*(x[k+2]-x[k+1]))

}

if ((x[i]<u&&u<=x[i+1])&&fabs(u-x[i])>fabs(u-x[i+1]))

{

k=i

v=y[k]*(u-x[k+1])*(u-x[k+2])/((x[k]-x[k+1])*(x[k]-x[k+2]))+y[k+1]*(u-x[k])*(u-x[k+2])/((x[k+1]-x[k])*(x[k+1]-x[k+2]))+y[k+2]*(u-x[k])*(u-x[k+1])/((x[k+2]-x[k])*(x[k+2]-x[k+1]))

}

if(u>x[4])

{

k=3

v=y[k]*(u-x[k+1])*(u-x[k+2])/((x[k]-x[k+1])*(x[k]-x[k+2]))+y[k+1]*(u-x[k])*(u-x[k+2])/((x[k+1]-x[k])*(x[k+1]-x[k+2]))+y[k+2]*(u-x[k])*(u-x[k+1])/((x[k+2]-x[k])*(x[k+2]-x[k+1]))

}

}

return v

}

void main()

{

double x[6] = {0.0, 0.1, 0.195, 0.3, 0.401, 0.5},y[6] = {0.39894,0.39695,0.39142,0.38138,0.36812,0.35206}

double u

scanf("%lf",&u)

printf("%f\n",Lagrange1(x,y,u))//拉格郎日插值

printf("%f\n",Lagrange2(x,y,u))//分段线性插值

printf("%f\n",Lagrange3(x,y,u))//分段二次插值

}

---