二次插值

二次插值,第1张

已知(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))//分段二次插值

}


欢迎分享,转载请注明来源:内存溢出

原文地址: http://outofmemory.cn/yw/11079909.html

(0)
打赏 微信扫一扫 微信扫一扫 支付宝扫一扫 支付宝扫一扫
上一篇 2023-05-13
下一篇 2023-05-13

发表评论

登录后才能评论

评论列表(0条)

保存