数学建模讲座之五---插值和拟合.ppt

1、2020/9/14,数学建模,MATLAB使用之二插值问题与拟合问题,应用数学学院 高崇山,2020/9/14,数学建模,插值与拟合的关系,插值和拟合都是函数逼近或者数值逼近的重要组成部分： 他们的共同点都是通过已知一些离散点集M上的约束，求取一个定义 在连续集合S(M包含于S)的未知连续函数，从而达到获取整体规律的目的，即通过“窥几斑”来达到“知全豹”。 简单的讲，所谓拟合是指已知某函数的若干离散函数值f1,f2,fn，通过调整该函数中若干待定系数f(1, 2,n), 使得该函数与已知点集的 差别(最小二乘意义)最小。如果待定函数是线性，就叫线性拟合或者线性回归(主要在统计中)，否则叫作非线

2、性拟合或者非线性回归。表 达式也可以是分段函数，这种情况下叫作样条拟合。 而插值是指已知某函数的在若干离散点上的函数值或者导数信息，通 过求解该函数中待定形式的插值函数以及待定系数，使得该函数在给定离散点上满足约束。插值函数又叫作基函数，如果该基函数定义在整个定义域上，叫作全域基，否则叫作分域基。如果约束条件中只有函数值的约束，叫作Lagrange插值，否则叫作Hermite插值。 从几何意义上将，拟合是给定了空间中的一些点，找到一个已知形式 未知参数的连续曲面来最大限度地逼近这些点；而插值是找到一个( 或几个分片光滑的)连续曲面来穿过这些点。,2020/9/14,数学建模,插值的主要内容,1

3、. 各种插值方法 1.1 Lagrange插值法 1.2 分段插值法 1.3 三次样条插值法 1.4 二维插值 2. 插值的Matlab实现 2.1 一维插值 2.2 二维插值 3. 建模实例：水塔流量的估计,2020/9/14,数学建模,1.1 Lagrange插值,已知y = f(x)（该函数未知）在互异的n+1个点x0,x1,x2,xn处的函数值y0, y1,y2,yn,则构造一个过n+1个点（xk,yk）k=0,1,2,n的次数不超过n的多项式 y = Ln(x)， （称为插值多项式） 使其满足 Ln(xk) = yk （称为插值条件） 然后用y = Ln(x)作为准确函数y = f(

4、x)的近似值。此方法称为插值法。 Theorem：满足插值条件的次数不超过n的多项式是唯一存在的。,2020/9/14,数学建模,Lagrange插值多项式的构造,显然，Ln（x）就是满足插值条件的n次多项式,上式称为基函数,2020/9/14,数学建模,Lagrange插值程序,function y=lagr(x0,y0,x) % lagrange插值法的程序 （x0，y0)表示已知的n个节点 % x表示m个插值点 y表示对应于x的m个插值 n=length(x0); m=length(x); for i=1:m z=x(i); s=0.0; for k=1:n p=1.0; for j=1

