【信息与计算科学专业毕业论文】三次样插值在Delphi与MATLAB中的实现.doc

毕 业 论 文学生姓名： 学 号： 学 院： 理学院 专 业： 信息与计算科学 题 目： 三次样插值在Delphi与MATLAB 中的实现 指导教师： 评阅教师： * 年 *月河北科技大学毕业设计（论文）成绩评定表姓 名李江娜学 号05104131成 绩专 业信息与计算科学题 目三次样插值在Delphi与MATLAB中的实现 指导教师评语及成绩 指导教师： 年 月 日评阅教师评语及成绩 评阅教师： 年 月 日答辩小组评语及成绩答辩小组组长： 年 月 日答辩委员会意见 学院答辩委员会主任： 年 月 日 注：该表一式两份，一份归档，一份装入学生毕业设计说明书（论文）中。毕 业 设 计（论文） 中 文 摘 要随着科学技术的发展，科学与工程计算已被推向科学活动的前沿，科学与工程计算的范围扩大到了所有科学领域，并与实验、理论三足鼎立，相辅相成，成为人类科学活动的三大方法之一。因此熟练的运用计算机进行科学计算，已成为科技工作者的一项基本技能。这就要求人们去研究和掌握适用于计算机上使用的数值计算方法。三次样条插值是科学与工程计算中的一类主要算法，本文主要是以数值分析中三次样条插值算法为基础，借助于Matlab数学软件进行处理，以Delphi开发软件为实现平台，将Delphi开发软件和Matlab数学工具相结合,利用Delphi灵活强大方便的编程功能,以及Matlab强大的科学计算功能来实现三次样条插值的可视化计算，从而使计算更为简单、实用且易操作。关键词 三次样条 插值函数 Matlab程序 Dephi程序毕 业 设 计（论文） 外 文 摘 要Title Realization of Cubic Spline Interpolation Fuction in Delphi and MATLAB AbstractAlong with the development of technology，science and engineering calculation have already occupied the forward position of scientific activities ，and its range has spreaded to all the scientific fields.It has the same status with experiment and theory, and supplements each other，becomes one of the three famous methods of human scientificActivities. So it is essential skill for engineer to master calculation,and it is necessary for people to research calculation and master calculation in computer Cubic spline is a very important arithmetic in the science and engineering calculation .This article is based on the cubic spline arithmetic of numerical analysis ,uses Matlab software to deal with the data, and Delphi as the actualizing flat, combine Delphi with Matlab，then make full use of the smart, strong, convenient programming function of Delphi and the powerful calculating capacity of Matlab, and carrys out the visible calculation of cubic spline, which is easier and more practical.Key Words Cubic spline Interpolation function Matlab program Delphi program目 录1 引言12 三次样条函数12.1 定义12.2 边界条件22.3三次样条函数插值的构造32.4 三次样条插值的算法63 用MATLAB做三次样条插值计算73.1 三次样条插值的源代码73.2 三次样条插值的各种边界条件在MATLAB中的命令124 程序设计164.1 设计步骤164.2程序代码18结 论22致 谢23参 考 文 献241 引言我们知道，高次插值不仅计算复杂，而且可能出现龙格现象。而分段低次插值有计算简单、插值曲线连续及能逼近函数的特点，但在分段插值多项式的分段点处有“尖点”出现，及光滑性比较差。例如，分段线性插值是在给定了插值节点上的函数值和微商值以后，构造一个整体上具有一阶连续微商的插值函数。为了提高整体光滑度，并且在只给出插值点上的函数值的情况下构造一个整体上具有二阶连续微商的插值函数，这就是三次样条插值。样条这一名词来源于工程中的样条曲线。绘图员为了将一些指定的点（或样点）连接成一条光滑曲线，往往用细长的木条（或样条）把相近的几点连接在一起，在逐步延伸连接起全部节点，使其形成一条光滑的样条曲线，它在连接处具有连续曲率。对绘图员的样条曲线进行数学模拟，得到的函数叫做样条函数，它在连接处具有一阶和二阶连续微商。本论文以数值分析中三次样条插值算法和Matlab数学软件为基础，以Delphi为实现平台来实现三次样条插值的可视化计算，从而使计算更为简单、更方便。Delphi作为一种功能强大的编程工具,具有易学易用、开发效率高,界面制作美观方便等优点，因此被很多程序员所青睐。Pascal作为历史上第一种结构化的高级语言，在从事复杂算法编写方面也有着诸多优点，可是在软件开发快速运作的今天，用Pascal原始开发一些复杂的算法，不仅编译效率不高而且也影响开发进度，因此Delphi如何与科学计算软件相结合，从而高效地完成编程任务成为一个困扰很多程序员的问题。而Matlab就正是一种高效率的科学工程计算语言,它在矩阵运算、数值计算、数字信号处理、系统识别、自动控制、神经网络、图形显示等方面比其它语言有难以比拟的优势。将Delphi和Matlab相结合,利用Delphi灵活强大方便的编程能力,Matlab强大的科学计算能力就可以开发出功能强大、操作灵活的软件。 2 三次样条函数 样条插值属于分段光滑插值。它的基本思想是，在由两个相邻结点所构成的每一个小区间内用低次多项式来逼近，并且，在各结点的连接处又保证是光滑的（即导数连续）。2.1 定义定义：设函数定义在区间上，给定个节点： 若满足（1）在上有二阶连续的导数（2）在每个小区间上是一个三次多项式则称为三次样条函数。若还满足 (3) 则称为三次样条插值函数。2.2 边界条件由样条插值函数的定义可知，在每个小区间上都是三次函数，即 （2.1） 其中系数、待定，它满足条件：（1） 插值条件 （2.2）（2） 连续条件 （2.3）式（2.2）、式（2.3）共给出了个条件，而待定系数有个，因此还需要个条件才能确定。通常是在区间端点、上个加一个条件，称为边界条件。（3）边界条件常用的边界条件有三种类型：.已知两端点的一阶导数值，即令 ， （2.4）.已知两端点的二阶导数值，即令 ， （2.5）.当是以为周期的函数时，则要求也是周期函数，这时边界条件应满足，当时， 。2.3三次样条函数插值的构造从原则上来讲，我们可以利用方程组（2.2）、（2.3）以及边界条件（2.4）或边界条件（2.5）、（2.6），求出个待定系数、，从而得到在每个小区间上的表达式。但是，由于这种做法的计算量相当大，在实际计算时很少用到。这里介绍一种简单的求法，因为三次样条插值函数可以有多种表达方法，既可以用节点处的二阶导数表示，也可以用节点处的一阶导数表示。所以为了简化计算，我们利用在节点处的二阶导数值表示；用在内节点上的连续性和边界条件来确定。具体步骤如下：1.的表达式由表示因为在每个小区间上是三次多项式，所以在此小区间上是一次多项式，并设，则的表达式为 （2.6）将上式积分两次并利用条件和，可得 （2.7）2.利用连续条件列出关于的方程组由式(2.7)可知，只要求出，就可以确定分段三次样条插值函数在每个小区间上的表达式。为求，将式(2.7)求导，得到 (2.8)这样 (2.9) (2.10)由可得 整理得 (2.11)其中 , 3.利用边界条件来求解对于型边界条件，把、分别代入式（2.9）和式（2.10），得 即 （2.12） （2.13）其中，。令、，将式（2.11）、式（2.12）、式（2.13）三式合成为矩阵形式： （2.14）对于型边界条件，把、代入，由式（2.11）得 相应的矩阵方程为 （2.15）对于型边界条件，因为，由（2.9）式、式（2.10）得 由得，代入上式整理得 其中 ， 相应的矩阵方程为 （2.16）方程组（2.14）-（2.16）对应的系数矩阵矩阵均为严格对角占优阵，因此都有唯一解，这样求这样插值函数的步骤为（1） 求、 ；（2） 对不同边界条件选用相应的矩阵方程，求；（3） 代入式（2.7），求出三次样条插值函数。2.4 三次样条插值的算法对于定义在点上的函数，构造三次样条插值，满足边界条件和。输入 数函数值；边界条件、及待求点输出 函数值Step1 对于，计算 Step2 对于，计算 Step3 计算 Step4 求解线性方程祖，解出： Step5 当时执行Step6；Step6 判定所在的区间；Step7 按照上述公式计算： Step8 输出。3 用MATLAB做三次样条插值计算3.1 三次样条插值的源代码在Matlab 环境下根据上述算法步骤进行编程，源程序如下：function output = spline(x,y,xx) %三次样条插值%yi=spline(x,y,xi)等价于YI=interp(x,y,xi,'spline'),% 根据数据(x,y)给出在xi的线性插值结果yi.% 使用"非扭结"端点条件, 即强迫第一二端多项式三次项系数相同,% 最后一段和倒数第二段三次项系数相同.%pp=spline(x,y)返回样条插值的分段多项式(pp形式).%breaks,coefs=unmkpp(pp)将pp形式展开,其中breaks为结点,coefs为各段多项式系数.%yi=ppval(pp,xi),pp形式在xi的函数值.%例 考虑数据% x | 1 2 4 5% -|-% y | 1 3 4 2 % clear;close;% x=1 2 4 5;y=1 3 4 2;% p=spline(x,y);% xi=1:0.1:5;yi=ppval(p,xi);% plot(x,y,'o',xi,yi,'k');% title('not-a-knot SPLINE');% b,c=unmkpp(p)%另一个例子见下列英文部分%SPLINE Cubic spline data interpolation.% YY = SPLINE(X,Y,XX) uses cubic spline interpolation% to find a vector YY corresponding to XX. X and Y are the% given data vectors and XX is the new abscissa vector.% The ordinates Y may be vector-valued, in which case Y(:,j) is% the j-th ordinate.% PP = SPLINE(X,Y) returns the pp-form of the cubic spline % interpolant instead, for later use with PPVAL, etc.% Ordinarily, the not-a-knot end conditions are used. However, if% Y contains two more ordinates than X has entries, then the first% and last ordinate in Y are used as the endslopes for the cubic spline.% Here's an example that generates a coarse sine curve, then% samples the spline over a finer mesh:% x = 0:10; y = sin(x);% xx = 0:.25:10;% yy = spline(x,y,xx);% plot(x,y,'o',xx,yy)% Here is an example that features a vector-valued spline, along with complete% spline interpolation, i.e., fitting to given end slopes (instead of using the% not-a-knot end condition); it uses SPLINE to generate a circle:% circle = spline( 0:4, 0 1 0 -1 0 1 0; pi/2 0 1 0 -1 0 pi/2 );% xx = 0:.1:4; cc = ppval(circle, xx); plot(cc(1,:), cc(2,:), axis equal% See also INTERP1, PPVAL, SPLINES (The Spline Toolbox).% Carl de Boor 7-2-86% Revised 11-24-87 JNL, 6-16-92 CBM, 10-14-97 CB.% Copyright (c) 1984-98 by The MathWorks, Inc.% $Revision: 5.11 $ $Date: 1997/12/03 19:22:33 $% Generate the cubic spline interpolant in pp form, depending on the% number of data points (and usually using the not-a-knot end condition).output=;n=length(x);if n<2, error('There should be at least two data points.'), endif any(diff(x)<0), x,ind=sort(x); else, ind=1:n; endx=x(:); dx = diff(x);if all(dx)=0, error('The data abscissae should be distinct.'), endyd,yn = size(y); % if Y happens to be a column matrix, change it to % the expected row matrix.if yn=1, yn=yd; y=reshape(y,1,yn); yd=1; endif yn=n notaknot = 1;elseif yn=n+2 notaknot = 0; endslopes = y(:,1 n+2).' y(:,1 n+2)=;else error('Abscissa and ordinate vector should be of the same length.')endyi=y(:,ind).' dd = ones(1,yd);dx = diff(x); divdif = diff(yi)./dx(:,dd);if n=2 if notaknot, % the interpolant is a straight line pp=mkpp(x.',divdif.' yi(1,:).',yd); else % the interpolant is the cubic Hermite polynomial divdif2 = diff(endslopes(1,:);divdif;endslopes(2,:)./dx(1 1,dd); pp = mkpp(x,. (diff(divdif2)./dx(1,dd).' (2 -1*divdif2).' . endslopes(1,:).' yi(1,:).',yd); endelseif n=3&notaknot, % the interpolant is a parabola yi(2:3,:)=divdif; yi(3,:)=diff(divdif)/(x(3)-x(1); yi(2,:)=yi(2,:)-yi(3,:)*dx(1); pp = mkpp(x(1),x(3),yi(3 2 1,:).',yd);else % set up the sparse, tridiagonal, linear system for the slopes at X . b=zeros(n,yd); b(2:n-1,:)=3*(dx(2:n-1,dd).*divdif(1:n-2,:)+dx(1:n-2,dd).*divdif(2:n-1,:); if notaknot x31=x(3)-x(1);xn=x(n)-x(n-2); b(1,:)=(dx(1)+2*x31)*dx(2)*divdif(1,:)+dx(1)2*divdif(2,:)/x31; b(n,:)=. (dx(n-1)2*divdif(n-2,:)+(2*xn+dx(n-1)*dx(n-2)*divdif(n-1,:)/xn; else x31 = 0; xn = 0; b(1 n,:) = dx(1 n-2,dd).*endslopes; end c = spdiags( dx(2:n-1);xn;0 . dx(2);2*dx(2:n-1)+dx(1:n-2);dx(n-2) . 0;x31;dx(1:n-2) ,-1 0 1,n,n); % sparse linear equation solution for the slopes mmdflag = spparms('autommd'); spparms('autommd',0); s=cb; spparms('autommd',mmdflag); % convert to pp form c4=(s(1:n-1,:)+s(2:n,:)-2*divdif(1:n-1,:)./dx(:,dd); c3=(divdif(1:n-1,:)-s(1:n-1,:)./dx(:,dd) - c4; pp=mkpp(x.',. reshape(c4./dx(:,dd).' c3.' s(1:n-1,:).' yi(1:n-1,:).',(n-1)*yd,4),yd);endif nargin=2 output=pp;else output=ppval(pp,xx);end设个结点以数组输入（程序中用个结点），个插值点以数组输入，输出数组为个插值。3.2 三次样条插值的各种边界条件在MATLAB中的命令MATLB中有现成的命令：Y1=spline(x,y,xi)等价于Y2=spline(x,y,xi,spline)其中Pp=spline(x,y)返回样条插值的多项式（pp）形式。Pp=cspe(x,y,边界类型,边界值)生成各种边界条件的三次样条插值。其中，边界类型可为not-a-kont，表示非扭条件，不用给定边界值（默认）；complete，表示给定边界的一阶导数；second，给定边界二阶导数；periodic，周期性边界条件，不用给定边界值；variational，自然样条（边界二阶导数为0）；yi=ppval(pp,xi) pp 样条在xi的函数；fnplt(pp)画出pp样条的图。3.2.1 非纽条件“非扭结条件”指分段三次样条插值的第一段、第二段多项式的三次项系数相同，最后一段和倒数第二段的三次项系数相同。>> x=0.1 0.2 0.15 0 -0.2 0.3;>> y=0.95 0.84 0.86 1.06 1.50 0.72;>> pp=spline(x,y)pp = form: 'pp' breaks: -0.2000 0 0.1000 0.1500 0.2000 0.3000 coefs: 5x4 double pieces: 5 order: 4 dim: 1>> pp.coefsans = -36.3850 21.8592 -5.1164 1.5000 -36.3850 0.0282 -0.7390 1.0600 227.6995 -10.8873 -1.8249 0.9500 -143.0047 23.2676 -1.2059 0.8600 -143.0047 1.8169 0.0484 0.8400显示样条函数的5个分段三次多项式；Spline使用“非扭结“端点条件，即强迫第一段、第二段多项式的三次项系数相同，最后一段和倒数第二段的三次项系数相同。>> x=0.1 0.2 0.15 0 -0.2 0.3;>> y=0.95 0.84 0.86 1.06 1.50 0.72;>> pp=csape(x,y); %等价于spline(x,y)>> xi=-0.2:0.01:0.3;>> yi=ppval(pp,xi);>> fnplt(pp) %画出样条的图，图如下：3.2.2 给定边界一阶导数若边界条件，则有，>> pp2=csape(x,y,'complete',0.2,1.5);pp2.coefsans = 73.7278 -26.7456 0.2000 1.5000 -119.8225 17.4911 -1.6509 1.0600 348.0473 -18.4556 -1.7473 0.9500 -442.0118 33.7515 -0.9825 0.8600 297.7515 -32.5503 -0.9225 0.84003.2.3 给定边界二阶导数若边界条件，则有，>> pp2=csape(x,y,'second',1.0,0.5);pp2.coefs ans = 11.9962 0.5000 -2.7798 1.5000 -72.9468 7.6977 -1.1403 1.0600 279.3923 -14.1863 -1.7892 0.9500 -269.5085 27.7225 -1.1124 0.8600 43.1792 -12.7038 -0.3614 0.84003.2.4 周期性边界条件>> pp2=csape(x,y,'periodic');pp2.coefsans = 28.1818 -6.6867 -1.9899 1.5000 -83.9448 10.2224 -1.2828 1.0600 281.8831 -14.9610 -1.7567 0.9500 -250.9740 27.3214 -1.1386 0.8600 12.1266 -10.3247 -0.2888 0.84003.2.5 自然样条>> pp2=csape(x,y,'variational');pp2.coefsans = 13.1233 0.0000 -2.7249 1.5000 -73.7265 7.8740 -1.1501 1.0600 279.7319 -14.2440 -1.7871 0.9500 -268.9008 27.7158 -1.1135 0.8600 42.0643 -12.6193 -0.3587 0.84004 程序设计4.1 设计步骤 步骤1：建立名为draw数据库； 字段名称 数据类型含义idtopicresultjpg自动编号文本备注OLE对象识别号名称插值结果插值图像 步骤2：打开Dephi，新建工程，保存工程。步骤3：在Dephi中建立界面，根据需要放入各种控件，并修改其属性，界面（程序运行后的界面）如下：.输入插值点和边界条件及类型；.显示出插值图像及分段函数的系数；.记录插值函数及图像；步骤4：建立数据库的链接，选中“ADOConnection1”组件，在对象查看器里，选择“ConnectionString”属性，在单击右面的图标，显示出“Form1->ADOConnection1 ConnectionString”对话框里，单击“Build”按钮，打开“数据链接热属性”对话框。在该对话框里，在“提供程序”页上选择连接的数据为“Microsoft Jet 4.0 OLE DB Provider”,然后单击“下一步”按钮，在“连接”页里，在“指定下列设置以连接Access数据”选项中的“选择或输入数据名称”文本框右边的图标上单击，在弹出的对话框里，选择要连接的Access数据库，完成后，单击“测试连接”按钮进行测试。若弹出“Microsoft数据连接”对话框并提示连接成功，则表明连接成功。若出现错误，则表明数据链接失败，这时就要对数据路径等进行检查。步骤5：选中“ADOQuery1”组件，在对象查看器里，选择“Connetion”属性，在单击右面的图标，与“ADOConnection1”关联；选中“DataSource1”组件，选作择“DataSource”与“ADOQuery1”关联。步骤6：将 “DBGrid1”, “DBImage1” ,“DBMemo1”与“DataSource”相关联；步骤7：规划菜单：设总菜单有“编辑”和“退出”；总菜单下“编辑”下有“清空”，用于清空数据库中的数据，“退出”用于退出程序。添加主菜单控件：双击“Standard”控件页中的“TMainMenu”控件，这样在“Form1”表单上出现名为“MainMenu1”控件，这个控件为不可见控件。双击“MainMenu1”控件，在对象查看器的“Caption”属性中，填写菜单的名称，按“Enter”键即可。4.2程序代码在implementation下填写 uses ComObj ;在“Button1”下编写代码：procedure TForm1.Button1Click(Sender: TObject);varMatlab:OleVariant;ReturnString:String;x,y,z,z1,pp,pp2:String;strList:TStrings;str:String;pic:TPicture;begin pic:=TPicture.Create(); x:='x='+Edit_X.Text+'' y:='y='+Edit_Y.Text+'' z:='z='+Edit_z.Text;if ComboBox1.ItemIndex=0 thenbegin pp:='pp=csape(x,y);' pp2:='pp2=pp.coefs' z1:='fnplt(pp)' ; endelse if ComboBox1.ItemIndex=1 thenbegin pp:='pp=csape(x,y,'+'''complete'''+','+Edit_z.Text+');' pp2:='pp2=pp.coefs' z1:='fnplt(pp)' ;endelse if ComboBox1.ItemIndex=2 thenbegin pp:='pp=csape(x,y,'''+'second'''+','+Edit_z.Text+');' pp2:='pp2=pp.coefs' z1:='fnplt(pp)' ; e