数值分析课程设计--三次样条插值.doc

数值分析课程设计三次样条插值算法院（系）名称 信息工程学院 专 业 班 级 09普本信计1班 学 号 090111073 学 生 姓 名 宣章然 指 导 教 师 孔繁民 2012年06月08日 数值分析 课程设计评阅书题目三次样条插值学生姓名宣章然学号090111073指导教师评语及成绩指导教师签名： 年 月 日答辩评语及成绩答辩教师签名： 年 月 日教研室意见 总成绩： 教研室主任签名：年 月 日课程设计任务书20082009学年第二学期专业班级： 09普本信计1班 学号： 060111060 姓名： 宣章然 课程设计名称： 数值分析 设计题目： 三次样条插值 完成期限：自 2012 年 6 月 8 日至 2012 年 6 月 13 日共 1 周设计依据、要求及主要内容：一、设计目的 熟练掌握三次样条插值算法的原理和推导过程，并且能够应用Matlab软件编写相应的程序和使用Matlab软件函数库软件。 二、设计要求 (1)用Matlab函数库中相应函数对选定的问题，求出具有一定精度的结果。 (2)使用所用的方法编写Matlab程序求解,对数值结果进行分析。 (3)对于使用多个方法解同一问题的，在界面上设计成菜单形式。 三、设计内容 首先构造三次样条插值函数的定义和一般特征,并对实例问题进行实例分析,并总结 四、参考文献 1 黄明游,冯果忱.数值分析M.北京:高等教育出版社,2008. 2 马东升,雷勇军.数值计算方法M.北京:机械工业出版社,2006. 3 石博强,赵金.MATLAB数学计算与工程分析范例教程M.北京:中国铁道出版社.2005. 4郝红伟，MATLAB 6，北京，中国电力出版社，2001 5姜健飞，胡良剑，数值分析及其MATLAB实验， 科学出版社，2004 6薛毅，数值分析实验，北京工业大学出版社，2005 计划答辩时间：2012年6月18日指导教师（签字）： 教研室主任（签字）： 批准日期： 年 月课程设计说明书（论文） 第XIX页三次样条插值摘 要分段低次样条插值虽然计算简单、稳定性好、收敛性有保证且易在电子计算机上实现，但只能保证各小段曲线在连接处的连续性，不能保证整件曲线的光滑性。利用样条插值，既可保持分段低次插值多项式，又可提高插值函数光滑性。故给出分段三次样条插值的构造过程、算法步骤，利用MATLAB软件编写三次样条插值函数通用程序，并通过数值算例证明程序的正确性。关键字:三次样条 插值函数 MATLAB编程 收敛性 算法步骤一 三次样条函数定义及特征定义1：若函数,且在每个小区间上上是三次多项式，其中 是给定节点，则称是节点上的三次样条函数。若节点上 给定函数值 ,且 （1.1）成立，则称 为三次样条差值函数。从定义知，要求出，在每个应小区间 上确定4个待定系数，共有 n个小区间，故应确定4n 个参数，根据 在 上二阶导数连续，在节点处应满足连续性条件 （1.2）共有 3n-3个条件，再加上满足插值条件（1.1），共有4n-2个条件，因此还需要2个条件才能确定。通常可在区间 端点上各加一个条件（称边界条件），边界条件可根据实际的问题要求给定。常见的三种：(1) 已知两端的一节导数值，即 （1.3）(2)两端的二阶导数已知，即 （1.4）特殊情况下的边界条件 （1.4） 称为自然边界条件（3）当是以 为周期函数时，则要求 也是周期函数，这时边界条件应满足sx0+0=sxn-0,sx0+0=sxn-0sx0+0=sxn-0,而此时式中y0=yn , 这样确定的样条函数sx 称为周期函数。二 函数推导原理及构造我们采用待定一阶导数的方法即设S(Xj)=Mj,j=0,1,.,n,因为分段三次Hermite插值多项式已经至少是一阶连续可导了，为了让它成为三次样条函数只需确定节点处的一阶导数使这些节点处的二阶导数连续即可！由于在内部节点处二阶导数连续条件：整理化简后得： 第一类三次样条插值问题方程组由于已知：基本方程组化为n-1阶方程组化为矩阵形式这是一个严格对角占优的三对角方程组,用追赶法可以求解！第二类三次样条插值问题的方程组,由于已知：故得：稍加整理得联合基本方程组得一个n+1阶三对角方程组，化成矩阵形式为：仍然是严格对角占优第三类样条插值问题的方程组，由于：立即可得下式：其中:联合基本方程得一个广义三对角或周期三对角方程组：求解这些不同类型的样条插值问题的方程组，我们可得所要待定的一阶导数：再代入S(x)的每一段表达式，就求得三次样条函数的表达式!利用插值（即求过已知有限个数据点的近似函数）的基本原理，用多项式作为研究插值的工具，进行代数插值。其基本问题是：已知函数f (x)在区间a,b上n +1个不同点x0，xn处的函数值 (i = 0,1,n)，求一个至多n 次多项式n(x)使其在给定点处与 f (x)同值，即满足插值条件: n(x)= = .许多工程技术中提出的计算问题对插值函数的光滑性有较高要求，如飞机的机翼外形，内燃机的进、排气门的凸轮曲线，都要求曲线具有较高的光滑程度，不仅要连续，而且要有连续的曲率，这就导致了样条插值的产生。 数学上将具有一定光滑性的分段多项式称为样条函数。具体地说，给定区间a,b的一个分划如果函数s(x) 满足：（i）在每个小区间 (i=0,1,n)上s(x)是k 次多项式；（ii）s(x)在a,b上具有k 1阶连续导数。则称s(x)为关于分划 的k 次样条函数，其图形称为k 次样条曲线。由于三次样条插值我、函数s(x)的插值节点处的二阶导数存在，因此令各节点处的二阶导数为 （1.01） 根据样条插值函数的定义，三次样条插值函数是s(x)在每一个小区间 上市不超过三次的多项式。在每一个小区间上，其二阶导数为线性函数，即 （1.02）对式（1.02）积分两次，则得到 （1.03）其中为任意常数。又根据样条插值函数定义中的条件（3），即可以确定与为= （1.04） 将式（1.04）中与的值代入表达式（1.03后，就可以得到样条插值函数在区间上的表达式为 （1.05）其中与分别为区间两端点处的二阶导数值。由此可以看处，只要能确定各点处的二阶导数值，则子渠道间上的三次样条插值函数也确定了。在区间a,b上的一阶导数连续，在各节点的左右两子区间上的s(x)虽然不同，但在连接点处的导数存在，即在连接点处的左右导数相等，有 （1.06）为了利用条件（2.18），在x属于时，县求为 （1.07）当x属于时， （1.08）整理得： （1.09）其中三 问题的提出上面讨论的分段低次插值函数都有一致收敛性，但光滑性较差，对于像高速飞机的机翼形线，船体放样等型值线往往要求有二阶光滑度，即有二阶连续导数，早期工程师制图时，把富有弹性的细长木条（所谓样条）用压铁固定在样点上，在其他地方让它自由弯曲，然后画下长条的曲线，称为样条曲线。它实际上是由分段三次曲线并接而成，在连接点即样点上要求二阶导数连续，从数学上加以概括就得到数学样条这一概念。下面我们讨论最常用的三次样条函数。四 实例应分析函数的MATLAB的程序设计例1：已知一组数据点，编写一程序求解三次样条插值函数满足 并针对下面一组具体实验数据0.250.30.390.450.530.50000.54770.62450.67080.7280求解，其中边界条件为.1）三次样条插值自然边界条件源程序：function s=spline3(x,y,dy1,dyn)%x为节点，y为节点函数值，dy1，dyn分别为x=0.25,0.53处的二阶导m=length(x);n=length(y);if m=n error(x or y输入有误) returnendh=zeros(1,n-1);h(n-1)=x(n)-x(n-1);for k=1:n-2 h(k)=x(k+1)-x(k); v(k)=h(k+1)/(h(k+1)+h(k); u(k)=1-v(k);endg(1)=3*(y(2)-y(1)/h(1)-h(1)/2*dy1;g(n)=3*(y(n)-y(n-1)/h(n-1)+h(n-1)/2*dyn;for i=2:n-1 g(i)=3*(u(i-1)*(y(i+1)-y(i)/h(i)+v(i-1)*(y(i)-y(i-1)/h(i-1);endfor i=2:n-1; A(i,i-1)=v(i-1); A(i,i+1)=u(i-1);endA(n,n-1)=1;A(1,2)=1;A=A+2*eye(n);M=zhuigf(A,g); %调用函数，追赶法求Mfprintf(三次样条（三对角）插值的函数表达式n);syms X;for k=1:n-1 fprintf(S%d-%d:n,k,k+1); s(k)=(h(k)+2*(X-x(k)./h(k).3.*(X-x(k+1).2.*y(k). +(h(k)-2*(X-x(k+1)./h(k).3.*(X-x(k).2.*y(k+1). +(X-x(k).*(X-x(k+1).2./h(k).2*M(k)+(X-x(k+1).*. (X-x(k).2./h(k).2*M(k+1);ends=s.;s=vpa(s,4); %画三次样条插值函数图像for i=1:n-1 X=x(i):0.01:x(i+1);st=(h(i)+2*(X-x(i)./(h(i)3).*(X-x(i+1).2.*y(i).+(h(i)-2.*(X-x(i+1)./(h(i)3).*(X-x(i).2.*y(i+1).+(X-x(i).*(X-x(i+1).2./h(i)2*M(i)+(X-x(i+1).*.(X-x(i).2./h(i)2*M(i+1); plot(x,y,o,X,st); hold on End plot(x,y); grid on%调用的函数：%追赶法function M=zhuigf(A,g)n=length(A);L=eye(n);U=zeros(n);for i=1:n-1 U(i,i+1)=A(i,i+1);endU(1,1)=A(1,1);for i=2:n L(i,i-1)=A(i,i-1)/U(i-1,i-1); U(i,i)=A(i,i)-L(i,i-1)*A(i-1,i);endY(1)=g(1);for i=2:n Y(i)=g(i)-L(i,i-1)*Y(i-1);endM(n)=Y(n)/U(n,n);for i=n-1:-1:1 M(i)=(Y(i)-A(i,i+1)*M(i+1)/U(i,i);end2)在命令窗口输入x，y,dy1,dyn,得到三次样条函数：x=0.25,0.3,0.39,0.45,0.53;y=0.5,0.5477,0.6245,0.6708,0.7280;dy1=0;dyn=0;s=spline3(x,y,dy1,dyn)运行结果：三次样条（三对角）插值的函数表达式S1-2:S2-3:S3-4:S4-5: .5000*(-3600.+.1600e5*X)*(X-.3000)2+.5477*(5200.-.1600e5*X)*(X-.2500)2+394.8*(X-.2500)*(X-.3000)2+355.2*(X-.3000)*(X-.2500)2 .5477*(-699.6+2743.*X)*(X-.3900)2+.6245*(1193.-2743.*X)*(X-.3000)2+109.6*(X-.3000)*(X-.3900)2+96.77*(X-.3900)*(X-.3000)2 .6245*(-3333.+9259.*X)*(X-.4500)2+.6708*(4444.-9259.*X)*(X-.3900)2+217.7*(X-.3900)*(X-.4500)2+207.6*(X-.4500)*(X-.3900)2 .6708*(-1602.+3906.*X)*(X-.5300)2+.7280*(2227.-3906.*X)*(X-.4500)2+116.8*(X-.4500)*(X-.5300)2+109.2*(X-.5300)*(X-.4500)2如将三次样条函数加以整理，可用如下程序：s=collect(s);则输出结果为： s= .4595000000000-13.200*X3+9.9000000*X2-1.48800000000*X .37459306800000-4.2924*X3+3.66332200*X2-.137976090000*X .5180681700000-3.3917*X3+3.90576600*X2-.733130370000*X -.408614400000e-1+2.5768*X3-4.03188800*X2+2.868770320000*X例2 ： 给出节点的数据如下：-1.00 -0.54 0.13 1.12 1.89 2.06 2.54 2.82 3.50-2.46 -5.26 -1.87 0.05 1.65 2.69 4.56 7.89 10.31分别求在下列条件下在插值点，处的压紧三次样条插值，并显示该样条函数的有关信息：（1）端点约束条件为，；（2）端点约束条件为，.解 （1）输入MATLAB程序 X=-1.00 -0.54 0.13 1.12 1.89 2.06 2.54 2.82 3.50;Y=-2.46 -5.26 -1.87 0.05 1.65 2.69 4.56 7.89 10.31; XI=-0.02 2.56; YI= spline (X, 5,Y,29.16,XI), PP = spline (X, 5,Y,29.16)运行后屏幕显示压紧样条分别在，处的插值和该样条函数的有关信息如下YI = -3.1058 4.7834PP = form: pp breaks: -1 -0.5400 0.1300 1.1200 1.8900 2.0600 2.5400 2.8200 3.5000 coefs: 8x4 double pieces: 8 order: 4 dim: 1（2）因为端点约束条件为，所以输入MATLAB程序 YI= spline (X, 0,Y,0,XI), PP= spline (X, 0,Y,0)运行后屏幕显示压紧三次样条分别在，的插值和该样条函数的有关信息如下YI = -3.0192 4.7501PP = form: pp breaks: -1 -0.5400 0.1300 1.1200 1.8900 2.0600 2.5400 2.8200 3.5000 coefs: 8x4 double pieces: 8 order: 4 dim: 1例3：求有关分段三次样条图形的MATLAB主程序（一）限定端点约束条件的作图程序function S=splinetx(x0,y0,x,x,y,dy1,dyn)S = spline(x0,dy1,y0,dyn,x);Sn = spline(x0,dy1,y0,dyn,x); plot(x0,y0,o,x,Sn,-,x,S,*,x,y,-.)legend(节点(xi,yi), 分段三次样条函数,插值点(x,S),被插值函数y)（二）不限定端点约束条件的作图程序function S=splinetx1(x0,y0,xi,x,y)S= interp1(x0,y0,xi, spline); Sn= interp1(x0,y0,x, spline); plot(x0,y0,o,x,Sn,-,xi,S,*,x,y,-.)legend(节点(xi,yi), 分段三次样条函数,插值点(x,S),被插值函数y)（三）自由作图程序直接在MATLAB工作窗口编程序，例如，subplot(2,2,1),x1=-8:4/3:-4,c1=sin(x1);xx1 = -8:0.1:-4;pp1 = interp1 (x1,c1,xx1,spline ); cc1 =sin(xx1);%pp1 = spline (x1,c1,xx1); plot(x1,c1,bo,xx1,pp1,k-,xx1,cc1,r-.)subplot(2,2,2)x2=-4:4/3:-0;c2=sin(x2); xx2 = -4:0.1:0; pp2 = interp1 (x2,c2,xx2,spline ); cc2=sin(xx2);plot(x2,c2,bo,xx2,pp2,k-,xx2,cc2,r-.)title(y=sinx及其三次样条插值函数,节点(xi,yi)的图形)subplot(2,1,2)x=-8:4/3:8;c=sin(x);xx = -8:0.1:8;pp = spline(x,c,xx);cc=sin(xx); plot(x,c,bo,xx,pp,k-,xx,cc,r-.)legend(节点(xi,yi),三次样条插值函数,y=sinx 的函数),y=cosx 的函数)例4：（机床加工） 待加工零件的外形根据工艺要求由一组数据（x, y）给出（在平面情况下），用程控铣床加工时每一刀只能沿x方向和y方向走非常小的一步，这就需要从已知数据得到加工所要求的步长很小的（x, y）坐标.表 615给出的x,y数据位于机翼断面的下轮廓线上（如图625），假设需要得到x坐标每改变0.1时的y坐标.试完成加工所需数据，画出曲线，并求出x=0处的曲线斜率和13 x 15范围内y的最小值.表 615 机翼断面下轮廓线上的部分数据X035791112131415Y01.21.72.02.12.01.81.21.01.6= = = =图625 机翼断面轮廓线（表 615数据用圆点表示）解 根据上述提出的加工要求，以所给数据为节点，在x=0到x=15范围内求步长为0.1的插值.用四种插值方法试验，编写并保存名为sancili6710.m程序为M文件x0=0 3 5 7 9 11 12 13 14 15 ; x=0:0.1:15;y0=0 1.2 1.7 2.0 2.1 2.0 1.8 1.2 1.0 1.6 ; yL=lagr1(x0,y0,x);yX=interp1(x0,y0,x); yS=interp1(x0,y0,x,spline); yH=interp1(x0,y0,x, pchip); CZ=x yL yX yS yHsubplot(4,1,1)plot(x0,y0,bo,x,yL,r), grid,title(拉格朗日插值)subplot(4,1,2)plot(x0,y0,bo,x,yX,r), grid,title(分段线性插值)subplot(4,1,3)plot(x0,y0,bo,x,yS,r), grid,title(三次样条)subplot(4,1,4)plot(x0,y0,bo,x,yH,r), grid,title(分段埃尔米特插值)在MATLAB工作窗口输入文件名 sancili6710运行后得到的拉格朗日插值、分段线性插值、三次样条插值和分段埃尔米特插值及其节点的图形，同时还得到拉格朗日插值、分段埃尔米特插值、分段线性插值和三次样条插值的结果五 结论 MATLAB环境下编写求解三次样条插值的通用程序，可直接显示各区间段三次样条函数的具体表达式，计算出已给点的插值，最后显示各区间分段曲线图，为三次样条插值函数的应用提供简便方法。随着计算机技术及硬件设施的不断发展，计算机语言的演化从最开始的机器语言到汇编语言，最后到支持面向对象技术的面向对象语言。这就要求提供的计算方法也不断发展。同时在实践也给旧的插值逼近方法不断提出问题，这就要求新的插值逼近方法要更复杂，误差精度更高，同时解决更多方面的问题。插值逼近方法的发展也需要新的理论指导。自然辩证法的科学理论中提到科学“范式”概括了插值逼近法得法的发展过程辩证唯物主义自然观、自然科学发展过程及其规律，分析与综合、归纳法与演绎法、想象和类比等科学逻辑思维方法的应用都在插值逼近理论的发展过程中起到了重要作用。用科学的逻辑思维方法认识事物才会清楚的了解其过去、现在和未来计算数学中的插值逼近方法发展同样遵循着科学技术、科学理论发展的一般规律以自然辩证法的观点来分析有助于我们更加深入地认识插值逼近以及整个数值计算方法发展的历史、现状和趋势。插值逼近方法及数学理论的进一步发展也必将为自然辩证法的发展提供基础。本次的课程设计的整个过程让我认识到了基础知识的欠缺，通过查阅大量的资料，从根源上了解，三样条插值的由来和应用范围等，是我受益良多。切切实实的认识的努力学习的重要性。参考资料1 黄明游,冯果忱.数值分析M.北京:高等教育出版社,2008. 2 马东升,雷勇军.数值计算方法M.北京:机械工业出版社,2006. 3 石博强,赵金.MATLAB数学计算与工程分析范例教程M.北京:中国铁道出版社.2005.4郝红伟，MATLAB 6，北京，中国电力出版社，20015姜健飞，胡良剑，数值分析及其MATLAB实验， 科学出版社，20046薛毅，数值分析实验，北京工业大学出版社，2005翻译Switch语句在几个case表达式基础上的转换。语法：case case_expr statement,.,statement case case_expr1,case_expr2,case_expr3,. statement,.,statement. otherwise statement,.,statementend用法讨论：switch表达式的语法足一种有条件的执行码。特别地，switch执行一系列任意的可供选强的数字下的语句，每一个可选项组成一个case，由以下3个部分组成：(1) case语句。(2) 一个或多个表达式。(3) 一个或多个语句。在基本语法中，switch执行的语句必须满足switch_exprcase_expr。当case表达式为单元数组时(如上面的第个case所示)，case_expr只有在单元数组元中有元素匹配switch表达式时才匹配上。如果没有case表达式匹配switch表达式，则控制转换到otherwise case(如果该项存在的话)。Case执行完后，程序必须从switch末尾开始继续执行。switch_expr是一个标量或者一个字符串。如果switch_exprcase_expr成立的，标量switch_expr匹配上case_expr。当且仅当strcmp(switch_expr,case_expr)=1(真)时 ，一个字符串switch_expr与case_expr匹配。应用实例：执行的列编码块基于what the string, method, is set to, methodmethod = Bilinear;switch lower(method) case linear,bilinear disp(Method is linear) case cubic disp(Method is cubic) case nearest disp(Method is nearest) otherwise disp(Unknown method.)endMethod is linear对照英文：Switch among several cases based on expression Syntaxswitch switch_expr case case_expr statement,.,statement case case_expr1,case_expr2,case_expr3,. statement,.,statement. otherwise statement,.,statementendDiscussionThe switch statement syntax is a means of conditionally executing code. In particular, switch executes one set of statements selected from an arbitrary number of alternatives. Each alternative is called a case, and consists of The case statement One or more case expressions One or more statements In its basic syntax, switch executes the statements associated with the first case whe