5、:n if j=k p=p*(z-x0(j)/(x0(k)-x0(j); end end s=p*y0(k)+s; end y(i)=s; end,2020/9/14,数学建模,Lagrange插值法的缺点,多数情况下，Lagrange插值法效果是不错的，但随着节点数n的增大，Lagrange多项式的次数也会升高，可能造成插值函数的收敛性和稳定性变差。如龙格（Runge）现象。 在-1,1上用n+1个等距节点作插值多项式Ln(x),使得它在节点处的值与函数y = 1/(1+25x2)在对应节点的值相等，当n增大时，插值多项式在区间的中间部分趋于y(x)，但对于满足条件0.728|x|1的x，

6、Ln(x)并不趋于y(x)在对应点的值，产生了Runge现象。,2020/9/14,数学建模,Runge现象的程序(1),clc;clf;clear all; m=21; x=-1:1/(m-1):1; y=1./(1+25*x.2);z=0*x; n=3; x0=-1:1/(n-1):1;y0=1./(1+25*x0.2);y1=lagr(x0,y0,x); subplot(2,2,1), plot(x,z,r-,x,y,m-),hold on %原曲线 plot(x,y1,b),gtext(L4(x),FontSize,12),pause %Lagrange曲线 n=5; x0=-1:1/

7、(n-1):1;y0=1./(1+25*x0.2);y1=lagr(x0,y0,x); subplot(2,2,2), plot(x,z,r-,x,y,m-),hold on %原曲线 plot(x,y1,b),gtext(L8(x),FontSize,12),pause %Lagrange曲线,2020/9/14,数学建模,Runge现象的程序(2),n=7; x0=-1:1/(n-1):1;y0=1./(1+25*x0.2);y1=lagr(x0,y0,x); subplot(2,2,3), plot(x,z,r-,x,y,m-),hold on, %原曲线 plot(x,y1,b),gt

8、ext(L12(x),FontSize,12),pause %Lagrange曲线 n=9; x0=-1:1/(n-1):1;y0=1./(1+25*x0.2);y1=lagr(x0,y0,x); subplot(2,2,4), plot(x,z,r-,x,y,m-),hold on, %原曲线 plot(x,y1,b),gtext(L16(x),FontSize,12) %Lagrange曲线,2020/9/14,数学建模,1.2 分段线性插值,可以证明：In(x) f(x),2020/9/14,数学建模,1.3 三次样条,设在区间a,b上，已给n+1个互不相同的节点 a=x0x1xn=b

9、而函数y = f(x)在这些节点的值f(xi)=yi,i=0,1,n.如果分段函数S(x)满足下列条件，就称S(x)为f(x)在点x0，x1，xn的三次样条插值函数. （1） S(x)在子区间xi,xi+1的表达式Si(x)都是次数为3的多项式； （2）S(xi) = yi; （3） S(x)在区间a,b上有连续的二阶导数。,2020/9/14,数学建模,三次样条,即 Si(x)=aix3+bix2+cix+di i=0,1,n xi-1x xi （4n个变量） 需要4n个方程 S(xi) = yi i=0,1,n (n+1个方程） Si(xi)= Si+1(xi) i=1,n-1 在xi连续

10、 (n-1个方程） Si/(xi)= Si+1/(xi) i=1,n-1 在xi连续(n-1个方程） Si/(xi)= Si+1 /(xi) i=1,n-1 在xi连续(n-1个方程） 再加两个条件 S/(x0)= S /(xn)=0 自然边界条件（2个方程） 可以证明：满足上述4n个线性方程组有唯一解。,2020/9/14,数学建模,1.4 二维插值,1.4.1 网格节点插值法 已知mn个节点（xi,yj,zij）(i=1,2,m,j=1,2,n），一般设a=x1x2,xm=b,c= y1y2yn=d,求任意一点（x*，y*）（ （xi,yj）处的插值z*. （1）最临近点插值 （2）分片线

11、性插值 （3）双线性插值 1.4.2 散乱点插值法 在T=a,b c,d上散乱分布n个点。一般采用反距离加权平均法。,2020/9/14,数学建模,2.插值的Matlab实现,2.1 一维插值 2.2 二维插值,2020/9/14,数学建模,2.1 一维插值的实现,基本格式： yc=interp1(x,y,cx,method) % x,y分别表示已知数据点的横、纵坐标向量，x必须单调； % cx为需要插值的横坐标数据 % method为插值方法，有 nearest 最临近点插值 linear 线性插值（默认） spline 三次样条插值 cubic 三次插值 注：Lagrange插值法需自编程

12、序，见前面。,2020/9/14,数学建模,2.2 二维插值的实现,2.2.1 插值节点为网格节点，即x,y向量是单调的。 调用格式：zi=interp2(x,y,z,xi,yi,method) Method的4种情况： nearest 最临近点插值 linear 线性插值（默认） spline 三次样条插值 cubic 三次插值 说明：这里x和y是两个独立的向量，它们必须是单调的。z是矩阵，是由x和y确定的点上的值。z和x,y之间的关系是z(i,:)=f(x,y(i) z(:,j)=f(x(j),y) 即：当x变化时，z的第i行与y的第i个元素相关，当y变化时z的第j列与x的第j个元素相关。

13、如果没有对x,y赋值，则默认x=1:n, y=1:m。n和m分别是矩阵z的行数和列数。,2020/9/14,数学建模,二维插值的实现,2.2.2 插值点为散乱节点 格式：cz=griddata(x,y,z,cx,cy,method) Method有： linear线性插值 (默认) bilinear 双线性插值 cubic 三次插值 bicubic双三次插值 nearest最近邻域插值 注意：cy必须为列向量,2020/9/14,数学建模,3.应用实例,3.1 数控机床加工零件 3.2 山区地形地貌图 3.3 海底曲面图,2020/9/14,数学建模,3.1 数控机床加工零件,图1 零件的轮廓

14、线（x间隔0.2）,表1 x间隔0.2的加工坐标x,y（图1右半部的数据）,加工时需要x每改变0.05时的y值,模型 将图1逆时针方向转90度，轮廓线上下对称，只需对上半部计算一个函数在插值点的值。,图2 逆时针方向转90度的结果,2020/9/14,数学建模,数控机床加工零件 程序,% 按照表1输入原始数据 x=0:0.2:5, 4.8:-0.2:0; y=5 4.71 4.31 3.68 3.05 2.5 2.05 1.69 1.4 1.18 1 0.86 0.74 0.64 0.57 0.5 . 0.44 0.4 0.36 0.32 0.29 0.26 0.24 0.2 0.15 0 -

15、1.4 -1.96 -2.37 -2.71 . -3 -3.25 -3.47 -3.67 -3.84 -4 -4.14 -4.27 -4.39 -4.49 -4.58 -4.66 . -4.74 -4.8 -4.85 -4.9 -4.94 -4.96 -4.98 -4.99 -5; % 逆时针方向转90度，节点（x, y）变为（u, v） v0=x; u0=-y; % 按0.05的间隔在u方向产生插值点 u=-5:0.05:5; % 在v方向计算分段线性插值 v1=interp1(u0,v0,u); % 在v方向计算三次样条插值 v2=spline(u0,v0,u); % 在（x, y）坐标系

16、输出结果 v1 v2 -u subplot(1,3,1),plot(x,y),axis(0 5 -5 5) gtext(原轮廓线,FontSize,12) subplot(1,3,2),plot(v1,-u),axis(0 5 -5 5) gtext(分段线性插值,FontSize,12) subplot(1,3,3),plot(v2,-u),axis(0 5 -5 5) gtext(三次样条插值,FontSize,12),2020/9/14,数学建模,数控机床加工零件 运行结果,2020/9/14,数学建模,3.2 山区地形地貌图,已知某处山区地形选点测量坐标数据为： x=0 0.5 1 1

17、.5 2 2.5 3 3.5 4 4.5 5 y=0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 5.5 6 海拔高度数据为： z=89 90 87 85 92 91 96 93 90 87 82 92 96 98 99 95 91 89 86 84 82 84 96 98 95 92 90 88 85 84 83 81 85 80 81 82 89 95 96 93 92 89 86 86 82 85 87 98 99 96 97 88 85 82 83 82 85 89 94 95 93 92 91 86 84 88 88 92 93 94 95 89 87 86 83 8

18、1 92 92 96 97 98 96 93 95 84 82 81 84 85 85 81 82 80 80 81 85 90 93 95 84 86 81 98 99 98 97 96 95 84 87 80 81 85 82 83 84 87 90 95 86 88 80 82 81 84 85 86 83 82 81 80 82 87 88 89 98 99 97 96 98 94 92 87,2020/9/14,数学建模,山区地形地貌图 程序,原始地貌图程序： x=0:.5:5; y=0:.5:6; xx,yy=meshgrid(x,y); z=89 90 87 85 92 91 9

19、6 93 90 87 82 92 96 98 99 95 91 89 86 84 82 84 96 98 95 92 90 88 85 84 83 81 85 80 81 82 89 95 96 93 92 89 86 86 82 85 87 98 99 96 97 88 85 82 83 82 85 89 94 95 93 92 91 86 84 88 88 92 93 94 95 89 87 86 83 81 92 92 96 97 98 96 93 95 84 82 81 84 85 85 81 82 80 80 81 85 90 93 95 84 86 81 98 99 98 97 9

20、6 95 84 87 80 81 85 82 83 84 87 90 95 86 88 80 82 81 84 85 86 83 82 81 80 82 87 88 89 98 99 97 96 98 94 92 87; mesh(xx,yy,z),加密后的地貌图 x=0:.5:5;y=0:.5:6; z=89 90 87 85 92 91 96 93 90 87 82 92 96 98 99 95 91 89 86 84 82 84 96 98 95 92 90 88 85 84 83 81 85 80 81 82 89 95 96 93 92 89 86 86 82 85 87 98 99

21、 96 97 88 85 82 83 82 85 89 94 95 93 92 91 86 84 88 88 92 93 94 95 89 87 86 83 81 92 92 96 97 98 96 93 95 84 82 81 84 85 85 81 82 80 80 81 85 90 93 95 84 86 81 98 99 98 97 96 95 84 87 80 81 85 82 83 84 87 90 95 86 88 80 82 81 84 85 86 83 82 81 80 82 87 88 89 98 99 97 96 98 94 92 87; xi=linspace(0,5,

22、50); %加密横坐标数据到50个 yi=linspace(0,6,80); %加密纵坐标数据到60个 xii,yii=meshgrid(xi,yi); %生成网格数据 zii=interp2(x,y,z,xii,yii,cubic); %插值 mesh(xii,yii,zii) %加密后的地貌图,2020/9/14,数学建模,山区地形地貌图 结果,2020/9/14,数学建模,3.3 海底曲面图,例：在某海域测得一些点（x,y)处的水深z由下表给出，在矩形区域(75,200) (-50,150)内画出海底曲面的图形.,2020/9/14,数学建模,海底曲面图 程序,clc;clf;clear

23、 all; x=129140103.5 88185.5195105157.5107.57781162162117.5; y=7.5141.523147 22.5137.585.5-6.5-81 3 56.5-66.584-33.5; z=- 4868688 9 988949; plot3(x,y,z,o), hold on %原始数据点 % 插值 cx=75:0.5:200; cy=-70:0.5:150; cz=griddata(x,y,z,cx,cy,cubic);% 三次插值 meshz(cx,cy,cz),2020/9/14,数学建模,海底曲面图 结果,2020/9/14,数学建模,曲

24、线拟合问题,2020/9/14,数学建模,各种拟合方法,1.线性拟合函数 2.多项式曲线拟合函数 3. 稳健回归函数 4.向自定义函数拟合 5.非线性曲线拟合,2020/9/14,数学建模,1.线性拟合,调用格式： b=regress(y,X) b,bint,r,rint,stats= regress(y,X) b,bint,r,rint,stats= regress(y,X,alpha) 说明： b=regress(y,X)返回X处y的最小二乘拟合值。该函数求解线性模型：y=X+ 是p1的参数向量； 是服从标准正态分布的随机干扰的n1的向量； y为n1的向量； X为np矩阵。 bint返回的

25、95%的置信区间。 r中为形状残差， rint中返回每一个残差的95%置信区间。 Stats向量包含R2统计量、回归的F值和p值。,2020/9/14,数学建模,线性拟合,例1：设y的值为给定的x的线性函数加服从标准正态分布的随机干扰值得到。即y=10+x+. 求线性拟合方程系数。 程序： x=ones(10,1) (1:10) y=x*10;1+normrnd(0,0.1,10,1) b,bint=regress(y,x,0.05) 回归方程为：y=9.9213+1.0143x,2020/9/14,数学建模,2.多项式曲线拟合函数,调用格式： p=polyfit(x,y,n) p,s= po

26、lyfit(x,y,n) 说明：x,y为数据点，n为多项式阶数，返回p为幂次从高到低的多项式系数向量p。矩阵s用于生成预测值的误差估计。,2020/9/14,数学建模,多项式曲线拟合函数,例2：由离散数据 x0.1.2.3.4.5.6.7.8.91y.3.511.41.61.9.6.4.81.52拟合出多项式。 程序： x=0:.1:1; y=.3 .5 1 1.4 1.6 1.9 .6 .4 .8 1.5 2 n=3; p=polyfit(x,y,n) xi=linspace(0,1,100); z=polyval(p,xi); %多项式求值 plot(x,y,o,xi,z,k:,x,y,b

27、) legend(原始数据,3阶曲线),2020/9/14,数学建模,多项式曲线拟合函数,也可由函数给出数据。 例3：x=1:20,y=x+3*sin(x) 程序： x=1:20; y=x+3*sin(x); p=polyfit(x,y,6) xi=1inspace(1,20,100); z=poyval(p,xi); %多项式求值函数 plot(x,y,o,xi,z,k:,x,y,b) legend(原始数据,6阶曲线),2020/9/14,数学建模,多项式曲线拟合函数,再用10阶多项式拟合 程序：x=1:20; y=x+3*sin(x); p=polyfit(x,y,10) xi=lins

28、pace(1,20,100); z=polyval(p,xi); plot(x,y,o,xi,z,k:,x,y,b) legend(原始数据,10阶多项式),2020/9/14,数学建模,多项式曲线求值函数： 调用格式： y=polyval(p,x) y,DELTA=polyval(p,x,s) 说明：y=polyval(p,x)为返回对应自变量x在给定系数p的多项式的值。 y,DELTA=polyval(p,x,s) 使用polyfit函数的选项输出s得出误差估计Y DELTA。它假设polyfit函数数据输入的误差是独立正态的，并且方差为常数。则Y DELTA将至少包含50%的预测值。 多

29、项式曲线拟合的评价和置信区间函数 调用格式： Y,DELTA=polyconf(p,x,s) Y,DELTA=polyconf(p,x,s,alpha) 说明：Y,DELTA=polyconf(p,x,s)使用polyfit函数的选项输出s给出Y的95%置信区间Y DELTA。它假设polyfit函数数据输入的误差是独立正态的，并且方差为常数。1-alpha为置信度。,2020/9/14,数学建模,3.稳健回归函数,稳健回归是指此回归方法相对于其他回归方法而言，受异常值的影响较小。 调用格式： b=robustfit(x,y) b,stats=robustfit(x,y) 说明：b返回系数估计

30、向量；stats返回各种参数估计。,2020/9/14,数学建模,稳健回归函数,例5：演示一个异常数据点如何影响最小二乘拟合值与稳健拟合。首先利用函数y=10-2x加上一些随机干扰的项生成数据集，然后改变一个y的值形成异常值。调用不同的拟合函数，通过图形观查影响程度。 程序：x=(1:10); y=10-2*x+randn(10,1); y(10)=0; bls=regress(y,ones(10,1) x) %线性拟合 brob=robustfit(x,y) %稳健拟合 scatter(x,y) hold on plot(x,bls(1)+bls(2)*x,:) plot(x,brob(1)

31、+brob(2)*x,r),2020/9/14,数学建模,稳健回归函数,分析：稳健拟合(实线)对数据的拟合程度好些，忽略了异常值。最小二乘拟合(点线)则受到异常值的影响，向异常值偏移。,2020/9/14,数学建模,4. 自定义函数拟合,对于给定的数据，根据经验拟合为带有待定常数的自定义函数。 所用函数：nlinfit( ) 调用格式： beta,r,J=nlinfit(X,y,fun,betao) 说明：beta返回函数fun中的待定常数；r表示残差；J表示雅可比矩阵。X,y为数据；fun自定义函数；beta0待定常数初值。,2020/9/14,数学建模,向自定义函数拟合,在化工生产中获得的

32、氯气的级分y随生产时间x下降，假定在x8时，y与x之间有如下形式的非线性模型： 现收集了44组数据，利用该数据通过拟合确定非线性模型中的待定常数。 x y x y x y 8 0.49 16 0.43 28 0.41 8 0.49 18 0.46 28 0.40 10 0.48 18 0.45 30 0.40 10 0.47 20 0.42 30 0.40 10 0.48 20 0.42 30 0.38 10 0.47 20 0.43 32 0.41 12 0.46 20 0.41 32 0.40 12 0.46 22 0.41 34 0.40 12 0.45 22 0.40 36 0.41

33、12 0.43 24 0.42 36 0.36 14 0.45 24 0.40 38 0.40 14 0.43 24 0.40 38 0.40 14 0.43 26 0.41 40 0.36 16 0.44 26 0.40 42 0.39 16 0.43 26 0.41,2020/9/14,数学建模,自定义函数拟合,首先定义非线性函数的m文件：model.m function yy=model(beta0,x) a=beta0(1); b=beta0(2); yy=a+(0.49-a)*exp(-b*(x-8); 程序： x=8.00 8.00 10.00 10.00 10.00 10.00 12.00 12.00 12.00 14.00 14.00 14.00 16.00 16.00 16.00 18.00 18.00 20.00 20.00 20.00 20.00