Matlab数据插值与拟合

1、数据插值与拟合在工程实践与科学实验中,常常需要从一组试验数据之中找到自变量与因变量之间的关系，一般可用一个近似函数表示。函数产生的办法因观测数据的要求不同而异，数据插值与拟合是两种常用的方法。4.1 MATLAB中的插值函数4.2 拉格朗日插值法4.3 利用均差的牛顿插值法4.4 利用差分的牛顿插值法4.5 Hermite插值4.6 spline三次样条插值4.7 多项式曲线拟合4.8 最小二乘拟合1精选PPT4.1 MATLAB中的插值函数函数插值来源于函数的以下问题：只知道函数在某区间有定义且已得到区间内一些离散点的值，希望用简单的表达式近似给出函数在此区间上的整体描述，并能与已知离散点上

2、的值相等。插值法按插值函数的形式主要分为以下几种形式： （1）代数多项式插值； （2）三角多项式插值； （3）有理分式插值。 2精选PPT代数多项式插值是最常用的插值方式，其内容也是最丰富的，它又可分为以下几种插值方式： （1）非等距节点插值，包括拉格朗日插值、利用均差的牛顿插值和埃特金插值； （2）非等距节点插值，包括利用差分的牛顿插值和高斯插值等； （3）在插值中增加了导数的Hermite（埃尔米特）插值； （4）分段插值，包括分段线性插值、分段Hermite（埃尔米特）插值和样条函数插值；（5）反插值。按被插值函数的变量个数还可把插值法分为一元插值和多元插值。3精选PPT4.1.1 一元

3、插值函数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

4、(x,Y,xi,method,extrap)对xi中超出已知点集的插值点用指定插值方法计算函数值y=interp1(x,Y,xi,method,extrap,extrapval)用指定方法插值xi上的函数值，超出已知点集处函数值取extrapvaly=interp1(x,Y,xi,method,pp)用指定方法插值，但返回结果为分段多项式4精选PPTmethod方法描述nearest最邻近插值：插值点处函数值取与插值点最邻近的已知点的函数值liner分段线性插值：插值点处函数值由连接其最邻近的两侧点的线性函数预测，MATLAB中interp1的默认方法spline样条插值：默认为三次样条插值。

5、可用spline函数代替pchip三次Hermite多项式插值。可用pchip函数代替cubic同pchip，三次Hermite多项式插值MATLAB中一维插值有多种算法，由interp1函数中的method指定。 MATLAB中一维插值的各种算法如表4-2所示。表4-2 一维插值算法(method)5精选PPT1.Linear（分段线性插值）它的算法是在每个小区间xi,xi+1上采用简单的线性插值。在区间xi,xi+1上的子插值多项式为： 由此整个区间xi,xi+1上的插值函数为： 其中 定义如下：6精选PPT分段线性插值方法在速度和误差之间取得了比较好的均衡，其插值函数具有连续性，但在已知

6、数据点处的斜率一般不会改变，因此不是光滑的。分段线性插值方法是MATLAB一维插值默认的方法。7精选PPT2.Spline（样条插值）样条插值是用分段低次多项式去逼近函数。样条函数可以给出光滑 的插值曲线，只要在插值区间端点提供某些导数信息，样条插值可以适应不同光滑需求。三次样条是使用最为广泛的样条插值，它在每个子区间xi,xi+1上都是有二阶连续导数的三次多项式，即 其中 都是三次多项式。 8精选PPT对于给定的离散的测量数据经x，y（称为断点），要寻找一个三次多项式y=p(x)， 以逼近每对数据(xi,yi)点间曲线。过两点(xi,yi)和(xi+1,yi+1)只能确定一条直线，而通过一点

7、的三次多项式曲线有无穷多条。为使通过中间断点的三次多项式曲线具有唯一性，要增加以下的连续条件和边界条件（因为三次多项式有4个系数）： （1）三次多项式在点(xi,yi)处有： ； （2）三次多项式在点(xi,yi)处有： ； （3）三次多项式在点(xi,yi)处有： ； （4）边界条件： 。9精选PPT表4-2中各种方法中： （1）nearest方法速度最快，占用内存最小，但一般来说误差最大，插值结果最不光滑； （2）spline三次样条插值是所有插值方法中运行耗时最长的，其插值函数以及插值函数的一阶、二阶导函数都连续，因此是最光滑的插值方法，占用内存上比cubic方法小，但当已知数据点不均匀

8、分布时可能出现异常结果。 （3）cubic三次多项式插值法中插值函数及其一阶导数都是连续的，因此其插值结果也比较光滑，运算速度比spline方法略快，但占用内存最多。在实际的使用中，应根据实际需求和运算条件选择合适的算法。10精选PPT例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画出插值结果。从图中

9、可以看出分段线性就是联结两个邻近的已知点的线性函数插值计算该区间内插值点上的函数值。11精选PPT例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);

