插值算法.docx

插值算法.docx

《插值算法.docx》由会员分享，可在线阅读，更多相关《插值算法.docx（26页珍藏版）》请在冰豆网上搜索。

插值算法

一插值算法简介:

1:

插值的涵义:

在离散数据的基础上补插连续函数，使得这条连续曲线通过全部给定的离散数据点。

插值是离散函数逼近的重要方法，利用它可通过函数在有限个点处的取值状况，估算出函数在其他点处的近似值。

早在6世纪，中国的刘焯已将等距二次插值用于天文计算。

17世纪之后，I.牛顿，J.-L.拉格朗日分别讨论了等距和非等距的一般插值公式。

在近代，插值法仍然是数据处理和编制函数表的常用工具，又是数值积分、数值微分、非线性方程求根和微分方程数值解法的重要基础，许多求解计算公式都是以插值为基础导出的。

插值问题的提法是：

假定区间[a，b］上的实值函数f（x）在该区间上n＋1个互不相同点x0，x1……xn处的值是f[x0］，……f（xn），要求估算f（x）在[a，b］中某点的值。

其做法是：

在事先选定的一个由简单函数构成的有n＋1个参数C0，C1，……Cn的函数类Φ（C0，C1，……Cn）中求出满足条件P（xi）＝f（xi）（i＝0，1，……n）的函数P（x），并以P（）作为f（）的估值。

此处f（x）称为被插值函数，c0，x1，……xn称为插值结（节）点，Φ（C0，C1，……Cn）称为插值函数类，上面等式称为插值条件，Φ（C0，……Cn）中满足上式的函数称为插值函数，R（x）＝f（x）－P（x）称为插值余项。

当估算点属于包含x0，x1……xn的最小闭区间时，相应的插值称为内插，否则称为外插。

2:

插值的种类

（1）多项式插值

这是最常见的一种函数插值。

在一般插值问题中，若选取Φ为n次多项式类，由插值条件可以唯一确定一个n次插值多项式满足上述条件。

从几何上看可以理解为：

已知平面上n＋1个不同点，要寻找一条n次多项式曲线通过这些点。

插值多项式一般有两种常见的表达形式，一个是拉格朗日插值多项式，另一个是牛顿插值多项式。

（2）埃尔米特插值

对于函数f（x），常常不仅知道它在一些点的函数值，而且还知道它在这些点的导数值。

这时的插值函数P（x），自然不仅要求在这些点等于f（x）的函数值，而且要求P（x）的导数在这些点也等于f（x）的导数值。

这就是埃尔米特插值问题，也称带导数的插值问题。

从几何上看，这种插值要寻求的多项式曲线不仅要通过平面上的已知点组，而且在这些点（或者其中一部分）与原曲线“密切”，即它们有相同的斜率。

可见埃尔米特插值多项式比起一般多项式插值有较高的光滑逼近要求。

（3）分段插值与样条插值

为了避免高次插值可能出现的大幅度波动现象，在实际应用中通常采用分段低次插值来提高近似程度，比如可用分段线性插值或分段三次埃尔米特插值来逼近已知函数，但它们的总体光滑性较差。

为了克服这一缺点，一种全局化的分段插值方法——三次样条插值成为比较理想的工具。

见样条函数。

3.1分段插值

在每个区间

上，用1阶多项式（直线）逼近f（x）:

即用折线代替曲线。

易证：

当

优点：

计算简单；

适用于光滑性要求不高的插值问题。

分段三次（Hermite）插值

不少实际插值问题不仅要求函数值相等，而且还要求导数值也相等。

这就导致下面的Hermite插值。

给定

（4）三角函数插值

当被插函数是以2π为周期的函数时，通常用n阶三角多项式作为插值函数，并通过高斯三角插值表出。

（5）其它提法

插值（Interpolation），有时也称为“重置样本”，是在不生成像素的情况下增加图像像素大小的一种方法，在周围像素色彩的基础上用数学公式计算丢失像素的色彩。

有些相机使用插值，人为地增加图像的分辨率。

插值：

用来填充图像变换时像素之间的空隙。

（6）三次样条插值:

上面介绍的分段线性插值,其总体光滑程度不够.在数学上,光滑程度的定量描述是：

函数（曲线）的k阶导数存在且连续,则称该曲线具有k阶光滑性.自然,阶数越高光滑程度越好.于是,分段线性插值具有零阶光滑性,也就是不光滑；分段三次埃尔米特插值具有一阶光滑性.仅有这些光滑程度,在工程设计和机械加工等实际中是不够的.提高分段函数如多项式函数的次数,可望提高整体曲线的光滑程度.但是,是否存在较低次多项式达到较高阶光滑性的方法？

三次样条插值就是一个很好的例子.

样条曲线本身就来源于飞机、船舶等外形曲线设计中所用的绘图工具.在工程实际中,要求这样的曲线应该具有连续的曲率,也就是连续的二阶导数.值得注意的是分段插值曲线的光滑性关键在于段与段之间的衔接点（节点）处的光滑性.

三次样条函数记为s（x）,它是定义在区间[a,b]上的函数,满足:

1）s（x）在每一个小区间[Xi-1,Xi]上是一个三次多项式函数；

2）在整个区间[a,b]上,其二阶导数存在且连续.即在每个节点处的二阶导数连续.

三次样条插值问题的提法:

给定函数f（x）在n+1个节点x0,x1,...,xn处的函数值为y0,y1,...,yn,求一个三次样条函数s（x）,使其满足：

s（xi）=yi,i=0,1,…,n.

如何确定三次样条函数在每一个小区间上的三次多项式函数的系数呢？

这是一个比较复杂的问题,这里只介绍确定系数的思想.

分段线性插值在每一段的线性函数的两个参数,是由两个方程（两个端点处的函数值为给定值）唯一确定；对于三次样条插值呢,每一个区间上的三次函数有四个参数,而在该区间上由两个端点的函数值为给定值只能够产生两个方程,仅此不足以唯一确定四个参数.注意到三次样条函数对整体光滑性的要求,其二阶导数存在且连续,从全局的角度上考虑参数个数与方程个数的关系如下：

参数：

每个小段上4个,n个小段共计4n个.

方程：

每个小段上由给定函数值得到2个,n个小段共计2n个；

光滑性要求每一个内部节点处的一阶、二阶导数连续,得出其左右导数相等,因此,每个节点产生2个方程,共计2（n-1）个.

现在得到了4n-2个方程,还差两个.为此,常用的方法是对边界节点除函数值外附加要求.这就是所谓的边界条件.需要两个,正好左右两个端点各一个.常用如下三类边界条件.m边界条件：

s'（X0）=m0,s'（Xn）=mn.即两个边界节点的一阶导数值为给定值：

m0,mn.

M边界条件：

s''（x0）=m0,s''（xn）=mn.即两个边界节点的二阶导数值为给定值：

m0,mn.特别地,当m0和mn都为零时,称为自然边界条件.

周期性边界条件：

s'（x0）=s'（xn）;s''（x0）=s''（xn）.

以上分析说明,理论上三次样条插值函数是确定的,具体如何操作,可以查阅有关资料.

例题，观测数据为

x=[012345678910];

y=[020-4040-2031];

待求的三次多项式函数s（x）在[010]上有连续的一阶，二阶导数.我们通过简单的讨论来认识问题。

在第一区间[01]、第二区间[12]上考虑两个三次多项式

s（x）=s1*x^3+s2*x^2+s3*x+s4

r（x）=r1*x^3+r2*x^2+r3*x+r4

示意图:

可以得到

s（0）=s1*0^3+s2*0^2+s3*0+s4=0

（1）

s

（1）=s1*1^3+s2*1^2+s3*1+s4=2

（2）

r

（1）=r1*2^3+r2*2^2+r3*2+r4=2（3）

r

（2）=r1*2^3+r2*2^2+r3*2+r4=0（4）

一阶导函数

s'（x）=3*s1*x^2+2*s2*x+s3

r'（x）=3*r1*x^2+2*r2*x+r3

由一阶导数的连续性且在1点处相等，有

3*s1*1^2+2*s2*1+s3=3*r1*1^2+2*r2*1+r3（5）

二阶导函数

s''（x）=6*s1*x+2*s2

r''（x）=6*r1*x+2*r2

由二阶导数的连续性且在1点处相等，有

6*s1*1+2*s2=6*r1*1+2*r2（6）

由m边界条件

s'（0）=1.6,r'

（2）=0.3

有

3*s1*0^2+2*s2*0+s3=1.6（7）

3*r1*2^2+2*r2*2+r3=0.3（8）

M边界条件

s''（0）=-1,r''

（2）=1

有

6*s1*0+2*s2=-1（7'）

6*r1*2+2*r2=1（8'）

由周期性边界条件

s'（0）=r'

（2）

s''（0）=r''

（2）

有

3*s1*0^2+2*s2*0+s3=-1（7''）

3*r1*2^2+2*r2*2+r3=1（8''）这样，对于两个多项式的8个未知量，我们给出了8个方程。

三次样条曲线的难点在于，我们不能分段去求解方程，完成绘图。

//三次样条插值函数代码：

/*

注:

所用的数据表如下:

x0701302103375787761012114214621841

y05778103135182214244256272275

边界条件为:

（1）y'（0）=1,y'（1841）=0;

（2）y"（0）=0,y"（1841）=0.

求解在以上边界条件下在x[k]=50k（k=1,2...,36）处的函数值程序源代码*/

二:

插值算法的代码实现：

1三次样条的代码实现及测试结果

//求三次样条插值函数

#include

#include

usingnamespacestd;

constintMAX=50;

floatx[MAX],y[MAX],h[MAX];//变量设置：

x为各点横坐标；y为各点纵坐标；h为步长

floatc[MAX],a[MAX],fxym[MAX];

floatf（intx1,intx2,intx3）/*****************求差分函数（含三个参数）****************************/

{

floata=（y[x3]-y[x2]）/（x[x3]-x[x2]）;

floatb=（y[x2]-y[x1]）/（x[x2]-x[x1]）;

return（a-b）/（x[x3]-x[x1]）;

}

voidcal_m（intn）/***********************用追赶法求解出弯矩向量M……***************************/

{

floatB[MAX];

B[0]=c[0]/2;

for（inti=1;i

B[i]=c[i]/（2-a[i]*B[i-1]）;

//fxym[0]=fxym[0]/2;

for（inti=1;i<=n;i++）

fxym[i]=（fxym[i]-a[i]*fxym[i-1]）/（2-a[i]*B[i-1]）;

for（inti=n-1;i>=0;i--）

fxym[i]=fxym[i]-B[i]*fxym[i+1];

}

voidprintout（intn）;

intmain（）

{

intn,i;charch;

do

{

cout<<"请输入已知断点个数:

";

cin>>n;

for（i=0;i<=n;i++）

{

cout<<"PleaseputinX"<

';

cin>>x[i];//cout<

cout<<"PleaseputinY"<

';

cin>>y[i];//cout<

}

for（i=0;i

h[i]=x[i+1]-x[i];

cout<<"Please输入边界条件\n1:

已知两端的一阶导数\n2:

两端的二阶导数已知\n默认:

自然边界条件\n";

intt;

floatf0,f1;

cin>>t;

switch（t）

{

case1:

cout<<"PleaseputinY0\'Y"<

cin>>f0>>f1;

c[0]=1;a[n]=1;

fxym[0]=6*（（y[1]-y[0]）/（x[1]-x[0]）-f0）/h[0];

fxym[n]=6*（f1-（y[n]-y[n-1]）/（x[n]-x[n-1]））/h[n-1];

break;

case2:

cout<<"PleaseputinY0\"Y"<

cin>>f0>>f1;

c[0]=a[n]=0;

fxym[0]=2*f0;fxym[n]=2*f1;

break;

default:

cout<<"不可用\n";//待定

};//switch

for（i=1;i

fxym[i]=6*f（i-1,i,i+1）;//调用差分函数（only!

）

for（i=1;i

{

a[i]=h[i-1]/（h[i]+h[i-1]）;

c[i]=1-a[i];

}

a[n]=h[n-1]/（h[n-1]+h[n]）;

cal_m（n）;//调用弯矩函数（only！

）

cout<<"\n输出三次样条插值函数：

\n";

printout（n）;//调用求解三次样条插值函数；函数输出

cout<<"Doyoutohaveanthertry?

y/n:

";

cin>>ch;

}

while（ch=='y'||ch=='Y'）;

return0;

}

voidprintout（intn）/***************求三次样条插值函数（因已知断点个数而异）***********************/

{

cout<所包含；括号内为参数。

for（inti=0;i

{

cout<

["<

floatt=fxym[i]/（6*h[i]）;

if（t>0）

cout<

else

cout<<-t<<"*（x-"<

t=fxym[i+1]/（6*h[i]）;

if（t>0）

cout<<"+"<

else

cout<<"-"<<-t<<"*（x-"<

cout<<"\n\t";

t=（y[i]-fxym[i]*h[i]*h[i]/6）/h[i];

if（t>0）

cout<<"+"<

else

cout<<"-"<<-t<<"*（"<

t=（y[i+1]-fxym[i+1]*h[i]*h[i]/6）/h[i];

if（t>0）

cout<<"+"<

else

cout<<"-"<<-t<<"*（x-"<

cout<

}

cout<

}

/****************提供测试数据：

（来自课本页例《数值分析》清华大学出版社第四版）***********/

/*输入

3

27.74.1

284.3

294.1

303.0

1

3.0-4.0*/

/*输出

输出三次样条插值函数：

1:

[27.7,28]

13.07*（x-28）^3+0.22*（x-27.7）^3

+14.84*（28-x）+14.31*（x-27.7）

2:

[28,29]

0.066*（29-x）^3+0.1383*（x-28）^3

+4.234*（29-x）+3.962*（x-28）

3:

[29,30]

0.1383*（30-x）^3-1.519*（x-29）^3

+3.962*（30-x）+4.519*（x-29）

*/

/*利用三次样条插值法求曲线在某个点处的函数值（本题用第二个边界条件）*/

＃include

＃include

#defineN11//N表示节点的个数

voidmain（）

{

doublex[N]={0,70,130,210,337,578,776,1012,1142,1462,1841},

y[N]={0,57,78,103,135,182,214,244,256,272,275},

h[N],a[N],b[N],A[N],B[N],m[N],s[37],xx[37];//s[N]表示曲线,其中数组xx[N]表示自变量x[k]

inti,k;

h[0]=x[1]-x[0];//初始化h0,a0,b0,A0,B0

a[0]=1;

b[0]=3*（y[1]-y[0]）/h[0];

A[0]=-a[0]/2;

B[0]=b[0]/2;

for（i=0;i

h[i]=x[i+1]-x[i];

for（i=1;i

a[i]=h[i-1]/（h[i-1]+h[i]）;

b[i]=3*（（1-a[i]）*（y[i]-y[i-1]）/h[i-1]+a[i]*（y[i+1]-y[i]）/h[i]）;

}

for（i=1;i

A[i]=-a[i]/（2+（1-a[i]）*A[i-1]）;

B[i]=（b[i]-（1-a[i]）*B[i-1]）/（2+（1-a[i]）*A[i-1]）;

}

m[N-1]=（b[N-1]-（1-a[N-1]）*B[N-2]）/（2+（1-a[N-1]）*A[N-2]）;//求Mn的值

for（i=N-2;i>=0;i--）//求m0,m1,-----mn-1的值

m[i]=A[i]*m[i+1]+B[i];

for（k=1;k<=36;k++）{//求曲线在x[k]处的函数值

xx[k]=50*k;

for（i=0;i

if（xx[k]>=x[i]&&xx[k]<=x[i+1]）

//s[k]即为x[k]所在区间的三次样条插值函数,以下即为求在x[k]处的函数值

s[k]=（1+2*（xx[k]-x[i]）/（x[i+1]-x[i]））*pow（（xx[k]-x[i+1]）/（x[i]-x[i+1]）,2）*y[i]+

（1+2*（xx[k]-x[i+1]）/（x[i]-x[i+1]））*pow（（xx[k]-x[i]）/（x[i+1]-x[i]）,2）*y[i+1]+

（xx[k]-x[i]）*pow（（xx[k]-x[i+1]）/（x[i]-x[i+1]）,2）*m[i]+

（xx[k]-x[i+1]）*pow（（xx[k]-x[i]）/（x[i+1]-x[i]）,2）*m[i+1];

}

printf（"所求的外形曲线在x[k]=50*k（k=1,2,...,36）处的函数值分别为:

\n"）;

for（k=1;k<=36;k++）//输出在x[k]处的函数值

printf（"s（%4d）=%13.8f\n",50*k,s[k]）;

printf（"\n"）;

printf（"m[i]的值分别为:

\n"）;

for（i=0;i

printf（"m（%2d）=%8f\n",i,m[i]）;

printf（"\n"）;

}

//程序代码实现2：

#include

#include

#include

usingnamespacestd;

#definen4//三次样条插值分段数=数据的组数减一

voidchase（doublea[n-1][n-1],doublef[n-1],doublem[n+1]）

{

//三对角阵的Crout分解

doubleL[n-1][n-1];

doubleU[n-1][n-1];

doubley[n-1];

doublex[n-1];

intt,r=0;

for（inti=0;i

{

for（intj=0;j

{

L[i][j]=U[i][j]=0;

}

}

for（inti=0;i

{

if（i==0）

{

if（a[i][i]==0||a[i][i+1]==0）

{

cout<<"a["<

break;

}

if（fabs（a[i][i]）<=fabs（a[i][i+1]））

{

cout<<"a["<

break;

}

}

elseif（i==n-1-1）

{

if（a[i][i]==0||a[i][i-1]==0）

{

cout<<"a["<

break;

}

if（fabs（a[i][i]）<=fabs（a[i][i-1]））

{

cout<<"a["<

break;

}

}

else

{

if（a[i][i]==0||a[i][i-1]==0||a[i][i+1]==0）

{

cout<<"a["<