Matlab数据插值与拟合.ppt

数据插值与拟合,在工程实践与科学实验中,常常需要从一组试验数据之中找到自变量与因变量之间的关系，一般可用一个近似函数表示。函数产生的办法因观测数据的要求不同而异，数据插值与拟合是两种常用的方法。4.1 MATLAB中的插值函数4.2 拉格朗日插值法4.3 利用均差的牛顿插值法4.4 利用差分的牛顿插值法4.5 Hermite插值4.6 spline三次样条插值4.7 多项式曲线拟合4.8 最小二乘拟合,4.1 MATLAB中的插值函数,函数插值来源于函数的以下问题：只知道函数在某区间有定义且已得到区间内一些离散点的值，希望用简单的表达式近似给出函数在此区间上的整体描述，并能与已知离散点上的值相等。插值法按插值函数的形式主要分为以下几种形式：（1）代数多项式插值；（2）三角多项式插值；（3）有理分式插值。,代数多项式插值是最常用的插值方式，其内容也是最丰富的，它又可分为以下几种插值方式：（1）非等距节点插值，包括拉格朗日插值、利用均差的牛顿插值和埃特金插值；（2）非等距节点插值，包括利用差分的牛顿插值和高斯插值等；（3）在插值中增加了导数的Hermite（埃尔米特）插值；（4）分段插值，包括分段线性插值、分段Hermite（埃尔米特）插值和样条函数插值；（5）反插值。按被插值函数的变量个数还可把插值法分为一元插值和多元插值。,4.1.1 一元插值函数,MATLAB中的一元插值函数为interp1()，它的功能是一维数据插值（表格查找）。该命令对数据点之间进行计算内插值，它出一元函数f(x)在中间点的数值，其中函数f(x)由所给数据决定。一元插值函数interp1()的几种调用格式如表4-1所示。表4-1 一维插值插值函数interp1的语法格式,MATLAB中一维插值有多种算法，由interp1函数中的method指定。MATLAB中一维插值的各种算法如表4-2所示。表4-2 一维插值算法(method),1.Linear（分段线性插值）,它的算法是在每个小区间xi,xi+1上采用简单的线性插值。在区间xi,xi+1上的子插值多项式为：由此整个区间xi,xi+1上的插值函数为：其中 定义如下：,分段线性插值方法在速度和误差之间取得了比较好的均衡，其插值函数具有连续性，但在已知数据点处的斜率一般不会改变，因此不是光滑的。分段线性插值方法是MATLAB一维插值默认的方法。,2.Spline（样条插值）,样条插值是用分段低次多项式去逼近函数。样条函数可以给出光滑 的插值曲线，只要在插值区间端点提供某些导数信息，样条插值可以适应不同光滑需求。三次样条是使用最为广泛的样条插值，它在每个子区间xi,xi+1上都是有二阶连续导数的三次多项式，即 其中 都是三次多项式。,对于给定的离散的测量数据经x，y（称为断点），要寻找一个三次多项式y=p(x)，以逼近每对数据(xi,yi)点间曲线。过两点(xi,yi)和(xi+1,yi+1)只能确定一条直线，而通过一点的三次多项式曲线有无穷多条。为使通过中间断点的三次多项式曲线具有唯一性，要增加以下的连续条件和边界条件（因为三次多项式有4个系数）：（1）三次多项式在点(xi,yi)处有：；（2）三次多项式在点(xi,yi)处有：；（3）三次多项式在点(xi,yi)处有：；（4）边界条件：。,表4-2中各种方法中：（1）nearest方法速度最快，占用内存最小，但一般来说误差最大，插值结果最不光滑；（2）spline三次样条插值是所有插值方法中运行耗时最长的，其插值函数以及插值函数的一阶、二阶导函数都连续，因此是最光滑的插值方法，占用内存上比cubic方法小，但当已知数据点不均匀分布时可能出现异常结果。（3）cubic三次多项式插值法中插值函数及其一阶导数都是连续的，因此其插值结果也比较光滑，运算速度比spline方法略快，但占用内存最多。在实际的使用中，应根据实际需求和运算条件选择合适的算法。,例4-1 用interp1对sin函数进行分段线性插值。,解：在MATLAB命令窗口中输入以下命令：x=0:2*pi;y=sin(x);xx=0:0.5:2*pi yy=interp1(x,y,xx);plot(x,y,s,xx,yy)注：例4-1中用默认的(分段线性插值的linear)对已知的7个sin函数的数据点进行插值，用plot画出插值结果。从图中可以看出分段线性就是联结两个邻近的已知点的线性函数插值计算该区间内插值点上的函数值。,例4-2 用其他一维插值方法对以下7个离散数据点(1,3.5)、(2,2.1)、(3,1.3)、(4.0.8)、(5,2.9)、(6,4.2)、(7,5.7)进行一维插值方法。,解：在MATLAB命令窗口中输入以下命令：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,:),4.2 拉格朗日插值法,拉格朗日插值法是基于基函数的插值方法，插值多项式可表示为 其中 称为i次基函数：,在MATLAB中编程实现拉格朗日插值法函数为：Language。功能：求已知数据点的拉格朗日多项式；调用格式：f=Language(x,y)或f=Language(x,y,x0)。其中，x为已知数据点的x 坐标向量；y为已知数据点的y 坐标向量；x0为插值点的x坐标；f为求得的拉格朗日多项式或x0处的插值。,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)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,例4-3 根据下表的数据点求出其拉格朗日插值多项式，并计算当x=1.6时y的值。,解：x=1 1.2 1.8 2.5 4;y=0.8415 0.9320 0.9738 0.5985-0.7568;f=language(x,y)f=1.05427*t-.145485e-1*t2-.204917*t3+.328112e-1*t4-.261189e-1 f=language(x,y,1.6)f=0.9992,4.3 利用均差的牛顿插值法,函数 f的零阶均差定义为，一阶定义均差为 一般地，函数f的k阶均差定义为：利用均差的牛顿插值法多项式为,系数的计算过程如表4-3所示,表4-3 均差计算表格,在MATLAB中编程实现利用均差牛顿插值法函数为：Newton。功能：求已知数据点的均差形式的牛顿插值多项式；调用格式：f=Newton(x,y)或f=Newton(x,y,x0)。其中，x为已知数据点的x 坐标向量；y为已知数据点的y 坐标向量；x0为插值点的x坐标；f为求得的牛顿插值法多项式或x0处的插值。,function f=Newton(x,y,x0)%求已知数据点的均差形式牛顿插值多项式%已知数据点的x 坐标向量:x%已知数据点的y 坐标向量:y%为插值点的x坐标:x0%求得的均差形式牛顿插值多项式或x0处的插值:fsyms t;if(length(x)=length(y)n=length(x);c(1:n)=0.0;else disp(x和y的维数不相等！);return;end,f=y(1);y1=0;l=1;for(i=1:n-1)for(j=i+1:n)y1(j)=(y(j)-y(i)/(x(j)-x(i);end c(i)=y1(i+1);l=l*(t-x(i);f=f+c(i)*l;simplify(f);y=y1;if(i=n-1)if(nargin=3)f=subs(f,t,x0);else f=collect(f);%将插值多项式展开 f=vpa(f,6);end endend,例4-4 根据下表的数据点求出其均差形式牛顿插值多项式，并计算当x=2.0时y的值。,解：x=1 1.2 1.8 2.5 4;y=1 1.44 3.24 6.25 16;f=Newton(x,y)f=.182711e-14-.482154e-14*t+1.00000*t2-.169177e-14*t3+.211471e-15*t4 f=Newton(x,y,2.0)f=4,4.4 利用差分的牛顿插值法,差分分为向前差分、后向差分和中心差分三种，它们的记法及定义如下为：其中：代表向前差分；代表向后差分；代表向后差分。,假设。为了方便，可构造如表4-4所示的差分表（）。表4-4 差分计算表格,向前牛顿插值,向前牛顿插值多项式可表示如下：其中 叫做步长，且 的取值范围为。,在MATLAB中编程实现向前牛顿插值法函数为：Newtonforward。功能：求已知数据点的向前牛顿插值法多项式；调用格式：f=Newtonforward(x,y)或 f=Newtonforward(x,y,x0)。其中，x为已知数据点的x 坐标向量；y为已知数据点的y 坐标向量；x0为插值点的x坐标；f为求得的向前牛顿插值法多项式或x0处的插值。,function f=Newtonforward(x,y,x0)%求已知数据点的向前差分牛顿插值多项式%已知数据点的x 坐标向量:x%已知数据点的y 坐标向量:y%为插值点的x坐标:x0%求得的向前差分牛顿插值多项式或x0处的插值:fsyms t;if(length(x)=length(y)n=length(x);c(1:n)=0.0;else disp(x和y的维数不相等！);return;end,f=y(1);y1=0;xx=linspace(x(1),x(n),(x(2)-x(1);if(xx=x)disp(节点之间不是等距的！);return;endfor(i=1:n-1)for(j=1:n-i)y1(j)=y(j+1)-y(j);end c(i)=y1(1);l=t;for(k=1:i-1)l=l*(t-k);end;f=f+c(i)*l/factorial(i);simplify(f);y=y1;if(i=n-1)if(nargin=3)f=subs(f,t,(x0-x(1)/(x(2)-x(1);else f=collect(f);f=vpa(f,6);end endend,向前牛顿插值,向后牛顿插值多项式可表示如下：其中 叫做步长，且 的取值范围为。,在MATLAB中编程实现向后牛顿插值法函数为：Newtonback。功能：求已知数据点的向前牛顿插值法多项式；调用格式：f=Newtonback(x,y)或 f=Newtonback(x,y,x0)。其中，x为已知数据点的x 坐标向量；y为已知数据点的y 坐标向量；x0为插值点的x坐标；f为求得的向前牛顿插值法多项式或x0处的插值。,function f=Newtonback(x,y,x0)%求已知数据点的向后差分牛顿插值多项式%已知数据点的x 坐标向量:x%已知数据点的y 坐标向量:y%为插值点的x坐标:x0%求得的向前差分牛顿插值多项式或x0处的插值:fsyms t;if(length(x)=length(y)n=length(x);c(1:n)=0.0;else disp(x和y的维数不相等！);return;end,f=y(n);y1=0;xx=linspace(x(1),x(n),(x(2)-x(1);if(xx=x)disp(节点之间不是等距的！);return;endfor(i=1:n-1)for(j=i+1:n)y1(j)=y(j)-y(j-1);end c(i)=y1(n);l=t;for(k=1:i-1)l=l*(t+k);end;f=f+c(i)*l/factorial(i);simplify(f);y=y1;if(i=n-1)if(nargin=3)f=subs(f,t,(x0-x(n)/(x(2)-x(1);else f=collect(f);f=vpa(f,6);end endend,例4-5 根据下表的数据点求出其差分形式的牛顿插值多项式，并计算当x=1.55时y的值。,解：x=1 1.2 1.4 1.6 1.8;y=0.8415 0.93200.9854 0.9996 0.9738;f=Newtonforward(x,y)f=.841500+.108025*t-.169042e-1*t2-.675000e-3*t3+.541667e-4*t4 f=Newtonforward(x,y,1.55)f=0.9998f=Newtonback(x,y)f=.973800-.457417e-1*t-.198042e-1*t2+.191667e-3*t3+.541667e-4*t4 f=Newtonback(x,y,1.55)f=0.9998,4.5 Hermite插值,Hermite插值满足在节点上等于给定函数值，而且在节点上的导数值也等于给定的导数值。对于高阶导数的情况，Hermite插值多项式比较复杂，在实际中，常常遇到的是函数值与一阶导数给定的情况。在此情况下，n个节点x1，x2，xn的Hermite插值多项式的表达式如下：其中，,在MATLAB中编程实现Hermite插值法函数为：Hermite。功能：求已知数据点的Hermite插值法多项式；调用格式：f=Hermite(x,y,y_1)或 f=Hermite(x,y,y_1,x0)。其中，x为已知数据点的x 坐标向量；y为已知数据点的y 坐标向量；y_1为已知数据点导数向量；x0为插值点的x坐标；f为求得的Hermite插值法多项式或x0处的插值。,function f=Hermite(x,y,y_1,x0)%求已知数据点的向后差分牛顿插值多项式%已知数据点的x 坐标向量:x%已知数据点的y 坐标向量:y%已知数据点的导数向量:y_1%求得的Hermite插值多项式或x0处的插值:fsyms t;f=0.0;if(length(x)=length(y)if(length(y)=length(y_1)n=length(x);else disp(y和y的导数的维数不相等！);return;endelse disp(x和y的维数不相等！);return;end,for i=1:n h=1.0;a=0.0;for j=1:n if(j=i)h=h*(t-x(j)2/(x(i)-x(j)2);a=a+1/(x(i)-x(j);end end f=f+h*(x(i)-t)*(2*a*y(i)-y_1(i)+y(i);if(i=n)if(nargin=4)f=subs(f,t,x0);else f=vpa(f,6);end endend,例4-6 根据下表的数据点求出其Hermite插值多项式，并计算当x=1.144时y的值。,解：x=1:0.2:1.8;y=1 1.0954 1.1832 1.2649 1.3416;y_1=0.5 0.4564 0.4226 0.3953 0.3727;f=Hermite(x,y,y_1)f=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.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.)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=1.2000,4.6 spline三次样条插值,Spline插值为分段三次插值，即：在每一个小区间上是不超过三次的多项式且具有二阶连续导数，利用具有一阶导数边界条件的源程序为：naspline.mfunction m=naspline(x,y,dy0,dyn,xx)%用途：三阶样条插值（一阶导数边界条件）%格式：x为节点向量；y为数据；dyo，dyn为左右两端点的一阶导数%如果xx缺省，则输出各节点的一阶导数值，否则m为xx的三阶样条插值n=length(x)-1;%计算小区间的个数h=diff(x);lemda=h(2:n)./(h(1:n-1)+h(2:n);mu=1-lemda;,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(n-1)=g(n-1)-mu(n-1)*dyn;%求解三对角方程dy=nachase(lemda,2*ones(1:n-1),mu,g);%若给插值节点，计算插值m=dy0;dy;dyn;if nargin=5 s=zeros(size(xx);for i=1:n if i=1 kk=find(xx=x(n);else kk=find(xxx(i)end,%追赶法function x=nachase(a,b,c,d)n=length(a);for k=2:n b(k)=b(k)-a(k)/b(k-1)*c(k-1);d(k)=d(k)-a(k)/b(k-1)*d(k-1);endx(n)=d(n)/b(n);for k=n-1:-1:1 x(k)=(d(k)-c(k)*x(k+1)/b(k);endx=x(:);%基函数function y=alpha0(x)y=2*x.3-3*x.2+1;function y=alpha1(x)y=-2*x.3+3*x.2;function y=beta0(x)y=x.3-2*x.2+x;function y=beta1(x)y=x.3-x.2;,4.7 多项式曲线拟合,对给定的试验数据点(xi,yi)(i=1,2,N)，可以构造m次多项式：由曲线拟合的定义，应该使得下式取极小值 通过简单运算可以得出系数是下面线性方程组的解,其中,在MATLAB中编程实现多项式曲线拟合函数为：multifit。功能：求已知数据点的多项式曲线拟合插值法多项式；调用格式：A=multifit(x,y,m)。其中，x为已知数据点的x 坐标向量；y为已知数据点的y 坐标向量；m为拟合多项式的次数；A为拟合多项式的系数向量。,function A=multifit(X,Y,m)%离散数据点的多项式曲线拟合%试验数据点x坐标向量:x%试验数据点y坐标向量:y%拟合多项式的次数:m%A-输出的拟合多项式的系数N=length(X);M=length(Y);if(N=M)disp(数据点坐标不匹配！);return;endc(1:(2*m+1)=0;b(1:(m+1)=0;for j=1:(2*m+1)%求出c和b for k=1:N c(j)=c(j)+X(k)(j-1);if(j(m+2)b(j)=b(j)+Y(k)*X(k)(j-1);end endendC(1,:)=c(1:(m+1);for s=2:(m+1)C(s,:)=c(s:(m+s);endA=bC;%直接求解法求出拟合系数,例4-7 用二次多项式拟合下表中所列的数据点。,解：x=1:3;y=2 5 10;A=multifit(x,y,2)A=0.1282 0.3235 0.8718即拟合的多项式为p=0.1282+0.3235 x+0.8718 x2。,4.8 最小二乘拟合,最小二乘法拟合在科学实验的统计方法中经常使用。它的具体操作过程是从一组实验数据(xi,yi)中拟合出函数关系y=f(x)，拟合的标准是使(f(xi)-yi)的平方取极小值。在MATLAB中可以使用polyfit函数对数据进行最小二乘法拟合，它的基本调用格式为：p=polyfit(X,Y,N)表示用N次多项式拟合数据点(xi,yi)。,例4-8 用四次多项式拟合下表中所列的数据点。,解：x=1:0.2:1.8;y=1 1.09541.18321.26491.3416;p=polyfit(x,y,4)p=-0.0104 0.0854-0.3121 0.9086 0.3285即拟合的多项式为p=-0.0104+0.0854x-0.3121x2+0.9086 x3+0.3285x4。,