数学建模讲座插值和拟合.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.41 网格节点插值法,已知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 041 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页,