1、2019/5/26,数学建模,MATLAB使用之二 插值问题与拟合问题,应用数学学院 高崇山,2019/5/26,数学建模,插值与拟合的关系,插值和拟合都是函数逼近或者数值逼近的重要组成部分： 他们的共同点都是通过已知一些离散点集M上的约束，求取一个定义 在连续集合S(M包含于S)的未知连续函数，从而达到获取整体规律的目的，即通过“窥几斑”来达到“知全豹”。 简单的讲，所谓拟合是指已知某函数的若干离散函数值f1,f2,fn，通过调整该函数中若干待定系数f(1, 2,n), 使得该函数与已知点集的 差别(最小二乘意义)最小。如果待定函数是线性，就叫线性拟合或者线性回归(主要在统计中)，否则叫作非
2、线性拟合或者非线性回归。表 达式也可以是分段函数，这种情况下叫作样条拟合。 而插值是指已知某函数的在若干离散点上的函数值或者导数信息，通 过求解该函数中待定形式的插值函数以及待定系数，使得该函数在给定离散点上满足约束。插值函数又叫作基函数，如果该基函数定义在整个定义域上，叫作全域基，否则叫作分域基。如果约束条件中只有函数值的约束，叫作Lagrange插值，否则叫作Hermite插值。 从几何意义上将，拟合是给定了空间中的一些点，找到一个已知形式 未知参数的连续曲面来最大限度地逼近这些点；而插值是找到一个( 或几个分片光滑的)连续曲面来穿过这些点。,2019/5/26,数学建模,插值的主要内容,
3、1. 各种插值方法 1.1 Lagrange插值法 1.2 分段插值法 1.3 三次样条插值法 1.4 二维插值 2. 插值的Matlab实现 2.1 一维插值 2.2 二维插值 3. 建模实例：水塔流量的估计,2019/5/26,数学建模,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的多项式是唯一存在的。,2019/5/26,数学建模,Lagrange插值多项式的构造,显然，Ln（x）就是满足插值条件的n次多项式,上式称为基函数,2019/5/26,数学建模,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:mz=x(i);s=0.0;for k=1:np=1.0;for j=1:nif
5、 j=kp=p*(z-x0(j)/(x0(k)-x0(j);endends=p*y0(k)+s;endy(i)=s; end,2019/5/26,数学建模,Lagrange插值法的缺点,多数情况下，Lagrange插值法效果是不错的，但随着节点数n的增大，Lagrange多项式的次数也会升高，可能造成插值函数的收敛性和稳定性变差。如龙格（Runge）现象。 在-1,1上用n+1个等距节点作插值多项式Ln(x),使得它在节点处的值与函数y = 1/(1+25x2)在对应节点的值相等，当n增大时，插值多项式在区间的中间部分趋于y(x)，但对于满足条件0.728|x|1的x， Ln(x)并不趋于y(
6、x)在对应点的值，产生了Runge现象。,2019/5/26,数学建模,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/(n-1):1;y0=1
7、./(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曲线,2019/5/26,数学建模,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),gtext(L12(x),F
8、ontSize,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曲线,2019/5/26,数学建模,1.2 分段线性插值,可以证明：In(x) f(x),2019/5/26,数学建模,1.3 三次样条,设在区间a,b上，已给n+1个互不相同的节点a=x0x1xn=b 而函数y = f(x)在这
9、些节点的值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上有连续的二阶导数。,2019/5/26,数学建模,三次样条,即 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连续 (n-1个方程） Si/
10、(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个线性方程组有唯一解。,2019/5/26,数学建模,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）分片线性插值 （3）双线性插值 1
11、.4.2 散乱点插值法 在T=a,b c,d上散乱分布n个点。一般采用反距离加权平均法。,2019/5/26,数学建模,2.插值的Matlab实现,2.1 一维插值 2.2 二维插值,2019/5/26,数学建模,2.1 一维插值的实现,基本格式： yc=interp1(x,y,cx,method) % x,y分别表示已知数据点的横、纵坐标向量，x必须单调； % cx为需要插值的横坐标数据 % method为插值方法，有nearest 最临近点插值 linear 线性插值（默认）spline 三次样条插值cubic 三次插值 注：Lagrange插值法需自编程序，见前面。,2019/5/26,
12、数学建模,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个元素相关。如果没有对x,y赋值，则默认x=1:n,
13、 y=1:m。n和m分别是矩阵z的行数和列数。,2019/5/26,数学建模,二维插值的实现,2.2.2 插值点为散乱节点 格式：cz=griddata(x,y,z,cx,cy,method) Method有： linear 线性插值 (默认)bilinear 双线性插值cubic 三次插值bicubic 双三次插值nearest 最近邻域插值 注意：cy必须为列向量,2019/5/26,数学建模,3.应用实例,3.1 数控机床加工零件 3.2 山区地形地貌图 3.3 海底曲面图,2019/5/26,数学建模,3.1 数控机床加工零件,图1 零件的轮廓线（x间隔0.2）,表1 x间隔0.2的加
14、工坐标x,y（图1右半部的数据）,加工时需要x每改变0.05时的y值,模型 将图1逆时针方向转90度，轮廓线上下对称，只需对上半部计算一个函数在插值点的值。,图2 逆时针方向转90度的结果,2019/5/26,数学建模,数控机床加工零件 程序,% 按照表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 -1.4 -1.96 -2.37 -2.71
15、.-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）坐标系输出结果 v1 v2 -u subplot(1,
16、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),2019/5/26,数学建模,数控机床加工零件 运行结果,2019/5/26,数学建模,3.2 山区地形地貌图,已知某处山区地形选点测量坐标数据为： x=0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 y
17、=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 8292 96 98 99 95 91 89 86 84 82 8496 98 95 92 90 88 85 84 83 81 8580 81 82 89 95 96 93 92 89 86 8682 85 87 98 99 96 97 88 85 82 8382 85 89 94 95 93 92 91 86 84 8888 92 93 94 95 89 87 86 83 81 9292 96 97 98 96 93 95 84 82
18、 81 8485 85 81 82 80 80 81 85 90 93 9584 86 81 98 99 98 97 96 95 84 8780 81 85 82 83 84 87 90 95 86 8880 82 81 84 85 86 83 82 81 80 8287 88 89 98 99 97 96 98 94 92 87,2019/5/26,数学建模,山区地形地貌图 程序,原始地貌图程序： x=0:.5:5; y=0:.5:6; xx,yy=meshgrid(x,y); z=89 90 87 85 92 91 96 93 90 87 8292 96 98 99 95 91 89 86
19、 84 82 8496 98 95 92 90 88 85 84 83 81 8580 81 82 89 95 96 93 92 89 86 8682 85 87 98 99 96 97 88 85 82 8382 85 89 94 95 93 92 91 86 84 8888 92 93 94 95 89 87 86 83 81 9292 96 97 98 96 93 95 84 82 81 8485 85 81 82 80 80 81 85 90 93 9584 86 81 98 99 98 97 96 95 84 8780 81 85 82 83 84 87 90 95 86 8880
20、82 81 84 85 86 83 82 81 80 8287 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 8292 96 98 99 95 91 89 86 84 82 8496 98 95 92 90 88 85 84 83 81 8580 81 82 89 95 96 93 92 89 86 8682 85 87 98 99 96 97 88 85 82 8382 85 89 94 95 93 92 91 86 84 8888
21、 92 93 94 95 89 87 86 83 81 9292 96 97 98 96 93 95 84 82 81 8485 85 81 82 80 80 81 85 90 93 9584 86 81 98 99 98 97 96 95 84 8780 81 85 82 83 84 87 90 95 86 8880 82 81 84 85 86 83 82 81 80 8287 88 89 98 99 97 96 98 94 92 87; xi=linspace(0,5,50); %加密横坐标数据到50个 yi=linspace(0,6,80); %加密纵坐标数据到60个 xii,yii=
22、meshgrid(xi,yi); %生成网格数据 zii=interp2(x,y,z,xii,yii,cubic); %插值 mesh(xii,yii,zii) %加密后的地貌图,2019/5/26,数学建模,山区地形地貌图 结果,2019/5/26,数学建模,3.3 海底曲面图,例：在某海域测得一些点（x,y)处的水深z由下表给出，在矩形区域(75,200) (-50,150)内画出海底曲面的图形.,2019/5/26,数学建模,海底曲面图 程序,clc;clf;clear all; x=129 140 103.5 88 185.5 195 105 157.5 107.5 77 81 162
23、 162 117.5; y=7.5 141.5 23 147 22.5 137.5 85.5 -6.5 -81 3 56.5 -66.5 84 -33.5; z=- 4 8 6 8 6 8 8 9 9 8 8 9 4 9; 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),2019/5/26,数学建模,海底曲面图 结果,2019/5/26,数学建模,曲线拟合问题,2019/5/26,数学建模,各种拟合方法
24、,1.线性拟合函数 2.多项式曲线拟合函数 3. 稳健回归函数 4. 向自定义函数拟合 5.非线性曲线拟合,2019/5/26,数学建模,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返回的95%的置信区间。 r中为形状残差， rint中返回每一
25、个残差的95%置信区间。 Stats向量包含R2统计量、回归的F值和p值。,2019/5/26,数学建模,线性拟合,例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,2019/5/26,数学建模,2.多项式曲线拟合函数,调用格式： p=polyfit(x,y,n)p,s= polyfit(x,y,n) 说明：x,y为数据点，n为多项式阶
26、数，返回p为幂次从高到低的多项式系数向量p。矩阵s用于生成预测值的误差估计。,2019/5/26,数学建模,多项式曲线拟合函数,例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 2n=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)legend(原始数据,3阶曲线),2019/5/26,数学建模,多项式
27、曲线拟合函数,也可由函数给出数据。 例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阶曲线),2019/5/26,数学建模,多项式曲线拟合函数,再用10阶多项式拟合程序：x=1:20; y=x+3*sin(x); p=polyfit(x,y,10) xi=linspace(1,20,100); z=polyval(p,xi); plot(x,y,o,xi
28、,z,k:,x,y,b) legend(原始数据,10阶多项式),2019/5/26,数学建模,多项式曲线求值函数： 调用格式： 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%的预测值。 多项式曲线拟合的评价和置信区间函数 调用格式： Y,DELTA=polyconf(p,x,s)Y
29、,DELTA=polyconf(p,x,s,alpha) 说明：Y,DELTA=polyconf(p,x,s)使用polyfit函数的选项输出s给出Y的95%置信区间Y DELTA。它假设polyfit函数数据输入的误差是独立正态的，并且方差为常数。1-alpha为置信度。,2019/5/26,数学建模,3.稳健回归函数,稳健回归是指此回归方法相对于其他回归方法而言，受异常值的影响较小。 调用格式： b=robustfit(x,y)b,stats=robustfit(x,y) 说明：b返回系数估计向量；stats返回各种参数估计。,2019/5/26,数学建模,稳健回归函数,例5：演示一个异常
30、数据点如何影响最小二乘拟合值与稳健拟合。首先利用函数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)+brob(2)*x,r),2019/5/26,数学建模,稳健回归函数,分析：稳健拟合(实线)对数
31、据的拟合程度好些，忽略了异常值。最小二乘拟合(点线)则受到异常值的影响，向异常值偏移。,2019/5/26,数学建模,4. 自定义函数拟合,对于给定的数据，根据经验拟合为带有待定常数的自定义函数。 所用函数：nlinfit( ) 调用格式： beta,r,J=nlinfit(X,y,fun,betao) 说明：beta返回函数fun中的待定常数；r表示残差；J表示雅可比矩阵。X,y为数据；fun自定义函数；beta0待定常数初值。,2019/5/26,数学建模,向自定义函数拟合,在化工生产中获得的氯气的级分y随生产时间x下降，假定在x8时，y与x之间有如下形式的非线性模型：现收集了44组数据，
32、利用该数据通过拟合确定非线性模型中的待定常数。 x y x y x y8 0.49 16 0.43 28 0.418 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 12 0.43 24 0.42 36 0.36 14 0.45 24 0.40 38 0.40 14 0
33、.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,2019/5/26,数学建模,自定义函数拟合,首先定义非线性函数的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
34、 18.00 18.00 20.00 20.00 20.00 20.00 22.00 22.00 24.00 24.00 24.00 26.00 26.00 26.00 28.00 28.00 30.00 30.00 30.00 32.00 32.00 34.00 36.00 36.00 38.00 38.00 40.00 42.00;y=0.49 0.49 0.48 0.47 0.48 0.47 0.46 0.46 0.45 0.43 0.45 0.43 0.43 0.44 0.43 0.43 0.46 0.42 0.42 0.43 0.41 0.41 0.40 0.42 0.40 0.40
35、 0.41 0.40 0.41 0.41 0.40 0.40 0.40 0.38 0.41 0.40 0.40 0.41 0.38 0.40 0.40 0.39 0.39;beta0=0.30 0.02; betafit = nlinfit(x,y,model,beta0) 结论：betafit =0.38960.1011 方程为： yy=0.3896+(0.49-0.3896)*exp(-0.1011*(x-8);,2019/5/26,数学建模,5.非线性曲线拟合,LSQCURVEFIT Solves non-linear least squares problems. 功能：根据输入数据x
36、data和得到的输出数据ydata，找到与方程F(x,xdata)最佳的拟合系数。 数学模型：,调用命令：x,resnorm,residual,exitflag,output,lambda,jacobian=lsqcurvefit(fun,x0,xdata,ydata,lb,ub,options,p1,p2,),2019/5/26,数学建模,举例,例、拟合函数： y(i)=a(1)*x(i)2+a(2)*sin(x(i)+a(3)*x(i)3 求解： function F = myfun5(a,x) F = a(1)*x.2 + a(2)*sin(x) + a(3)*x.3; % 假设通过试验得到数据x和y x = 3.6 7.7 9.3 4.1 8.6 2.8 1.3 7.9 10.0 5.4; y = 16.5 150.6 263.1 24.7 208.5 9.9 2.7 163.9 325.0 54.3; x0 = 10, 10, 10 %初值 a,resnorm = lsqcurvefit(myfun5,x0,x,y) xx=1:10; yy= = a(1)*x.2 + a(2)*sin(x) + a(3)*x.3; plot(x,y,o,xx,yy),
