插值算法Word文件下载.docx

1、不少实际插值问题不仅要求函数值相等，而且还要求导数值也相等。这就导致下面的Hermite插值。给定（4）三角函数插值 当被插函数是以2为周期的函数时，通常用n阶三角多项式作为插值函数，并通过高斯三角插值表出。（5）其它提法插值（Interpolation），有时也称为“重置样本”，是在不生成像素的情况下增加图像像素大小的一种方法，在周围像素色彩的基础上用数学公式计算丢失像素的色彩。有些相机使用插值，人为地增加图像的分辨率。插值：用来填充图像变换时像素之间的空隙。（6）三次样条插值:上面介绍的分段线性插值,其总体光滑程度不够.在数学上,光滑程度的定量描述是：函数（曲线）的k阶导数存在且连续,则称

2、该曲线具有k阶光滑性.自然,阶数越高光滑程度越好. 于是,分段线性插值具有零阶光滑性,也就是不光滑；分段三次埃尔米特插值具有一阶光滑性.仅有这些光滑程度,在工程设计和机械加工等实际中是不够的.提高分段函数如多项式函数的次数,可望提高整体曲线的光滑程度. 但是,是否存在较低次多项式达到较高阶光滑性的方法？三次样条插值就是一个很好的例子. 样条曲线本身就来源于飞机、船舶等外形曲线设计中所用的绘图工具.在工程实际中,要求这样的曲线应该具有连续的曲率,也就是连续的二阶导数.值得注意的是分段插值曲线的光滑性关键在于段与段之间的衔接点（节点）处的光滑性. 三次样条函数记为s（x）,它是定义在区间a,b上的

3、函数,满足: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、,而在该区间上由两个端点的函数值为给定值只能够产生两个方程,仅此不足以唯一确定四个参数.注意到三次样条函数对整体光滑性的要求,其二阶导数存在且连续,从全局的角度上考虑参数个数与方程个数的关系如下：参数：每个小段上4个,n个小段共计4n个. 方程：每个小段上由给定函数值得到2个,n个小段共计2n个；光滑性要求每一个内部节点处的一阶、二阶导数连续,得出其左右导数相等,因此,每个节点产生2个方程,共计2（n-1）个. 现在得到了4n-2个方程,还差两个.为此,常用的方法是对边界节点除函数值外附加要求.这就是所谓的边界条件.需要两个,正好左右两个端点各一个.常用如下三类边界条件. m边界条件：s（X0

5、）=m0,s（Xn）=mn.即两个边界节点的一阶导数值为给定值：m0,mn. M边界条件：（x0）=m0, s（xn）=mn.即两个边界节点的二阶导数值为给定值：m0,mn. 特别地,当m0和mn都为零时,称为自然边界条件. 周期性边界条件：（x0）=s（xn）; 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

6、 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）一阶导函数（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

7、）二阶导函数（x）=6*s1*x+2*s2（x）=6*r1*x+2*r2由二阶导数的连续性且在1点处相等，有 6*s1*1+2*s2=6*r1*1+2*r2 （6） 由m边界条件（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边界条件（0）=-1,r（2）=1有 6*s1*0+2*s2=-1 （7） 6*r1*2+2*r2=1 （8由周期性边界条件（0）=r（2） 3*s1*02+2*s2*0+s3=-1 （7） 3*r1*22+2*r2*2+r3=1 （8）这样，对于两个多项式的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三次样条的代码实现及测试结果/求三次样条插值函数#includeiomanipusing namespace

9、 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 =

10、c0 / 2; for（int i = 1; i = 0; i-） fxymi = fxymi - Bi*fxymi+1;void printout（int n）;int main（） int n,i; char ch; do coutn; for（i = 0; coutPlease put in Xixi; /coutendl;Please put in Yyi; for（i = 0; i+） /求步长；其数组值较之x，y个数少一 hi = xi+1 - xi;Please 输入边界条件n 1: 已知两端的一阶导数n 2:两端的二阶导数已知n 默认:自然边界条件n int t; float

11、f0, f1;t; switch（t） case 1:coutPlease put in Y0 Ynf1; 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:Please put in Y0/显示数据为Y0至Yn，即断点的二阶导数 c0 = an = 0; fxym0 = 2*f0; fxymn = 2*f1; default:不可用n/待定 ;/switch for（i = 1; fxymi =

12、6 * f（i-1, i, i+1）;/调用差分函数（only! ai = hi-1 / （hi + hi-1）; ci = 1 - ai; an = hn-1 / （hn-1 + hn）; cal_m（n）;/调用弯矩函数（only！n输出三次样条插值函数： printout（n）;/调用求解三次样条插值函数；函数输出Do you to have anther try ? y/n :ch; while（ch = y | ch = Y）; return 0;void printout（int n）/*求三次样条插值函数（因已知断点个数而异）*/ coutsetprecision（6）;/通过操

13、作器setprecision（）设置有效位数；其为头文件所包含；括号内为参数。 for（int i = 0; i+）/所输出函数个数由所设断点个数而定 i+1 xi , xi+1 0）t*（ - x）3 else -t*（x - ）3 t = fxymi+1/（6*hi）; + - nt t = （yi - fxymi*hi*hi/6）/hi;+ - x） else - t = （yi+1 - fxymi+1*hi*hi/6）/hi;）endl/*提供测试数据：（来自课本页例 数值分析清华大学出版社第四版）*/*输入327.7 4.128 4.329 4.130 3.013.0 -4.0*/*

14、输出输出三次样条插值函数： 27.7 , 2813.07*（x - 28）3 + 0.22*（x - 27.7）3+ 14.84*（28 - x） + 14.31*（x - 27.7） 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 ncludemath.h#def

15、ine 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表示自变量xk int i,k; h0 = x1-x0; /初始化h0,a0,b0,A0,B0 a0 = 1; b0 = 3*（y1-y0）/h0; A0 = -a0/2; B0 = b0/2; for（i=0;N;i+） /求hi hi=xi+1-xi

16、; for（i=1;N-1;i+） /求ai,bi ai=hi-1/（hi-1+hi）; bi = 3 * （1-ai）*（yi-yi-1）/hi-1 + ai*（yi+1-yi）/hi）;i+） /求Ai,Bi Ai = -ai/（2+（1-ai）*Ai-1）; Bi = （bi-（1-ai）*Bi-1）/（2+（1-ai）*Ai-1）; mN-1=（bN-1-（1-aN-1）*BN-2）/（2+（1-aN-1）*AN-2）;/求Mn的值 for（i=N-2;i=0;i-） /求m0,m1,-mn-1的值 mi = Ai * mi+1 + Bi; for（k=1;k= xi & xxk =

17、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）处的函数值分别为:k+） /输出在xk处的函

18、数值 printf（s（%4d） = %13.8fn,50*k,sk）; printf（mi的值分别为:i+）m（%2d） = %8fn,i,mi）;/程序代码实现2：#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;n-1; for（int j=0;jj+） Lij=Uij=0;i+）/判断是否符合Crout分解定理 if（i=0） if（aii=0|aii+1=0） coutan-1可约 break; if（fabs（aii）=fabs（aii+1）强超 break; else if（i=n-1-1） if（aii=0|aii-1=0）=fabs（aii-1） if（aii=0|aii-1=0|aii+1=0） cout