《实验四数据插值与拟合幻灯片.ppt》由会员分享，可在线阅读，更多相关《实验四数据插值与拟合幻灯片.ppt（49页珍藏版）》请在taowenge.com淘文阁网|工程机械CAD图纸|机械工程制图|CAD装配图下载|SolidWorks_CaTia_CAD_UG_PROE_设计图分享下载上搜索。
1、实验四数据插值与实验四数据插值与拟合拟合第1页，共49页，编辑于2022年，星期五4.1 MATLAB中的插值函数中的插值函数n n函数插值来源于函数的以下问题：只知道函数在某区间有定义且已得到区间内一些离散点的值，希望用简单的表达式近似给出函数在此区间上的整体描述，并能与已知离散点上的值相等。n n插值法按插值函数的形式主要分为以下几种形式：（1）代数多项式插值；（2）三角多项式插值；（3）有理分式插值。第2页，共49页，编辑于2022年，星期五n n代数多项式插值是最常用的插值方式，其内容也是最丰代数多项式插值是最常用的插值方式，其内容也是最丰富的，它又可分为以下几种插值方式：富的，它又可
2、分为以下几种插值方式：（1 1）非等距节点插值非等距节点插值，包括拉格朗日插值、利用均差的牛，包括拉格朗日插值、利用均差的牛顿插值和埃特金插值；顿插值和埃特金插值；（2 2）非等距节点插值非等距节点插值，包括利用差分的牛顿插值和高斯，包括利用差分的牛顿插值和高斯插值等；插值等；（3 3）在插值中增加了导数的）在插值中增加了导数的HermiteHermite（埃尔米特）插值（埃尔米特）插值；（4 4）分段插值，包括分段线性插值、分段，包括分段线性插值、分段HermiteHermite（埃（埃尔米特）插值和样条函数插值；尔米特）插值和样条函数插值；（5）反插值反插值。n n按被插值函数的变量个数还
3、可把插值法分为按被插值函数的变量个数还可把插值法分为一元插值一元插值和和多多元插值元插值。第3页，共49页，编辑于2022年，星期五4.1.1 一元插值函数一元插值函数n nMATLABMATLAB中的一元插值函数为中的一元插值函数为interp1()interp1()，它的功能是一维数据插，它的功能是一维数据插值（表格查找）。该命令对数据点之间进行计算内插值，它出值（表格查找）。该命令对数据点之间进行计算内插值，它出一元函数一元函数f(x)f(x)在中间点的数值，其中函数在中间点的数值，其中函数f(x)f(x)由所给数据决定。由所给数据决定。n n一元插值函数一元插值函数interp1()i
4、nterp1()的几种调用格式如表的几种调用格式如表4-14-1所示。所示。表表4-1 4-1 一维插值插值函数一维插值插值函数interp1interp1的语法格式的语法格式语法形式语法形式说明说明 y=interp1(y=interp1(x x,Y Y,x,xi i)由已知点集由已知点集(x x,Y Y）插值计算）插值计算x xi i上的函数值上的函数值y=interp1(y=interp1(x x,Y Y,x,xi i)相当于相当于x x=1:length(=1:length(Y Y)的的interp(interp(x x,Y Y,x,xi i)y=interp1(y=interp1(x
5、 x,Y Y,x,xi i,method),method)用指定插值方法计算插值点用指定插值方法计算插值点x xi i上的函数值上的函数值y=interp1(y=interp1(x x,Y Y,x,xi i,method,method,extrapextrap)对对x xi i中超出已知点集的插值点用指定插值中超出已知点集的插值点用指定插值方法计算函数值方法计算函数值y=interp1(y=interp1(x x,Y Y,x,xi i,method,method,extrapextrap,extrapval),extrapval)用指定方法插值用指定方法插值x xi i上的函数值，超出已知上的
6、函数值，超出已知点集处函数值取点集处函数值取extrapvalextrapvaly=interp1(y=interp1(x x,Y Y,x,xi i,method,method,pppp)用指定方法插值，但返回结果为分段多用指定方法插值，但返回结果为分段多项式项式第4页，共49页，编辑于2022年，星期五method方法描述方法描述 nearestnearest 最邻近插值：插值点处函数值取与插值点最邻近的已知点的函数值最邻近插值：插值点处函数值取与插值点最邻近的已知点的函数值 linerliner 分段线性插值：插值点处函数值由连接其最邻近的两侧点的线性函分段线性插值：插值点处函数值由连接其
7、最邻近的两侧点的线性函数预测，数预测，MATLABMATLAB中中interp1interp1的默认方法的默认方法 splinespline 样条插值：默认为三次样条插值。可用样条插值：默认为三次样条插值。可用splinespline函数代替函数代替 pchippchip 三次三次HermiteHermite多项式插值。可用多项式插值。可用pchippchip函数代替函数代替 cubiccubic 同同 pchippchip，三次，三次HermiteHermite多项式插值多项式插值n nMATLABMATLAB中一维插值有多种算法，由interp1函数中的methodmethod指定。指定。
8、MATLABMATLAB中一维插值的各种算法如表中一维插值的各种算法如表4-24-2所示。所示。表表4-2 4-2 一维插值算法一维插值算法(method)第5页，共49页，编辑于2022年，星期五1.Linear（分段线性插值）（分段线性插值）n n它的算法是在每个小区间它的算法是在每个小区间x xi i,x,xi+1i+1 上采用简单的线性插值。在上采用简单的线性插值。在区间区间x xi i,x,xi+1i+1 上的子插值多项式为：上的子插值多项式为：由此整个区间由此整个区间x xi i,x,xi+1i+1 上的插值函数为：上的插值函数为：其中其中 定义如下：定义如下：第6页，共49页，编
9、辑于2022年，星期五n n分段线性插值方法在速度和误差之间取得了比较好的均衡，其插值函数具有连续性，但在已知数据点处的斜率一般不会改变，因此不是光滑的。分段线性插值方法是MATLAB一维插值默认的方法。第7页，共49页，编辑于2022年，星期五2.Spline（样条插值）（样条插值）n n样条插值是用分段低次多项式去逼近函数。样条函数可以给出光滑 的插值曲线，只要在插值区间端点提供某些导数信息，样条插值可以适应不同光滑需求。三次样条是使用最为广泛的样条插值，它在每个子区间xi,xi+1上都是有二阶连续导数的三次多项式，即 其中 都是三次多项式。第8页，共49页，编辑于2022年，星期五n n
10、对于给定的离散的测量数据经x，y（称为断点），要寻找一个三次多项式y=p(x)，以逼近每对数据(xi,yi)点间曲线。过两点(xi,yi)和(xi+1,yi+1)只能确定一条直线，而通过一点的三次多项式曲线有无穷多条。为使通过中间断点的三次多项式曲线具有唯一性，要增加以下的连续条件和边界条件（因为三次多项式有4个系数）：（1）三次多项式在点(xi,yi)处有：；（2）三次多项式在点(xi,yi)处有：；（3）三次多项式在点(xi,yi)处有：；（4）边界条件：。第9页，共49页，编辑于2022年，星期五n n表表4-24-2中各种方法中：中各种方法中：（1 1）nearestnearest方法
11、速度最快，占用内存最小，但一般来说误差最大，插值结果最不光滑；（2 2）spline三次样条插值是所有插值方法中运行耗时最长的，其插值函数以及插值函数的一阶、二阶导函数都连续，因此是最光滑的插值方法，占用内存上比cubic方法小，但当已知数据点不均匀分布时可能出现异常结果。（3 3）cubiccubic三次多项式插值法中插值函数及其一阶导数都三次多项式插值法中插值函数及其一阶导数都是连续的，因此其插值结果也比较光滑，运算速度比是连续的，因此其插值结果也比较光滑，运算速度比splinespline方法略快，但占用内存最多。在实际的使用中，应根方法略快，但占用内存最多。在实际的使用中，应根据实际需
12、求和运算条件选择合适的算法。据实际需求和运算条件选择合适的算法。第10页，共49页，编辑于2022年，星期五例例4-1 用用interp1对对sin函数进行分段线性插值。函数进行分段线性插值。解：解：在在MATLABMATLAB命令窗口中输入以下命令：命令窗口中输入以下命令：x=0:2*pi;x=0:2*pi;y=sin(x);y=sin(x);xx=0:0.5:2*pi xx=0:0.5:2*pi yy=interp1(x,y,xx);yy=interp1(x,y,xx);plot(x,y,s,xx,yy)plot(x,y,s,xx,yy)注：注：例例4-14-1中用默认的中用默认的(分段线
13、性插值的分段线性插值的linear)linear)对已知的对已知的7 7个个sinsin函数的函数的数据点进行插值，用数据点进行插值，用plotplot画出插值结果。从图中可以看出分段线性就是联结两个画出插值结果。从图中可以看出分段线性就是联结两个邻近的已知点的线性函数插值计算该区间内插值点上的函数邻近的已知点的线性函数插值计算该区间内插值点上的函数值。值。第11页，共49页，编辑于2022年，星期五例例4-2 用其他一维插值方法对以下用其他一维插值方法对以下7个离散数据点个离散数据点(1,3.5)、(2,2.1)、(3,1.3)、(4.0.8)、(5,2.9)、(6,4.2)、(7,5.7)
14、进行一维插值方法。进行一维插值方法。解：解：在在MATLABMATLAB命令窗口中输入以下命令：命令窗口中输入以下命令：x=1 2 3 4 5 6 7;y=3.5 2.1 1.3 0.8 2.9 4.2 5.7;xx=1:0.5:7;y1=interp1(x,y,xx,nearest);y2=interp1(x,y,xx,spline);y3=interp1(x,y,xx,cubic);plot(x,y,o,xx,y1,-,xx,y2,-.,xx,y3,:)第12页，共49页，编辑于2022年，星期五第13页，共49页，编辑于2022年，星期五4.2 拉格朗日插值法拉格朗日插值法n n拉格朗日
15、插值法拉格朗日插值法是基于基函数的插值方法，插值多项式可表示是基于基函数的插值方法，插值多项式可表示为为 其中其中 称为称为i i次基函数：次基函数：第14页，共49页，编辑于2022年，星期五n n在MATLAB中编程实现拉格朗日插值法拉格朗日插值法函数为：LanguageLanguage。n n功能：求已知数据点的拉格朗日多项式；功能：求已知数据点的拉格朗日多项式；n n调用格式：f=Language(x,y)Language(x,y)或或f=f=Language(x,y,x0)Language(x,y,x0)。其中，其中，x为已知数据点的为已知数据点的x x 坐标向量；坐标向量；y为已知
16、数据点的y y 坐标向量；坐标向量；x0 x0为插值点的x x坐标；坐标；f f为求得的拉格朗日多项式或x0 x0处的插值。第15页，共49页，编辑于2022年，星期五 function f=Language(x,y,x0)%求已知数据点的拉格朗日多项式%已知数据点的x 坐标向量:x%已知数据点的y 坐标向量:y%为插值点的x坐标:x0%求得的拉格朗日多项式或x0处的插值:fsyms t;if(length(x)=length(y)n=length(x);else disp(x和y的维数不相等！);return;end%检错f=0.0;for(i=1:n)l=y(i);for(j=1:i-1)
17、l=l*(t-x(j)/(x(i)-x(j);end;for(j=i+1:n)l=l*(t-x(j)/(x(i)-x(j);%计算拉格朗日基函数 end;f=f+l;%计算拉格朗日插值函数 simplify(f);%化简 if(i=n)if(nargin=3)f=subs(f,t,x0);%计算插值点的函数值 else f=collect(f);%将插值多项式展开 f=vpa(f,6);%将插值多项式的系数化成6位精度的小数 end endend第16页，共49页，编辑于2022年，星期五例例4-3 根据下表的数据点求出其拉格朗日插根据下表的数据点求出其拉格朗日插值多项式，并计算当值多项式，并
18、计算当x=1.6时时y的值。的值。解：解：x=1 1.2 1.8 2.5 4;x=1 1.2 1.8 2.5 4;y=0.8415 0.9320 0.9738 0.5985-0.7568;y=0.8415 0.9320 0.9738 0.5985-0.7568;f=language(x,y)f=language(x,y)f=f=1.05427*t-.145485e-1*t2-.204917*t3+.328112e-1*t4-.261189e-1 1.05427*t-.145485e-1*t2-.204917*t3+.328112e-1*t4-.261189e-1 f=language(x,y,
19、1.6)f=language(x,y,1.6)f=f=0.9992 0.9992x x1 11.21.21.81.82.52.54 4 y y0.84150.84150.93200.93200.97380.97380.59850.5985-0.7568-0.7568第17页，共49页，编辑于2022年，星期五4.3 利用均差的牛顿插值法利用均差的牛顿插值法n n函数函数 f f的零阶均差定义为的零阶均差定义为 ，一阶定义均差为，一阶定义均差为 一般地，函数一般地，函数f f的的k k阶均差定义为：阶均差定义为：利用均差的牛顿插值法利用均差的牛顿插值法多项式为多项式为第18页，共49页，编辑于2
20、022年，星期五系数的计算过程如表系数的计算过程如表4-3所示所示表表4-3 4-3 均差计算表格均差计算表格一阶均差一阶均差二阶均差二阶均差三阶均差三阶均差n n阶均差阶均差第19页，共49页，编辑于2022年，星期五n n在MATLAB中编程实现利用均差牛顿插值法利用均差牛顿插值法函数为：NewtonNewton。n n功能：求已知数据点的均差形式的牛顿插值多项式；功能：求已知数据点的均差形式的牛顿插值多项式；n n调用格式：调用格式：f=f=Newton(x,y)或或f=Newton(x,y,x0)。其中，其中，x x为已知数据点的为已知数据点的x 坐标向量；坐标向量；y y为已知数据点
21、的为已知数据点的y y 坐标向量；坐标向量；x0为插值点的x x坐标；坐标；f f为求得的牛顿插值法多项式或x0 x0处的插值。处的插值。第20页，共49页，编辑于2022年，星期五function f=Newton(x,y,x0)function f=Newton(x,y,x0)%求已知数据点的均差形式牛顿插值多项式求已知数据点的均差形式牛顿插值多项式%已知数据点的已知数据点的x x 坐标向量坐标向量:x:x%已知数据点的已知数据点的y y 坐标向量坐标向量:y:y%为插值点的为插值点的x x坐标坐标:x0:x0%求得的均差形式牛顿插值多项式或求得的均差形式牛顿插值多项式或x0 x0处的插值
22、处的插值:f:fsyms t;syms t;if(length(x)=length(y)if(length(x)=length(y)n=length(x);n=length(x);c(1:n)=0.0;c(1:n)=0.0;elseelse disp(x disp(x和和y y的维数不相等！的维数不相等！););return;return;endend第21页，共49页，编辑于2022年，星期五f=y(1);f=y(1);y1=0;y1=0;l =1;l =1;for(i=1:n-1)for(i=1:n-1)for(j=i+1:n)for(j=i+1:n)y1(j)=(y(j)-y(i)/(x
23、(j)-x(i);y1(j)=(y(j)-y(i)/(x(j)-x(i);end end c(i)=y1(i+1);c(i)=y1(i+1);l=l*(t-x(i);l=l*(t-x(i);f=f+c(i)*l;f=f+c(i)*l;simplify(f);simplify(f);y=y1;y=y1;if(i=n-1)if(i=n-1)if(nargin=3)if(nargin=3)f=subs(f,t,x0);f=subs(f,t,x0);else else f=collect(f);%f=collect(f);%将插值多项式展开将插值多项式展开 f=vpa(f,6);f=vpa(f,6);
24、end end end endendend第22页，共49页，编辑于2022年，星期五例例4-4 根据下表的数据点求出其均差形式牛顿插值根据下表的数据点求出其均差形式牛顿插值多项式，并计算当多项式，并计算当x=2.0时时y的值。的值。解：解：x=1 1.2 1.8 2.5 4;x=1 1.2 1.8 2.5 4;y=1 1.44 3.24 6.25 16;y=1 1.44 3.24 6.25 16;f=Newton(x,y)f=Newton(x,y)f=f=.182711e-14-.482154e-14*t+1.00000*t2-.169177e-14*t3+.211471e-15*t4 .1
25、82711e-14-.482154e-14*t+1.00000*t2-.169177e-14*t3+.211471e-15*t4 f=Newton(x,y,2.0)f=Newton(x,y,2.0)f=f=4 4x x1 11.21.21.81.82.52.54 4 y y1 11.441.443.243.246.256.251616第23页，共49页，编辑于2022年，星期五4.4 利用差分的牛顿插值法利用差分的牛顿插值法n n差分分为向前差分、后向差分和中心差分三种，它们的记法及定义如下为：其中：代表向前差分；代表向后差分；代表向后差分。第24页，共49页，编辑于2022年，星期五 假设假
26、设 。为了方便，可构造如表。为了方便，可构造如表4-44-4所所示的差分表（示的差分表（）。）。表表4-4 4-4 差分计算表格差分计算表格第25页，共49页，编辑于2022年，星期五4.4.1向前牛顿插值向前牛顿插值n n向前牛顿插值向前牛顿插值多项式可表示如下：多项式可表示如下：其中其中 叫做步长，叫做步长，且且 的取值范围为的取值范围为 。第26页，共49页，编辑于2022年，星期五n n在MATLAB中编程实现向前牛顿插值法向前牛顿插值法函数为：NewtonforwardNewtonforward。n n功能：求已知数据点的向前牛顿插值法多项式；n n调用格式：f=f=Newtonfo
27、rward(x,y)Newtonforward(x,y)或或 f=f=NewtonforwardNewtonforward (x,y,x0)(x,y,x0)。其中，其中，x为已知数据点的为已知数据点的x x 坐标向量；坐标向量；y为已知数据点的为已知数据点的y y 坐标向量；坐标向量；x0 x0为插值点的为插值点的x x坐标；坐标；f f为求得的向前牛顿插值法多项式或为求得的向前牛顿插值法多项式或x0 x0处的插值。处的插值。第27页，共49页，编辑于2022年，星期五function f=Newtonforward(x,y,x0)function f=Newtonforward(x,y,x0
28、)%求已知数据点的向前差分牛顿插值多项式求已知数据点的向前差分牛顿插值多项式%已知数据点的已知数据点的x x 坐标向量坐标向量:x:x%已知数据点的已知数据点的y y 坐标向量坐标向量:y:y%为插值点的为插值点的x x坐标坐标:x0:x0%求得的向前差分牛顿插值多项式或求得的向前差分牛顿插值多项式或x0 x0处的插值处的插值:f:fsyms t;syms t;if(length(x)=length(y)if(length(x)=length(y)n=length(x);n=length(x);c(1:n)=0.0;c(1:n)=0.0;elseelse disp(x disp(x和和y y的
29、维数不相等！的维数不相等！););return;return;endend第28页，共49页，编辑于2022年，星期五f=y(1);f=y(1);y1=0;y1=0;xx=linspace(x(1),x(n),(x(2)-x(1);xx=linspace(x(1),x(n),(x(2)-x(1);if(xx=x)if(xx=x)disp(disp(节点之间不是等距的！节点之间不是等距的！););return;return;endendfor(i=1:n-1)for(i=1:n-1)for(j=1:n-i)for(j=1:n-i)y1(j)=y(j+1)-y(j);y1(j)=y(j+1)-y(
30、j);end end c(i)=y1(1);c(i)=y1(1);l=t;l=t;for(k=1:i-1)for(k=1:i-1)l=l*(t-k);l=l*(t-k);end;end;f=f+c(i)*l/factorial(i);f=f+c(i)*l/factorial(i);simplify(f);simplify(f);y=y1;y=y1;if(i=n-1)if(i=n-1)if(nargin=3)if(nargin=3)f=subs(f,t,(x0-x(1)/(x(2)-x(1);f=subs(f,t,(x0-x(1)/(x(2)-x(1);else else f=collect(f
31、);f=collect(f);f=vpa(f,6);f=vpa(f,6);end end end endendend第29页，共49页，编辑于2022年，星期五4.4.2向前牛顿插值向前牛顿插值n n向后牛顿插值向后牛顿插值多项式可表示如下：多项式可表示如下：其中其中 叫做步长，叫做步长，且且 的取值范围为的取值范围为 。第30页，共49页，编辑于2022年，星期五n n在MATLAB中编程实现向后牛顿插值法向后牛顿插值法函数为：函数为：NewtonbackNewtonback。n n功能：求已知数据点的向前牛顿插值法多项式；功能：求已知数据点的向前牛顿插值法多项式；n n调用格式：调用格式：
32、f=f=NewtonbackNewtonback (x,y)或 f=f=Newtonback(x,y,x0)Newtonback(x,y,x0)。其中，其中，x为已知数据点的x x 坐标向量；y y为已知数据点的为已知数据点的y 坐标向量；x0 x0为插值点的为插值点的x x坐标；f为求得的向前牛顿插值法多项式或为求得的向前牛顿插值法多项式或x0 x0处的插值。第31页，共49页，编辑于2022年，星期五function f=Newtonback(x,y,x0)function f=Newtonback(x,y,x0)%求已知数据点的向后差分牛顿插值多项式求已知数据点的向后差分牛顿插值多项式%
33、已知数据点的已知数据点的x x 坐标向量坐标向量:x:x%已知数据点的已知数据点的y y 坐标向量坐标向量:y:y%为插值点的为插值点的x x坐标坐标:x0:x0%求得的向前差分牛顿插值多项式或求得的向前差分牛顿插值多项式或x0 x0处的插值处的插值:f:fsyms t;syms t;if(length(x)=length(y)if(length(x)=length(y)n=length(x);n=length(x);c(1:n)=0.0;c(1:n)=0.0;elseelse disp(x disp(x和和y y的维数不相等！的维数不相等！););return;return;endend第3
34、2页，共49页，编辑于2022年，星期五f=y(n);f=y(n);y1=0;y1=0;xx=linspace(x(1),x(n),(x(2)-x(1);xx=linspace(x(1),x(n),(x(2)-x(1);if(xx=x)if(xx=x)disp(disp(节点之间不是等距的！节点之间不是等距的！););return;return;endendfor(i=1:n-1)for(i=1:n-1)for(j=i+1:n)for(j=i+1:n)y1(j)=y(j)-y(j-1);y1(j)=y(j)-y(j-1);end end c(i)=y1(n);c(i)=y1(n);l=t;l=
35、t;for(k=1:i-1)for(k=1:i-1)l=l*(t+k);l=l*(t+k);end;end;f=f+c(i)*l/factorial(i);f=f+c(i)*l/factorial(i);simplify(f);simplify(f);y=y1;y=y1;if(i=n-1)if(i=n-1)if(nargin=3)if(nargin=3)f=subs(f,t,(x0-x(n)/(x(2)-x(1);f=subs(f,t,(x0-x(n)/(x(2)-x(1);else else f=collect(f);f=collect(f);f=vpa(f,6);f=vpa(f,6);en
36、d end end endendend第33页，共49页，编辑于2022年，星期五例例4-5 根据下表的数据点求出其差分形式的牛顿根据下表的数据点求出其差分形式的牛顿插值多项式，并计算当插值多项式，并计算当x=1.55时时y的值。的值。解：解：x=1 1.2 1.4 1.6 1.8;x=1 1.2 1.4 1.6 1.8;y=0.8415 0.9320 y=0.8415 0.93200.9854 0.9996 0.9738;0.9854 0.9996 0.9738;f=Newtonforward(x,y)f=Newtonforward(x,y)f=f=.841500+.108025*t-.16
37、9042e-1*t2-.675000e-3*t3+.541667e-4*t4.841500+.108025*t-.169042e-1*t2-.675000e-3*t3+.541667e-4*t4 f=Newtonforward(x,y,1.55)f=Newtonforward(x,y,1.55)f=f=0.9998 0.9998f=Newtonback(x,y)f=Newtonback(x,y)f=f=.973800-.457417e-1*t-.198042e-1*t2+.191667e-3*t3+.541667e-4*t4.973800-.457417e-1*t-.198042e-1*t2+
38、.191667e-3*t3+.541667e-4*t4 f=Newtonback(x,y,1.55)f=Newtonback(x,y,1.55)f=f=0.99980.9998x x1 11.21.21.41.41.61.61.81.8 y y0.84150.84150.93200.93200.98540.98540.99960.99960.97380.9738第34页，共49页，编辑于2022年，星期五4.5 Hermite插值插值n nHermiteHermite插值插值满足在节点上等于给定函数值，而且在节点上的满足在节点上等于给定函数值，而且在节点上的导数值也等于给定的导数值。对于高阶导
39、数的情况，导数值也等于给定的导数值。对于高阶导数的情况，HermiteHermite插值插值多项式比较复杂，在实际中，常常遇到的是函数多项式比较复杂，在实际中，常常遇到的是函数值与一阶导数给定的情况。在此情况下，值与一阶导数给定的情况。在此情况下，n n个节点个节点x x1 1，x x2 2，x xn n的的HermiteHermite插值插值多项式的表达式如下：多项式的表达式如下：其中其中 ，第35页，共49页，编辑于2022年，星期五n n在MATLAB中编程实现HermiteHermite插值法插值法函数为：函数为：Hermite。n n功能：求已知数据点的功能：求已知数据点的Hermi
40、te插值法多项式；插值法多项式；n n调用格式：f=f=Hermite (x,y,y_1)(x,y,y_1)或或 f=f=HermiteHermite (x,y,y_1,x0)(x,y,y_1,x0)。其中，其中，x为已知数据点的为已知数据点的x x 坐标向量；坐标向量；y y为已知数据点的为已知数据点的y y 坐标向量；坐标向量；y_1y_1为已知数据点导数向量；为已知数据点导数向量；x0 x0为插值点的为插值点的x x坐标；坐标；f f为求得的为求得的HermiteHermite插值法多项式或插值法多项式或x0 x0处的插值。第36页，共49页，编辑于2022年，星期五function f
41、=Hermite(x,y,y_1,x0)function f=Hermite(x,y,y_1,x0)%求已知数据点的向后差分牛顿插值多项式求已知数据点的向后差分牛顿插值多项式%已知数据点的已知数据点的x x 坐标向量坐标向量:x:x%已知数据点的已知数据点的y y 坐标向量坐标向量:y:y%已知数据点的导数向量已知数据点的导数向量:y_1:y_1%求得的求得的HermiteHermite插值多项式或插值多项式或x0 x0处的插值处的插值:f:fsyms t;syms t;f=0.0;f=0.0;if(length(x)=length(y)if(length(x)=length(y)if(len
42、gth(y)=length(y_1)if(length(y)=length(y_1)n=length(x);n=length(x);else else disp(y disp(y和和y y的导数的维数不相等！的导数的维数不相等！););return;return;end endelseelse disp(x disp(x和和y y的维数不相等！的维数不相等！););return;return;endend第37页，共49页，编辑于2022年，星期五for i=1:nfor i=1:n h=1.0;h=1.0;a=0.0;a=0.0;for j=1:n for j=1:n if(j=i)if(j
43、=i)h=h*(t-x(j)2/(x(i)-x(j)2);h=h*(t-x(j)2/(x(i)-x(j)2);a=a+1/(x(i)-x(j);a=a+1/(x(i)-x(j);end end end end f=f+h*(x(i)-t)*(2*a*y(i)-y_1(i)+y(i);f=f+h*(x(i)-t)*(2*a*y(i)-y_1(i)+y(i);if(i=n)if(i=n)if(nargin=4)if(nargin=4)f=subs(f,t,x0);f=subs(f,t,x0);else else f=vpa(f,6);f=vpa(f,6);end end end endendend
44、第38页，共49页，编辑于2022年，星期五例例4-6 根据下表的数据点求出其根据下表的数据点求出其Hermite插值多项插值多项式，并计算当式，并计算当x=1.144时时y的值。的值。解：解：x=1:0.2:1.8;x=1:0.2:1.8;y=1 1.0954 1.1832 1.2649 1.3416;y=1 1.0954 1.1832 1.2649 1.3416;y_1=0.5 0.4564 0.4226 0.3953 0.3727;y_1=0.5 0.4564 0.4226 0.3953 0.3727;f=Hermite(x,y,y_1)f=Hermite(x,y,y_1)f=f=678
45、.168*(t-1.20000)2*(t-1.40000)2*(t-1.60000)2*(t-1.80000)2*(-20.3333+21.3333*t)+10850.7*(t-1.)2*(t-1.40000)2*(t-678.168*(t-1.20000)2*(t-1.40000)2*(t-1.60000)2*(t-1.80000)2*(-20.3333+21.3333*t)+10850.7*(t-1.)2*(t-1.40000)2*(t-1.60000)2*(t-1.80000)2*(-10.4063+9.58473*t)+24414.1*(t-1.)2*(t-1.20000)2*(t-1
46、.60000)2*(t-1.60000)2*(t-1.80000)2*(-10.4063+9.58473*t)+24414.1*(t-1.)2*(t-1.20000)2*(t-1.60000)2*(t-1.80000)2*(.591560+.422600*t)+10850.7*(t-1.)2*(t-1.20000)2*(t-1.40000)2*(t-1.80000)2*(17.4978-10.1455*t)+678.168*(t-1.80000)2*(.591560+.422600*t)+10850.7*(t-1.)2*(t-1.20000)2*(t-1.40000)2*(t-1.80000)
47、2*(17.4978-10.1455*t)+678.168*(t-1.)2*(t-1.20000)2*(t-1.40000)2*(t-1.60000)2*(50.9807-27.5773*t)1.)2*(t-1.20000)2*(t-1.40000)2*(t-1.60000)2*(50.9807-27.5773*t)f=Hermite(x,y,y_1,1.44)f=Hermite(x,y,y_1,1.44)f=f=1.2000 1.2000 x x1 11.21.21.41.41.61.61.81.8 y y1 11.09541.09541.18321.18321.26491.26491.34
48、161.3416y y 0.50000.50000.45640.45640.42260.42260.39530.39530.37270.3727第39页，共49页，编辑于2022年，星期五4.6 spline三次样条插值三次样条插值n nSpline插值为分段三次插值，即：在每一个小区间上是不超过三次的多项式且具有二阶连续导数，利用具有一阶导数边界条件的源程序为：naspline.mfunction m=naspline(x,y,dy0,dyn,xx)function m=naspline(x,y,dy0,dyn,xx)%用途：三阶样条插值（一阶导数边界条件）用途：三阶样条插值（一阶导数边界条
49、件）%格式：格式：x x为节点向量；为节点向量；y y为数据；为数据；dyodyo，dyndyn为左右两端点的一阶导数为左右两端点的一阶导数%如果如果xxxx缺省，则输出各节点的一阶导数值，否则缺省，则输出各节点的一阶导数值，否则mm为为xxxx的三阶样条插值的三阶样条插值n=length(x)-1;%n=length(x)-1;%计算小区间的个数计算小区间的个数h=diff(x);h=diff(x);lemda=h(2:n)./(h(1:n-1)+h(2:n);lemda=h(2:n)./(h(1:n-1)+h(2:n);mu=1-lemda;mu=1-lemda;第40页，共49页，编辑于
50、2022年，星期五g=3*(lemda.*diff(y(1:n)./h(1:n-1)+mu.*diff(y(2:n+1)./h(2:n);g=3*(lemda.*diff(y(1:n)./h(1:n-1)+mu.*diff(y(2:n+1)./h(2:n);g(1)=g(1)-lemda(1)*dy0;g(1)=g(1)-lemda(1)*dy0;g(n-1)=g(n-1)-mu(n-1)*dyn;g(n-1)=g(n-1)-mu(n-1)*dyn;%求解三对角方程求解三对角方程dy=nachase(lemda,2*ones(1:n-1),mu,g);dy=nachase(lemda,2*on