数学建模讲座插值和拟合.pptx

单击此处编辑母版标题样式,单击此处编辑母版文本样式,第二级,第三级,第四级,第五级,2021/10/3,#,MATLAB使用之二,插值问题与拟合问题,应用数学学院 高崇山,第1页,插值与拟合关系,插值,和,拟合,都是函数迫近或者数值迫近主要组成部分：他们共同点都是经过已知一些离散点集M上约束，求取一个定义 在连续集合S(M包含于S)未知连续函数，从而到达获取整体规律目标，即经过“窥几斑”来到达“知全豹”。简单讲，所谓拟合是指已知某函数若干离散函数值f,1,f,2,f,n,，经过调整该函数中若干待定系数f(,1,2,n,),使得该函数与已知点集 差异(最小二乘意义)最小。假如待定函数是线性，就叫线性拟合或者线性回归(主要在统计中)，不然叫作非线性拟合或者非线性回归。表 达式也能够是分段函数，这种情况下叫作样条拟合。而插值是指已知某函数在若干离散点上函数值或者导数信息，通 过求解该函数中待定形式插值函数以及待定系数，使得该函数在给定离散点上满足约束。插值函数又叫作基函数，假如该基函数定义在整个定义域上，叫作全域基，不然叫作分域基。假如约束条件中只有函数值约束，叫作Lagrange插值，不然叫作Hermite插值。从几何意义上将，拟合是给定了空间中一些点，找到一个已知形式 未知参数连续曲面来最大程度地迫近这些点；而插值是找到一个(或几个分片光滑)连续曲面来穿过这些点。,第2页,插值主要内容,1.各种插值方法,1.1 Lagrange插值法,1.2 分段插值法,1.3 三次样条插值法,1.4 二维插值,2.插值Matlab实现,2.1 一维插值,2.2 二维插值,3.建模实例：水塔流量预计,第3页,1.1 Lagrange插值,已知,y,=,f,(,x,)（该函数未知）在互异n+1个点,x,0,x,1,x,2,x,n,处函数值,y,0,y,1,y,2,y,n,则结构一个过n+1个点（,x,k,y,k,）,k,=0,1,2,n次数不超出n多项式,y,=,L,n,(,x,)，（称为插值多项式）,使其满足,L,n,(,x,k,)=,y,k,（称为插值条件）,然后用,y,=,L,n,(,x,)作为准确函数,y,=,f,(,x,)近似值。此方法称为,插值法,。,Theorem：满足插值条件次数不超出n多项式是唯一存在。,第4页,Lagrange插值多项式结构,显然，,L,n,（,x,）就是满足插值条件n次多项式,上式称为基函数,第5页,Lagrange插值程序,function y=lagr(x0,y0,x),%lagrange插值法程序 （x0，y0,),表示已知n个节点,%x表示m个插值点 y表示对应于xm个插值,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: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,第6页,Lagrange插值法缺点,多数情况下，Lagrange插值法效果是不错，但伴随节点数n增大，Lagrange多项式次数也会升高，可能造成插值函数收敛性和稳定性变差。如龙格（Runge）现象。,在-1,1上用n+1个等距节点作插值多项式,L,n,(,x,),使得它在节点处值与函数,y,=1/(1+25,x,2,)在对应节点值相等，当n增大时，插值多项式在区间中间部分趋于,y,(,x,)，但对于满足条件0.728|,x,|1,x,，,L,n,(,x,)并不趋于,y,(,x,)在对应点值，产生了,Runge现象。,第7页,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./(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曲线,第8页,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),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曲线,第9页,1.2 分段线性插值,x,j,x,j-1,x,j+1,x,0,x,n,能够证实：,I,n,(,x,),f,(,x,),第10页,1.3 三次样条,设在区间,a,b,上，已给,n,+1个互不相同节点,a,=,x,0,x,1,x,n,=,b,而函数,y=f,(,x,)在这些节点值,f,(,x,i,)=,y,i,i,=0,1,n,.假如分段函数,S,(,x,)满足以下条件，就称,S,(,x,)为,f,(,x,)在点,x,0,，,x,1,，,x,n,三次样条插值函数.,（1）,S,(,x,)在子区间,x,i,x,i+,1,表示式,S,i,(,x,)都是次数为3多项式；,（2）,S,(,x,i,)=,y,i,;,（3）,S,(,x,)在区间,a,b,上有连续二阶导数。,第11页,三次样条,即,S,i,(,x,)=,a,i,x,3,+,b,i,x,2,+,c,i,x,+,d,i,i,=0,1,n x,i,-1,x,x,i,（4n个变量）,需要4,n,个方程,S,(,x,i,)=,y,i,i,=0,1,n,(,n+,1个方程）,S,i,(,x,i,)=,S,i+,1,(,x,i,),i,=1,n-,1 在,x,i,连续(,n-,1个方程）,S,i,/,(,x,i,)=,S,i+,1,/,(,x,i,),i,=1,n-,1 在,x,i,连续(,n-,1个方程）,S,i,/,(,x,i,)=,S,i+,1,/,(,x,i,),i,=1,n-,1 在,x,i,连续(,n-,1个方程）,再加两个条件,S,/,(,x,0,)=,S,/,(,x,n,)=0 自然边界条件（2个方程）,能够证实：,满足上述4n个线性方程组有唯一解,。,第12页,1.4 二维插值,1.4.1 网格节点插值法,已知mn个节点（,x,i,y,j,z,ij,）(,i,=1,2,m,j,=1,2,n,），普通设,a=x,1,x,2,x,m,=b,c=y,1,y,2,y,n,=,d,求任意一点（,x,*,，y,*,）（,（,x,i,y,j,）处插值,z,*,.,（1）最临近点插值,（2）分片线性插值,（3）双线性插值,1.4.2 散乱点插值法,在T=a,b c,d上散乱分布n个点。普通采取反距离加权平均法。,第13页,2.插值Matlab实现,2.1 一维插值,2.2 二维插值,第14页,2.1 一维插值实现,基本格式：,yc=interp1(x,y,cx,method),%x,y分别表示已知数据点横、纵坐标向量，x必须单调；,%cx为需要插值横坐标数据,%method为插值方法，有,nearest 最临近点插值,linear 线性插值（默认）,spline 三次样条插值,cubic 三次插值,注：Lagrange插值法需自编程序，见前面。,第15页,2.2 二维插值实现,2.2.1 插值节点为网格节点，即x,y向量是单调。,调用格式：zi=interp2(x,y,z,xi,yi,method),Method4种情况：,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,y=1:m。n和m分别是矩阵z行数和列数。,第16页,二维插值实现,2.2.2 插值点为散乱节点,格式：cz=griddata(x,y,z,cx,cy,method),Method有：,linear线性插值(默认),bilinear 双线性插值,cubic 三次插值,bicubic双三次插值,nearest最近邻域插值,注意,：cy必须为列向量,第17页,3.应用实例,3.1 数控机床加工零件,3.2 山区地形地貌图,3.3 海底曲面图,第18页,3.1 数控机床加工零件,图1,零件轮廓线（,x,间隔0.2）,表1,x,间隔0.2加工坐标,x,y,（图1右半部数据）,0.0,5.00,0.2,4.71,0.4,4.31,0.6,3.68,0.8,3.05,1.0,2.50,1.2,2.05,1.4,1.69,1.6,1.40,1.8,1.18,2.0,1.00,2.2,0.86,2.4,0.74,2.6,0.64,加工时需要,x,每改变,0.05,时,y,值,模型,将图1逆时针方向转90度，轮廓线上下对称，只需对上半部计算一个函数在插值点值。,图2 逆时针方向转90度结果,第19页,数控机床加工零件 程序,%按照表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.,-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,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),第20页,数控机床加工零件 运行结果,第21页,3.2 山区地形地貌图,已知某处山区地形选点测量坐标数据为：,x=0 0.5 1 1.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 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,第22页,山区地形地貌图 程序,原始地貌图程序：,x=0:.5:5;,y=0:.5:6;,xx,yy=meshgrid(x,y);,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 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;,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 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,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)%加密后地貌图,第23页,山区地形地貌图 结果,第24页,3.3 海底曲面图,例：在某海域测得一些点（,x,y,)处水深z由下表给出，在矩形区域(75,200),(-50,150)内画出海底曲面图形.,X,129,140,103.5,88,185.5,195,105,Y,7.5,141.5,23,147,22.5,137.5,85.5,Z,4,8,6,8,6,8,8,X,157.5,107.5,77,81,162,162,117.5,Y,-6.5,-81,3,56.5,-66.5,84,-33.5,z,9,9,8,8,9,4,9,第25页,海底曲面图 程序,clc;clf;clear 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),第26页,海底曲面图 结果,第27页,曲线拟合问题,第28页,各种拟合方法,1.线性拟合函数,2.多项式曲线拟合函数,3.稳健回归函数,4.向自定义函数拟合,5.非线性曲线拟合,第29页,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+,是p,1参数向量；,是服从标准正态分布随机干扰n1向量；,y为n1向量；,X为np矩阵。,bint返回95%置信区间。,r中为形状残差，,rint中返回每一个残差95%置信区间。,Stats向量包含R2统计量、回归F值和p值。,第30页,线性拟合,例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,第31页,2.多项式曲线拟合函数,调用格式：,p=polyfit(x,y,n),p,s=polyfit(x,y,n),说明：x,y为数据点，n为多项式阶数，返回p为幂次从高到低多项式系数向量p。矩阵s用于生成预测值误差预计。,第32页,多项式曲线拟合函数,例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),legend(原始数据,3阶曲线),第33页,多项式曲线拟合函数,也可由函数给出数据。,例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阶曲线),第34页,多项式曲线拟合函数,再用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,z,k:,x,y,b),legend(原始数据,10阶多项式),第35页,多项式曲线求值函数,：,调用格式：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,DELTA=polyconf(p,x,s,alpha),说明：Y,DELTA=polyconf(p,x,s)使用polyfit函数选项输出s给出Y95%置信区间Y DELTA。它假设polyfit函数数据输入误差是独立正态，而且方差为常数。1-alpha为置信度。,第36页,3.稳健回归函数,稳健回归是指此回归方法相对于其它回归方法而言，受异常值影响较小。,调用格式：,b=robustfit(x,y),b,stats=robustfit(x,y),说明：b返回系数预计向量；stats返回各种参数预计。,第37页,稳健回归函数,例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)+brob(2)*x,r),第38页,稳健回归函数,分析：稳健拟合(实线)对数据拟合程度好些，忽略了异常值。最小二乘拟合(点线)则受到异常值影响，向异常值偏移。,第39页,4.自定义函数拟合,对于给定数据，依据经验拟合为带有待定常数自定义函数。,所用函数：nlinfit(),调用格式：beta,r,J=nlinfit(X,y,fun,betao),说明：beta返回函数fun中待定常数；r表示残差；J表示雅可比矩阵。X,y为数据；fun自定义函数；beta0待定常数初值。,第40页,向自定义函数拟合,在化工生产中取得氯气级分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,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,第41页,自定义函数拟合,首先定义非线性函数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 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 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.3896,0.1011,方程为：,yy=0.3896+(0.49-0.3896)*exp(-0.1011*(x-8);,第42页,5.非线性曲线拟合,LSQCURVEFIT,Solves non-linear least squares problems.,功效：依据输入数据xdata和得到输出数据ydata，找到与方程F(x,xdata)最正确拟合系数。,数学模型：,调用命令：,x,resnorm,residual,exitflag,output,lambda,jacobian=,lsqcurvefit(fun,x0,xdata,ydata,lb,ub,options,p1,p2,.),第43页,举例,例、拟合函数：,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),第44页,