MATLAB数值积分与拟合

单击此处编辑母版标题样式,单击此处编辑母版文本样式,第二级,第三级,第四级,第五级,*,数值积分,数值积分原理：,原函数不存在，采用数值积分：,1,函数,:,quad,功能,:,数值定积分，,Quad :,自适应,Simpleson,积分法。,格式,q = quad(fun,a,b) %,近似地从,a,到,b,计算函数,fun,的数值积分，误差为,10-6,。若给,fun,输入向量,x,，应返回向量,y,，即,fun,是一单值函数。,q = quad(fun,a,b,tol) %,用指定的绝对误差,tol,代替缺省误差。,tol,越大，函数计算的次数越少，速度越快，但结果精度变小。,2,程序：,fun = inline(3*x.2./(x.3-2*x.2+3);,Q1 = quad(fun,0,2),计算结果为：,Q1 =,3.7224,例,1.,计算,0,2,上如下函数的积分,3,梯形法数值积分,T = trapz(Y) %,用等距梯形法近似计算,Y,的积分。若,Y,是一向量，则,trapz(Y),为,Y,的积分；若,Y,是一矩阵，则,trapz(Y),为,Y,的每一列的积分；若,Y,是一多维阵列，则,trapz(Y),沿着,Y,的第一个非单元集的方向进行计算。,T = trapz(X,Y) %,用梯形法计算,Y,在,X,点上的积分。若,X,为一列向量，,Y,为矩阵，且,size(Y,1) = length(X),，则,trapz(X,Y),通过,Y,的第一个非单元集方向进行计算。,T = trapz(,dim) %,沿着,dim,指定的方向对,Y,进行积分。若参量中包含,X,，则应有,length(X)=size(Y,dim),。,4,例,2.,用梯形公式计算,-1,1,上,1/(1+25*X2),的积分。,X = -1:.1:1;,Y = 1./(1+25*X.2);,T = trapz(X,Y),计算结果为：,T =,0.5492,5,二元函数重积分的数值计算,1.,矩形区域上的二重积分的数值计算,格式,q = dblquad(fun,xmin,xmax,ymin,ymax),调用函数,quad,在区域,xmin,xmax, ymin,ymax,上计算二元函数,z=f(x,y),的二重积分。输入向量,x,，标量,y,，则,f(x,y),必须返回一用于积分的向量。,q = dblquad(fun,xmin,xmax,ymin,ymax,tol) %,用指定的精度,tol,代替缺省精度,10-6,，再进行计算。,q = dblquad(fun,xmin,xmax,ymin, ymax, tol,method) %,用指定的算法,method,代替缺省算法,quad,。,method,的取值有,quadl,或用户指定的、与命令,quad,与,quadl,有相同调用次序的函数句柄。,6,例如：,fun = inline(y./sin(x)+x.*exp(y);,Q = dblquad(fun,1,3,5,7),输出结果：,Q =,3.8319e+003,7,插值与拟合,Matlab,中,Lagrange,插值函数命令：,Lagrange(x,y,x0),8,例,1,：已知数据如下表，试用,Lagrange,插值多项式求,x,分别为,0.5626,0.5635,0.5645,函数的近似值。,xi 0.56160 0.56280 0.56401 0.56521,Yi,0.82741 0.82659 0.82577 0.81495,在命令窗口输入如下执行命令：,x=,0.56160;0.56280;0.56401;0.56521,;,y=0.82741;0.82659;0.82577;0.81495;,x0=,0.5626;0.5635;0.5645,;,y0=lagrang(x,y,x0),Y0=0.8265 0.8268 0.8231,plot(x,y,o,x0,y0,k*),9,2.,分段三次埃尔米特插值,Matlab,中所用命令,pchip(x,y,x0),x=,0.56160;0.56280;0.56401;0.56521,;,y=0.82741;0.82659;0.82577;0.81495;,x0=,0.5626;0.5635;0.5645,;,y0=pchip(x,y,x0),Y0=0.8267 0.8262 0.8232,10,2.,分段三次埃尔米特插值,Matlab,中所用命令,pchip(x,y,x0),x=,0.56160;0.56280;0.56401;0.56521,;,y=0.82741;0.82659;0.82577;0.81495;,x0=,0.5626;0.5635;0.5645,;,y0=spline(x,y,x0),Y0=0.8265 0.8268 0.8231,11,拟 合 问 题,已知热敏电阻数据：,温度,t(,0,C) 20.5 32.7 51.0 73.0 95.7,电阻,R(,) 765 826 873 942 1032,求,60,0,C,时的电阻,R,。,设,R=at+b,a,b,为待定系数,12,曲 线 拟 合 问 题 的 提 法,已知一组（二维）数据，即平面上,n,个点,（,x,i,y,i,) i=1,n,寻求一个函数（曲线）,y=f(x),使,f(x),在某种准则下与所有数据点最为接近，即曲线拟合得最好。,+,+,+,+,+,+,+,+,+,x,y,y=f(x),(x,i,y,i,),i,i,为点,（,x,i,y,i,),与,曲线,y=f(x),的距离,13,拟合与插值的关系,说明：,函数插值与曲线拟合都是要根据一组数据构造一个函数作为近似，由于近似的要求不同，二者的数学方法上完全不同。,实例：,下面数据是某次实验所得，希望得到,x,和,f,之间的关系？,问题：,给定一批数据点，需确定满足特定要求的曲线或曲面,解决方案：,若不要求曲线（面）通过所有数据点，而是要求它反映对象整体的变化趋势，就是,数据拟合,，又称曲线拟合或曲面拟合。,若要求所求曲线（面）通过所给所有数据点，就是,插值问题,；,14,曲线拟合问题最常用的解法,线性最小二乘法的基本思路,第一步,:,先选定一组函数,r,1,(x), r,2,(x), r,m,(x), mn,令,f(x)=a,1,r,1,(x)+a,2,r,2,(x)+ +a,m,r,m,(x),（,1,）,其中,a,1,a,2, a,m,为待定系数。,第二步,:,确定,a,1,a,2, a,m,的准则（最小二乘准则）：,使,n,个点,（,x,i,y,i,),与,曲线,y=f(x),的距离,i,的平方和最小,。,记,问题归结为，求,a,1,a,2, a,m,使,J(a,1,a,2, a,m,),最小。,15,线性最小二乘法的求解：预备知识,超定方程组,：方程个数大于未知量个数的方程组,即,Ra=y,其中,超定方程一般是不存在解的矛盾方程组。,如果有向量,a,使得 达到最小，,则称,a,为上述,超定方程的最小二乘解,。,16,线性最小二乘法的求解,定理：,当,R,T,R,可逆时，超定方程组（,3,）存在最小二乘解， 且即为方程组,R,T,Ra=R,T,y -,正则（正规）方程组,的解：,a=(R,T,R),-1,R,T,y,所以，曲线拟合的最小二乘法要解决的问题，实际上就是求以下超定方程组的最小二乘解的问题。,其中,Ra=y,（,3,）,17,线性最小二乘拟合,f(x)=a,1,r,1,(x)+ +a,m,r,m,(x),中函数,r,1,(x),r,m,(x),的选取,1.,通过机理分析建立数学模型来确定,f(x),；,+,+,+,+,+,+,+,+,+,+,+,+,+,+,+,+,+,+,+,+,+,+,+,+,+,+,+,+,+,+,f=a,1,+a,2,x,f=a,1,+a,2,x+a,3,x,2,f=a,1,+a,2,x+a,3,x,2,f=a,1,+a,2,/x,f=ae,bx,f=ae,-bx,2.,将数据,(x,i,y,i,) i=1, n,作图，通过直观判断确定,f(x),：,18,线性最小二乘拟合,1.,作多项式,f(x)=a,1,x,m,+ +a,m,x+a,m+1,拟合,可利用已有程序,:,a=polyfit(x,y,m),2.,对超定方程组,可得最小二乘意义下的解。,，用,3.,多项式在,x,处的值,y,的计算命令：,y=polyval,（,a,，,x,）,输出拟合多项式系数,a=,a,1,a,m,a,m+1,（,数组）,输入同长度,数组,X,，,Y,拟合多项式次数,19,即要求 出二次多项式,:,中 的,使得,:,例 对下面一组数据作二次多项式拟合,20,1,）输入命令,：,x=0:0.1:1;,y=-0.447,1.978,3.28,6.16,7.08,7.34,7.66,9.56,9.48,9.30,11.2;,R=(x.2), x, ones(11,1),；,A=Ry,解法,1,解超定方程的方法,2,）计算结果,： ,= -9.8108, 20.1293, -0.0317,21,2,）计算结果：,= -9.8108, 20.1293, -0.0317,解法,2,用多项式拟合的命令,MATLAB(zxec2),1,）输入命令：,x=0:0.1:1;,y=-0.447,1.978,3.28,6.16,7.08,7.34,7.66,9.56,9.48,9.30,11.2;,A=polyfit(x,y,2),z=polyval(A,x);,plot(x,y,k+,x,z,r) %,作出数据点和拟合曲线的图形,22,1. lsqcurvefit,已知,数据点,：,xdata=,（,xdata,1,，,xdata,2,，,，,xdata,n,）,ydata=,（,ydata,1,，,ydata,2,，,，,ydata,n,）,用,MATLAB,作非线性最小二乘拟合,两个求非线性最小二乘拟合的函数：,lsqcurvefit,、,lsqnonlin,。,相同点和不同点：两,个命令都要先建立,M-,文件,fun.m,，定义函数,f(x),。,lsqcurvefit,用以求含参量,x,（向量）的向量值函数,F(x,xdata)=,（,F,（,x,，,xdata,1,），,，,F,（,x,，,xdata,n,）,T,中的参变量,x(,向量,),使得,23,fun,是一个事先建立的定义函数,F(x,xdata),的,M-,文件,自变量为,x,和,xdata,x = lsqcurvefit (fun,x0,xdata,ydata);,迭代初值,已知数据点,24,lsqnonlin,用以求含参量,x,（向量）的向量值函数,f(x),=(f,1,(x),f,2,(x),f,n,(x),T,中的参量,x,，使得,最小。,其中,f,i,（,x,）,=,f,（,x,，,xdata,i,，,ydata,i,）,=F(x,xdata,i,)-ydata,i,2. lsqnonlin,已知数据点：,xdata=,（,xdata,1,，,xdata,2,，,，,xdata,n,）,ydata=,（,ydata,1,，,ydata,2,，,，,ydata,n,）,25,x=,lsqnonlin,（,fun,，,x0,）；,fun,是一个事先建立的定义函数,f(x),的,M-,文件，,自变量为,x,迭代初值,26,例,2,用下面一组数据拟合,中的参数,a,，,b,，,k,该问题即解最优化问题：,27,F(x,，,tdata)=,，,x=(a,，,b,，,k),解法,1,.,用命令,lsqcurvefit,2,）输入命令,tdata=100:100:1000,cdata=,1e-3,*4.54,4.99,5.35,5.65,5.90,6.10,6.26,6.39,6.50,6.59;,x0=0.2,0.05,0.05;,x=lsqcurvefit (curvefun1,x0,tdata,cdata),f=,curvefun1(x,tdata),1,）编写,M-,文件,curvefun1.m,function,f=curvefun1(,x,tdata,),f=x(1)+x(2)*exp(-0.02*x(3)*tdata),%,其中,x(1)=a; x(2)=b,；,x(3)=k;,28,3,）运算结果,：,f =0.0043 0.0051 0.0056 0.0059 0.0061,0.0062 0.0062 0.0063 0.0063 0.0063,x =0.0063 -0.0034 0.2542,4,）结论,：,a=0.0063, b=-0.0034, k=0.2542,29,1,）,编写,M-,文件,curvefun2.m,function,f=curvefun2(,x,),tdata=100:100:1000;,cdata=1e-03*4.54,4.99,5.35,5.65,5.90,6.10,6.26,6.39,6.50,6.59;,f=x(1)+x(2)*exp(-0.02*x(3)*tdata)- cdata,2,）输入命令,:,x0=0.2,0.05,0.05;,x=lsqnonlin(,curvefun2,x0),f= curvefun2(x),函数,curvefun2,的自变量是,x,，,cdata,和,tdata,是已知参数，故应将,cdata tdata,的值写在,curvefun2.m,中,解法,2,用命令,lsqnonlin,x=,（,a,，,b,，,k,）,30,3,）运算结果为,f =1.0e-003 *,（,0.2322 -0.1243 -0.2495 -0.2413 -0.1668 -0.0724 0.0241 0.1159 0.2030 0.2792,x =0.0063 -0.0034 0.2542,可以看出，两个命令的计算结果是相同的,。,4,）结论,：,即拟合得,a=,0.0063 b=-0.0034 k=0.2542,31,谢谢,32,