插值算法.doc

一 插值算法简介:1:插值的涵义:在离散数据的基础上补插连续函数，使得这条连续曲线通过全部给定的离散数据点。插值是离散函数逼近的重要方法，利用它可通过函数在有限个点处的取值状况，估算出函数在其他点处的近似值。 早在6世纪，中国的刘焯已将等距二次插值用于天文计算。17世纪之后，I.牛顿，J.-L.拉格朗日分别讨论了等距和非等距的一般插值公式。在近代，插值法仍然是数据处理和编制函数表的常用工具，又是数值积分、数值微分、非线性方程求根和微分方程数值解法的重要基础，许多求解计算公式都是以插值为基础导出的。 插值问题的提法是：假定区间a，b上的实值函数f（x）在该区间上 n1个互不相同点x0，x1xn 处的值是f x0，f（xn），要求估算f（x）在a，b中某点的值。其做法是：在事先选定的一个由简单函数构成的有n1个参数C0，C1，Cn的函数类(C0，C1，Cn)中求出满足条件P（xi）f（xi）（i0，1， n）的函数P(x)，并以P()作为f()的估值。此处f（x）称为被插值函数，c0，x1，xn称为插值结(节)点，(C0，C1，Cn)称为插值函数类，上面等式称为插值条件，（C0，Cn）中满足上式的函数称为插值函数，R（x） f（x）P（x）称为插值余项。当估算点属于包含x0，x1xn的最小闭区间时，相应的插值称为内插，否则称为外插。2:插值的种类(1)多项式插值 这是最常见的一种函数插值。在一般插值问题中，若选取为n次多项式类，由插值条件可以唯一确定一个n次插值多项式满足上述条件。从几何上看可以理解为：已知平面上n1个不同点，要寻找一条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=0 1 2 3 4 5 6 7 8 9 10; y=0 2 0 -4 0 4 0 -2 0 3 1;待求的三次多项式函数s(x)在0 10上有连续的一阶，二阶导数.我们通过简单的讨论来认识问题。在第一区间0 1、第二区间1 2上考虑两个三次多项式 s(x)=s1*x3+s2*x2+s3*x+s4 r(x)=r1*x3+r2*x2+r3*x+r4示意图:可以得到 s(0)=s1*03+s2*02+s3*0+s4=0 (1) s(1)=s1*13+s2*12+s3*1+s4=2 (2) r(1)=r1*23+r2*22+r3*2+r4=2 (3) r(2)=r1*23+r2*22+r3*2+r4=0 (4)一阶导函数 s(x)=3*s1*x2+2*s2*x+s3 r(x)=3*r1*x2+2*r2*x+r3由一阶导数的连续性且在1点处相等，有 3*s1*12+2*s2*1+s3=3*r1*12+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*02+2*s2*0+s3=1.6 (7) 3*r1*22+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*02+2*s2*0+s3=-1 (7) 3*r1*22+2*r2*2+r3=1 (8)这样，对于两个多项式的8个未知量，我们给出了8个方程。三次样条曲线的难点在于，我们不能分段去求解方程，完成绘图。/三次样条插值函数代码：/*注:所用的数据表如下:x 0 70 130 210 337 578 776 1012 1142 1462 1841y 0 57 78 103 135 182 214 244 256 272 275边界条件为:(1)y(0)=1,y(1841)=0;(2)y(0)=0,y(1841)=0.求解在以上边界条件下在xk=50k(k=1,2.,36)处的函数值程序源代码*/二:插值算法的代码实现：1三次样条的代码实现及测试结果/求三次样条插值函数#include#includeusing namespace std;const int MAX = 50;float xMAX, yMAX, hMAX;/变量设置：x为各点横坐标；y为各点纵坐标；h为步长float cMAX, aMAX, fxymMAX;float f(int x1, int x2, int x3)/*求差分函数(含三个参数)*/float a = (yx3 - yx2) / (xx3 - xx2);float b = (yx2 - yx1) / (xx2 - xx1);return (a - b)/(xx3 - xx1); void cal_m(int n)/*用追赶法求解出弯矩向量M*/ float BMAX;B0 = c0 / 2;for(int i = 1; i n; i+)Bi = ci / (2 - ai*Bi-1);/fxym0 = fxym0 / 2;for(int i = 1; i = 0; i-)fxymi = fxymi - Bi*fxymi+1;void printout(int n);int main()int n,i; char ch;docoutn;for(i = 0; i = n; i+)coutPlease put in Xixi; /coutendl;coutPlease put in Yiyi; /coutendl;for(i = 0; i n; i+) /求步长；其数组值较之x，y个数少一hi = xi+1 - xi;coutt;switch(t)case 1:coutPlease put in Y0 Ynf0f1;c0 = 1; an = 1;fxym0 = 6*(y1 - y0) / (x1 - x0) - f0) / h0;fxymn = 6*(f1 - (yn - yn-1) / (xn - xn-1) / hn-1;break;case 2:coutPlease put in Y0 Ynf0f1;c0 = an = 0;fxym0 = 2*f0; fxymn = 2*f1;break;default:cout不可用n;/待定;/switchfor(i = 1; i n; i+)fxymi = 6 * f(i-1, i, i+1);/调用差分函数(only!)for(i = 1; i n; i+)ai = hi-1 / (hi + hi-1);ci = 1 - ai;an = hn-1 / (hn-1 + hn);cal_m(n);/调用弯矩函数（only！）coutn输出三次样条插值函数：n;printout(n);/调用求解三次样条插值函数；函数输出coutch;while(ch = y | ch = Y);return 0;void printout(int n)/*求三次样条插值函数（因已知断点个数而异）*/coutsetprecision(6);/通过操作器setprecision()设置有效位数；其为头文件所包含；括号内为参数。for(int i = 0; i n; i+)/所输出函数个数由所设断点个数而定couti+1: xi , xi+1n 0)coutt*(xi+1 - x)3;else cout-t*(x - xi+1 0)cout + t*(x - xi)3;else cout - -t*(x - xi)3;cout 0)cout+ t*(xi+1 - x);else cout- -t*(xi+1 0)cout + t*(x - xi);else cout - -t*(x - xi);coutendlendl;coutendl;/*提供测试数据：（来自课本页例 数值分析清华大学出版社第四版）*/*输入327.7 4.128 4.329 4.130 3.013.0 -4.0*/*输出输出三次样条插值函数：1: 27.7 , 2813.07*(x - 28)3 + 0.22*(x - 27.7)3+ 14.84*(28 - x) + 14.31*(x - 27.7)2: 28 , 290.066*(29 - x)3 + 0.1383*(x - 28)3+ 4.234*(29 - x) + 3.962*(x - 28)3: 29 , 300.1383*(30 - x)3 - 1.519*(x - 29)3+ 3.962*(30 - x) + 4.519*(x - 29)*/*利用三次样条插值法求曲线在某个点处的函数值(本题用第二个边界条件)*/i ncludei nclude#define N 11 /N表示节点的个数void main()double xN=0,70,130,210,337,578,776,1012,1142,1462,1841,yN=0,57,78,103,135,182,214,244,256,272,275,hN,aN,bN,AN,BN,mN,s37,xx37; /sN表示曲线,其中数组xxN表示自变量xkint i,k; h0 = x1-x0; /初始化h0,a0,b0,A0,B0a0 = 1;b0 = 3*(y1-y0)/h0;A0 = -a0/2;B0 = b0/2;for(i=0;iN;i+) /求hihi=xi+1-xi;for(i=1;iN-1;i+) /求ai,biai=hi-1/(hi-1+hi); bi = 3 * (1-ai)*(yi-yi-1)/hi-1 + ai*(yi+1-yi)/hi);for(i=1;i=0;i-) /求m0,m1,-mn-1的值mi = Ai * mi+1 + Bi;for(k=1;k=36;k+) /求曲线在xk处的函数值xxk = 50 * k;for(i=0;i= xi & xxk = xi+1)/sk即为xk所在区间的三次样条插值函数,以下即为求在xk处的函数值sk = (1+2*(xxk-xi)/(xi+1-xi)*pow(xxk-xi+1)/(xi-xi+1),2)*yi +(1+2*(xxk-xi+1)/(xi-xi+1)*pow(xxk-xi)/(xi+1-xi),2)*yi+1 +(xxk-xi)*pow(xxk-xi+1)/(xi-xi+1),2)*mi +(xxk-xi+1)*pow(xxk-xi)/(xi+1-xi),2)*mi+1 ; printf(所求的外形曲线在xk=50*k(k=1,2,.,36)处的函数值分别为:n);for(k=1;k=36;k+) /输出在xk处的函数值printf(s(%4d) = %13.8fn,50*k,sk); printf(n);printf(mi的值分别为:n);for(i=0;iN;i+)printf(m(%2d) = %8fn,i,mi);printf(n);/程序代码实现2：#include#include#includeusing namespace std; #define n 4 /三次样条插值分段数=数据的组数减一void chase(double an-1n-1,double fn-1,double mn+1) /三对角阵的Crout分解double Ln-1n-1;double Un-1n-1;double yn-1;double xn-1;int t,r=0;for(int i=0;in-1;i+)for(int j=0;jn-1;j+)Lij=Uij=0;for(int i=0;in-1;i+)/判断是否符合Crout分解定理if(i=0)if(aii=0|aii+1=0)coutan-1n-1可约endl;break;if(fabs(aii)=fabs(aii+1)coutan-1n-1强超endl;break;else if(i=n-1-1)if(aii=0|aii-1=0)coutan-1n-1可约endl;break;if(fabs(aii)=fabs(aii-1)coutan-1n-1强超endl;break;else if(aii=0|aii-1=0|aii+1=0)coutan-1n-1可约endl;break;if(fabs(aii)fabs(aii-1)+fabs(aii+1)coutan-1n-1强超endl;break;t=i;if(t=n-1-1)coutendl该三角方程系数矩阵存在唯一的两个三角阵endl;int i=0;for(int i=0;in-1;i+)if(i=0)Lii=aii;Uii+1=aii+1/Lii;Uii=1;continue;if(i=n-1-1)Lii-1=Lii=aii-aii-1*Ui-1i;Uii=1;elseLii-1=Lii=aii-aii-1*Ui-1i;Uii+1=aii+1/Lii;Uii=1;coutendl下三角阵L为:endl;for(int i=0;in-1;i+)for(int j=0;jn-1;j+)coutsetw(10)Lij; /设定输出宽度coutendl;coutendl上三角阵U为:endl;for(int i=0;in-1;i+)for(int j=0;jn-1;j+)coutsetw(10)Uij;coutendl;y0=f0/a00; /解下三角方程组for(int i=1;i=0;i-)xi=yi-Uii+1*xi+1;for(int i=0;in-1;i+)mi+1=xi; void main()double xn+1=0.25,0.30,0.39,0.45,0.53;double yn+1=0.5000,0.5477,0.6245,0.6708,0.7280;double dy0=1,dyn=0.6868;double h_in;double wn-1;double vn-1;double dn-1;double An-1n-1;double bn-1;double m(n-1)+2; /有n+1个单元int i,j;for(i=0;in-1;i+)bi=0;for(j=0;jn-1;j+)if(i=j)Aii=2;elseAij=0;cout第一边界条件下样条插值算法endl;cout插值数据：endlsetw(8)xi;for(i=0;in+1;i+)coutsetw(8)xi;coutendlsetw(8)yi;for(i=0;in+1;i+)coutsetw(8)yi;for(i=0;in;i+)h_ii=xi+1-xi;for(i=0;in-1;i+)wi=h_ii+1/(h_ii+h_ii+1);/miuvi=h_ii/(h_ii+h_ii+1);/na m dadi=3*(wi*(yi+1-yi)/h_ii+vi*(yi+2-yi+1)/h_ii+1);for(i=0;in-2;i+)Aii+1=vi;Ai+1i=wi+1;for(i=0;in-1;i+)bi=di;if(i=0)bi-=wi*dy0;else if(i=n-2)bi-=vi*dyn;coutendl解得对角方程组endlA=endl;for(i=0;in-1;i+)for(j=0;jn-1;j+)coutsetw(10)Aij;coutendl;coutendlb=endl;for(i=0;in-1;i+)coutsetw(10)bi;coutendlendl使用“追赶法”解三角方程组：;chase(A,b,m);cout解得：endl;for(i=1;in;i+)coutsetw(10)mi=miendl;m0=dy0;mn=dyn;/*for(i=0;in+1;i+)coutsetw(10)mi=miendl;*/double s=0; /具有局限性，分多少段就定义多少个变量/分段函数系数double ss4n; /第i+1行为-i次方次数，i=0,1,2,3double xx;int flag=0,tap;char ans;/y or nfor(i=0;in;i+)ss0i=2*(yi-yi+1)/pow(xi+1-xi),3)+(mi+mi+1)/pow(xi+1-xi),2);ss1i=3*(xi+1+xi)*(yi-yi+1)/pow(xi-xi+1),3)-(2*xi+1+xi)*mi+(2*xi+xi+1)*mi+1)/pow(xi+1-xi),2);ss2i=6*xi+1*xi*(yi-yi+1)/pow(xi+1-xi),3)+(pow(xi+1,2)+2*xi+1*xi)*mi+(pow(xi,2)+2*xi*xi+1)*mi+1)/pow(xi+1-xi),2);ss3i=(pow(xi+1,2)*(xi+1-3*xi)*yi-pow(xi,2)*(xi-3*xi+1)*yi+1)/pow(xi+1-xi),3)-xi+1*xi*(xi+1*mi+xi*mi+1)/pow(xi-xi+1),2);cout三次样条插值函数：endl;couts(x)=endl;for(i=0;in;i+)coutsi+1(x)=ss0ix3;cout.setf(ios:showpos); /强迫显示正负号coutss1ix2ss2ixss3isetw(15);coutresetiosflags(ios:showpos); /除去正负号格式coutx-xi,xi+1endl;loop:cout请输入x属于x0,xn，并在x上插值：xx;/注意输入这个数据时不能是字母if(xxxn)flag=0;/xx不在定义域内for(i=1;in+1;i+)if(xx=xi)flag=1;/tap=i;break;if(flag=1)cout用第tap个方程：endl;for(i=0;i4;i+)s+=ssitap-1*pow(xx,3-i);/ss4ncouts(tap)=sendl;s=0;flag=0;cout若继续插值，请输入y,否则输入任意键结束:ans;/字母或数字都可以if(ans=121|ans=89)goto loop;2线性分段插值代码实现/*拉格朗日和牛顿插值算法实现 */#include#includeusing namespace std;typedef strUCt datafloat x;float y;Data;/变量x和函数值y的结构Data d20;/最多二十组数据float f(int s,int t)/牛顿插值法，用以返回插商if(t=s+1)return (dt.y-ds.y)/(dt.x-ds.x);elsereturn (f(s+1,t)-f(s,t-1)/(dt.x-ds.x); float Newton(float x,int count)int n;while(1)coutn;if(n=count-1)/ 插值次数不得大于count次break;elsesystem(cls);/初始化t，y，yt。float t=1.0;float y=d0.y;float yt=0.0;/计算y值for(int j=1;j=n;j+)t=(x-dj-1.x)*t;yt=f(0,j)*t;/coutf(0,j)endl;y=y+yt;return y;float lagrange(float x,int count)float y=0.0;for(int k=0;kcount;k+)/这儿默认为count次插值float p=1.0;/初始化pfor(int j=0;jcount;j+)/计算p的值if(k=j)