Matlab数据插值与拟合

数据插值与拟合在工程实践与科学实验中,常常需要从一组试验数据之中找到自变量与因变量之间的关系，一般可用一个近似函数表示。函数产生的方法因观测数据的要求不同而异，数据插值与拟合是两种常用的方法。4.1MATLAB中的插值函数4.2拉格朗日插值法4.3利用均差的牛顿插值法4.4利用差分的牛顿插值法4.5Hermite插值4.6spline三次样条插值4.7多项式曲线拟合4.8最小二乘拟合4.1MATLAB中的插值函数函数插值来源于函数的以下问题：只知道函数在某区间有定义且已得到区间内一些离散点的值，希望用简单的表达式近似给出函数在此区间上的整体描述，并能与离散点上的值相等。插值法按插值函数的形式主要分为以下几种形式：〔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的语法格式语法形式说明y=interp1(x,Y,xi)由已知点集(x,Y）插值计算xi上的函数值y=interp1(x,Y,xi)相当于x=1:length(Y)的interp(x,Y,xi)y=interp1(x,Y,xi,method)用指定插值方法计算插值点xi上的函数值y=interp1(x,Y,xi,method,’extrap’)对xi中超出已知点集的插值点用指定插值方法计算函数值y=interp1(x,Y,xi,method,’extrap’,extrapval)用指定方法插值xi上的函数值，超出已知点集处函数值取extrapvaly=interp1(x,Y,xi,method,’pp’)用指定方法插值，但返回结果为分段多项式method方法描述‘nearest’最邻近插值：插值点处函数值取与插值点最邻近的已知点的函数值‘liner’分段线性插值：插值点处函数值由连接其最邻近的两侧点的线性函数预测，MATLAB中interp1的默认方法‘spline’样条插值：默认为三次样条插值。可用spline函数代替‘pchip’三次Hermite多项式插值。可用pchip函数代替‘cubic’同‘pchip’，三次Hermite多项式插值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=[1234567];>>y=[3.52.11.30.82.94.25.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处的插值。functionf=Language(x,y,x0)%求数据点的拉格朗日多项式%数据点的x坐标向量:x%数据点的y坐标向量:y%为插值点的x坐标:x0%求得的拉格朗日多项式或x0处的插值:fsymst;if(length(x)==length(y))n=length(x);elsedisp('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);%计算插值点的函数值elsef=collect(f);%将插值多项式展开f=vpa(f,6);%将插值多项式的系数化成6位精度的小数endendend例4-3根据下表的数据点求出其拉格朗日插值多项式，并计算当x=1.6时y的值。解：>>

x=[11.21.82.54];>>y=[0.84150.93200.97380.5985-0.7568];>>f=language(x,y)f=1.05427*t-.145485e-1*t^2-.204917*t^3+.328112e-1*t^4-.261189e-1>>f=language(x,y,1.6)f=0.9992x11.21.82.54y0.84150.93200.97380.5985-0.75684.3利用均差的牛顿插值法函数f的零阶均差定义为，一阶定义均差为一般地，函数f的k阶均差定义为：利用均差的牛顿插值法多项式为

系数的计算过程如表4-3所示表4-3均差计算表格一阶均差二阶均差三阶均差……n阶均差……………在MATLAB中编程实现利用均差牛顿插值法函数为：Newton。功能：求数据点的均差形式的牛顿插值多项式；调用格式：f=Newton(x,y)或f=Newton(x,y,x0)。其中，x为数据点的x坐标向量；y为数据点的y坐标向量；x0为插值点的x坐标；f为求得的牛顿插值法多项式或x0处的插值。functionf=Newton(x,y,x0)%求数据点的均差形式牛顿插值多项式%数据点的x坐标向量:x%数据点的y坐标向量:y%为插值点的x坐标:x0%求得的均差形式牛顿插值多项式或x0处的插值:fsymst;if(length(x)==length(y))n=length(x);c(1:n)=0.0;elsedisp('x和y的维数不相等！');return;endf=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));endc(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);elsef=collect(f);%将插值多项式展开

f=vpa(f,6);endendend例4-4根据下表的数据点求出其均差形式牛顿插值多项式，并计算当x=2.0时y的值。解：>>x=[11.21.82.54];>>y=[11.443.246.2516];>>f=Newton(x,y)f=.182711e-14-.482154e-14*t+1.00000*t^2-.169177e-14*t^3+.211471e-15*t^4>>f=Newton(x,y,2.0)f=4x11.21.82.54y11.443.246.25164.4利用差分的牛顿插值法差分分为向前差分、后向差分和中心差分三种，它们的记法及定义如下为：

其中：代表向前差分；代表向后差分；代表向后差分。假设。为了方便，可构造如表4-4所示的差分表〔〕。表4-4差分计算表格………………向前牛顿插值向前牛顿插值多项式可表示如下：其中叫做步长，，且的取值范围为。在MATLAB中编程实现向前牛顿插值法函数为：Newtonforward。功能：求数据点的向前牛顿插值法多项式；调用格式：f=Newtonforward(x,y)或f=Newtonforward(x,y,x0)。其中，x为数据点的x坐标向量；y为数据点的y坐标向量；x0为插值点的x坐标；f为求得的向前牛顿插值法多项式或x0处的插值。functionf=Newtonforward(x,y,x0)%求数据点的向前差分牛顿插值多项式%数据点的x坐标向量:x%数据点的y坐标向量:y%为插值点的x坐标:x0%求得的向前差分牛顿插值多项式或x0处的插值:fsymst;if(length(x)==length(y))n=length(x);c(1:n)=0.0;elsedisp('x和y的维数不相等！');return;endf=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);endc(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)));elsef=collect(f);f=vpa(f,6);endendend向前牛顿插值向后牛顿插值多项式可表示如下：其中叫做步长，，且的取值范围为。在MATLAB中编程实现向后牛顿插值法函数为：Newtonback。功能：求数据点的向前牛顿插值法多项式；调用格式：f=Newtonback(x,y)或f=Newtonback(x,y,x0)。其中，x为数据点的x坐标向量；y为数据点的y坐标向量；x0为插值点的x坐标；f为求得的向前牛顿插值法多项式或x0处的插值。functionf=Newtonback(x,y,x0)%求数据点的向后差分牛顿插值多项式%数据点的x坐标向量:x%数据点的y坐标向量:y%为插值点的x坐标:x0%求得的向前差分牛顿插值多项式或x0处的插值:fsymst;if(length(x)==length(y))n=length(x);c(1:n)=0.0;elsedisp('x和y的维数不相等！');return;endf=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);endc(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)));elsef=collect(f);f=vpa(f,6);endendend例4-5根据下表的数据点求出其差分形式的牛顿插值多项式，并计算当x=1.55时y的值。解：>>x=[11.21.41.61.8];>>y=[0.84150.9320 0.98540.99960.9738];>>f=Newtonforward(x,y)f=.841500+.108025*t-.169042e-1*t^2-.675000e-3*t^3+.541667e-4*t^4>>f=Newtonforward(x,y,1.55)f=0.9998f=Newtonback(x,y)f=.973800-.457417e-1*t-.198042e-1*t^2+.191667e-3*t^3+.541667e-4*t^4>>f=Newtonback(x,y,1.55)f=0.9998x11.21.41.61.8y0.84150.93200.98540.99960.97384.5Hermite插值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处的插值。functionf=Hermite(x,y,y_1,x0)%求数据点的向后差分牛顿插值多项式%数据点的x坐标向量:x%数据点的y坐标向量:y%数据点的导数向量:y_1%求得的Hermite插值多项式或x0处的插值:fsymst;f=0.0;if(length(x)==length(y))if(length(y)==length(y_1))n=length(x);elsedisp('y和y的导数的维数不相等！');return;endelsedisp('x和y的维数不相等！');return;endfori=1:nh=1.0;a=0.0;forj=1:nif(j~=i)h=h*(t-x(j))^2/((x(i)-x(j))^2);a=a+1/(x(i)-x(j));endend

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);elsef=vpa(f,6);endendend例4-6根据下表的数据点求出其Hermite插值多项式，并计算当x=1.144时y的值。解：