10、plot(x,y,o,xx,y1,-,xx,y2,-.,xx,y3,:)12精选PPT13精选PPT4.2 拉格朗日插值法拉格朗日插值法是基于基函数的插值方法，插值多项式可表示为 其中 称为i次基函数： 14精选PPT在MATLAB中编程实现拉格朗日插值法函数为：Language。功能：求已知数据点的拉格朗日多项式；调用格式：f= Language(x,y)或f= Language(x,y,x0)。 其中，x为已知数据点的x 坐标向量； y为已知数据点的y 坐标向量； x0为插值点的x坐标； f为求得的拉格朗日多项式或x0处的插值。15精选PPT function f = Language(x

11、,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;

12、 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 endend16精选PPT例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.054

13、27*t-.145485e-1*t2-.204917*t3+.328112e-1*t4-.261189e-1 f=language(x,y,1.6)f = 0.9992x11.21.82. 54 y0.84150.93200.97380.5985-0.756817精选PPT4.3 利用均差的牛顿插值法函数 f的零阶均差定义为 ，一阶定义均差为 一般地，函数f的k阶均差定义为：利用均差的牛顿插值法多项式为18精选PPT系数的计算过程如表4-3所示表4-3 均差计算表格一阶均差二阶均差三阶均差n阶均差19精选PPT在MATLAB中编程实现利用均差牛顿插值法函数为：Newton。功能：求已知数据点的

14、均差形式的牛顿插值多项式；调用格式：f= Newton(x,y)或f= Newton(x,y,x0)。 其中，x为已知数据点的x 坐标向量； y为已知数据点的y 坐标向量； x0为插值点的x坐标； f为求得的牛顿插值法多项式或x0处的插值。20精选PPTfunction 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;e

15、lse disp(x和y的维数不相等！); return;end21精选PPTf = 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 endend22精选PPT例4-

16、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 = 4x11.21.82. 54 y11.443.246.251623精选PPT4.4 利用差分的牛顿插值法差分分为向前差分、后向差分和中心差分三种，它们的记法及定义如下为： 其中： 代表向前差分； 代表向后差分； 代表向后差

17、分。24精选PPT 假设 。为了方便，可构造如表4-4所示的差分表（ ）。 表4-4 差分计算表格25精选PPT4.4.1向前牛顿插值向前牛顿插值多项式可表示如下： 其中 叫做步长， ， 且 的取值范围为 。26精选PPT在MATLAB中编程实现向前牛顿插值法函数为：Newtonforward。功能：求已知数据点的向前牛顿插值法多项式；调用格式：f= Newtonforward(x,y)或 f= Newtonforward (x,y,x0)。 其中，x为已知数据点的x 坐标向量； y为已知数据点的y 坐标向量； x0为插值点的x坐标； f为求得的向前牛顿插值法多项式或x0处的插值。27精选PP

18、Tfunction 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;end28精选PPTf = y(1);y1 = 0;xx =linspace(x(1),x(n),(x(2)-x(1);if(xx = x) disp(节点之间不是等距的！); r

19、eturn;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 endend29精选PPT4.4.2向前牛顿插值向后牛顿插值多项式可表示如下：

20、 其中 叫做步长， ， 且 的取值范围为 。30精选PPT在MATLAB中编程实现向后牛顿插值法函数为：Newtonback。功能：求已知数据点的向前牛顿插值法多项式；调用格式：f= Newtonback (x,y)或 f = Newtonback (x,y,x0)。 其中，x为已知数据点的x 坐标向量； y为已知数据点的y 坐标向量； x0为插值点的x坐标； f为求得的向前牛顿插值法多项式或x0处的插值。31精选PPTfunction f = Newtonback(x,y,x0)%求已知数据点的向后差分牛顿插值多项式%已知数据点的x 坐标向量:x%已知数据点的y 坐标向量:y%为插值点的x坐

21、标:x0%求得的向前差分牛顿插值多项式或x0处的插值:fsyms t;if(length(x) = length(y) n = length(x); c(1:n) = 0.0;else disp(x和y的维数不相等！); return;end32精选PPTf = 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

22、) 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 endend33精选PPT例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=Newtonforwa

23、rd(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.9998x11.21.41.61.8 y0.84150.93200.98540.99960.973834精选PPT4.5 Hermite插值Hermite插值满足在

24、节点上等于给定函数值，而且在节点上的导数值也等于给定的导数值。对于高阶导数的情况， Hermite插值多项式比较复杂，在实际中，常常遇到的是函数值与一阶导数给定的情况。在此情况下，n个节点x1，x2，xn的Hermite插值多项式的表达式如下： 其中 ， ， ，35精选PPT在MATLAB中编程实现Hermite插值法函数为： Hermite。功能：求已知数据点的Hermite插值法多项式；调用格式：f= Hermite (x,y,y_1)或 f = Hermite (x,y, y_1,x0)。 其中，x为已知数据点的x 坐标向量； y为已知数据点的y 坐标向量； y_1为已知数据点导数向量；

25、 x0为插值点的x坐标； f为求得的Hermite插值法多项式或x0处的插值。36精选PPTfunction 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的

26、维数不相等！); return;end37精选PPTfor 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 endend38精选PPT例4-6 根据下表的数据点求出其Hermite插值多项式，并计算当x=1.144时y的值。解：

27、 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.60

28、000)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 x11.21.41.61.8 y11.09541.18321.26491.3416y0.50000.45640.42260.39530.372739精选P

29、PT4.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;40精选PPTg=

30、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)&xx=x(i+1); end xbar=(xx(kk)-x(i)/h(i);

31、s(kk)=alpha0(xbar)*y(i)+alpha1(xbar)*y(i+1)+. h(i)*beta0(xbar)*m(i)+h(i)*beta1(xbar)*m(i+1); end m=s;end41精选PPT%追赶法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;42精选PPT4.7 多项式曲线拟合对给定