>>x=1:0.2:1.8;>>y=[11.09541.18321.26491.3416];>>y_1=[0.50.45640.42260.39530.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.2000x11.21.41.61.8y11.09541.18321.26491.3416y’0.50000.45640.42260.39530.37274.6spline三次样条插值Spline插值为分段三次插值，即：在每一个小区间上是不超过三次的多项式且具有二阶连续导数，利用具有一阶导数边界条件的源程序为：naspline.mfunctionm=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];ifnargin>=5s=zeros(size(xx));fori=1:nifi==1kk=find(xx<=x(2));elseifi==nkk=find(xx>=x(n));elsekk=find(xx>x(i)&xx<=x(i+1));endxbar=(xx(kk)-x(i))/h(i);s(kk)=alpha0(xbar)*y(i)+alpha1(xbar)*y(i+1)+...h(i)*beta0(xbar)*m(i)+h(i)*beta1(xbar)*m(i+1);endm=s;end%追赶法functionx=nachase(a,b,c,d)n=length(a);fork=2:nb(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);fork=n-1:-1:1x(k)=(d(k)-c(k)*x(k+1)/b(k));endx=x(:);%基函数functiony=alpha0(x)y=2*x.^3-3*x.^2+1;functiony=alpha1(x)y=-2*x.^3+3*x.^2;functiony=beta0(x)y=x.^3-2*x.^2+x;functiony=beta1(x)y=x.^3-x.^2;4.7多项式曲线拟合对给定的试验数据点(xi,yi)(i=1,2,…,N)，可以构造m次多项式：由曲线拟合的定义，应该使得下式取极小值通过简单运算可以得出系数是下面线性方程组的解