Matlab8数值积分与数值微分

1Matlab数值积分与数值微分Matlab数值积分1.一重数值积分的实现方法变步长辛普森法、高斯-克朗罗德法、梯形积分法1.1变步长辛普森法Matlab提供了quad函数和quadl函数用于实现变步长辛普森法求数值积分.调用格式为:[I,n]=Quad(@fname,a,b,tol,trace)[I,n]=Quadl(@fname,a,b,tol,trace)Fname是函数文件名，a,b分别为积分下限、积分上限；tol为精度控制，默认为1.0×10-6，trace控制是否展开积分过程，若为0则不展开，非0则展开，默认不展开.返回值I为积分数值；n为调用函数的次数.---------------------------------------------------------------------例如：求∫𝐞𝟎.𝟓𝐱𝐬𝐢𝐧(𝒙+𝝅𝟔)𝑑𝒙3π0的值.先建立函数文件fesin.mfunctionf=fesin(x)f=exp(-0.5*x).*sin(x+(pi/6));再调用quad函数[I,n]=quad(@fesin,0,3*pi,1e-10)I=0.9008n=365---------------------------------------------------------------------例如：分别用quad函数和quadl函数求积分∫𝐞𝟎.𝟓𝐱𝐬𝐢𝐧(𝒙+𝝅𝟔)𝑑𝒙3π0的近似值，比较函数调用的次数.先建立函数文件fesin.mfunctionf=fesin(x)f=exp(-0.5*x).*sin(x+(pi/6));formatlong2[I,n]=quadl(@fesin,0,3*pi,1e-10)I=0.900840787818886n=198[I,n]=quad(@fesin,0,3*pi,1e-10)I=0.900840787826926n=365---------------------------------------------------------------------可以发现quadl函数调用原函数的次数比quad少，并且比quad函数求得的数值解更精确.1.2高斯-克朗罗德法Matlab提供了自适应高斯-克朗罗德法的quadgk函数来求震荡函数的定积分，函数的调用格式为：[I,err]=quadgk(@fname,a,b)Err返回近似误差范围，其他参数的意义与quad函数相同，积分上下限可以是-Inf或Inf，也可以是复数，若为复数则在复平面上求积分.---------------------------------------------------------------------例如：求积分∫𝒙𝐬𝐢𝐧𝒙𝟏+𝐜𝐨𝐬𝟐𝒙𝒅𝒙π0的数值.先编写被积函数的m文件fsx.mfunctionf=fsx(x)f=x.*sin(x)./(1+cos(x).^2);再调用quadgk函数I=quadgk(@fsx,0,pi)I=2.4674---------------------------------------------------------------------例如：求积分∫𝒙𝐬𝐢𝐧𝒙𝟏+𝐜𝐨𝐬𝟐𝒙𝒅𝒙+∞−∞的值.先编写被积函数的m文件3fsx.mfunctionf=fsx(x)f=x.*sin(x)./(1+cos(x).^2);再调用quadgk函数I=quadgk(@fsx,-Inf,Inf)I=-9.0671e+017---------------------------------------------------------------------1.3梯形积分法对于一些不知道函数关系的函数问题，只有实验测得的一组组样本点和样本值，由表格定义的函数关系求定积分问题用梯形积分法，其函数是trapz函数，调用格式为：I=Traps(X,Y)X,Y为等长的两组向量，对应着函数关系Y=f(X)X=(x1,x2,…,xn)(x1x2…xn),Y=(y1,y2,…,yn)，积分区间是[x1,xn]---------------------------------------------------------------------例如：已知某次物理实验测得如下表所示的两组样本点.x1.381.562.213.975.517.799.1911.1213.39y3.353.965.128.9811.4617.6324.4129,8332.21现已知变量x和变量y满足一定的函数关系，但此关系未知，设y=f(x)，求积分∫𝒇(𝒙)𝒅𝒙13.391.38的数值.X=[1.38,1.56,2.21,3.97,5.51,7.79,9.19,11.12,13.39];Y=[3.35,3.96,5.12,8.98,11.46,17.63,24.41,29.83,32.21];I=trapz(X,Y)I=217.1033---------------------------------------------------------------------例如：用梯形积分法求积分：∫𝒆−𝒙𝒅𝒙2.51的数值.x=1:0.01:2.5;y=exp(-x);I=trapz(x,y)I=0.2858---------------------------------------------------------------------2.多重数值积分的实现4重积分的积分函数一般是二元函数f(x,y)或三元函数f(x,y,z)；形如：∫∫𝒇(𝒙,𝒚)𝒅𝒙𝒅𝒚𝑏𝑎𝑑𝑐∫∫∫𝒇(𝒙,𝒚,𝒛)𝒅𝒙𝒅𝒚𝒅𝒛𝑏𝑎𝑑𝑐𝑓𝑒Matlab中有dblquad函数和triplequad函数来对上述两个积分实现.调用格式为：I=dblquad(@fun,a,b,c,d,tol)I=triplequad(@fun,a,b,c,d,e,f,tol)Fun为被积函数，[a,b]为x的积分区间；[c,d]为y的积分区间；[e,f]为z的积分区间.Dblquad函数和triplequad函数不允许返回调用的次数，如果需要知道函数调用的次数，则在定义被积函数的m文件中增加一个计数变量，统计出被积函数被调用的次数.---------------------------------------------------------------------例如:计算二重积分I=∫∫√𝒙𝟐+𝒚𝟐𝒅𝒙𝒅𝒚𝝅𝟐−𝝅𝟐𝝅𝟐−𝝅𝟐的值.先编写函数文件fxy.mfunctionf=fxy(x,y)globalk;k=k+1;f=sqrt(x.^2+y.^2);再调用函数dblquadglobalk;k=0;I=dblquad(@fxy,-pi/2,pi/2,-pi/2,pi/2,1.0e-10)I=11.8629kk=37656---------------------------------------------------------------------例如：求三重积分∫∫∫𝒙𝒛𝒆−𝒛𝟐𝒚−𝒙𝟐𝒅𝒙𝒅𝒚𝒅𝒛𝝅𝟎𝝅𝟎𝟏𝟎的值.编写函数文件fxyz1.mfunctionf=fxyz1(x,y,z)globalj;j=j+1;5f=4*x.*z.*exp(-z.*z.*y-x.*x);调用triplequad函数editglobalj;j=0;I=triplequad(@fxyz1,0,pi,0,pi,0,1,1.0e-10)I=1.7328jj=1340978---------------------------------------------------------------------Matlab数值微分1.数值微分与差商导数的三种极限定义𝒇(𝒙)=𝐢𝟎𝒇(𝒙+)𝒇(𝒙)𝒇(𝒙)=𝐢𝟎𝒇(𝒙)𝒇(𝒙)𝒇(𝒙)=𝐢𝟎𝒇(𝒙+𝟐)𝒇(𝒙𝟐)上述公式中假设h0，引进记号：𝒇(𝒙)=𝒇(𝒙+)𝒇(𝒙)𝒇(𝒙)=𝒇(𝒙)𝒇(𝒙)𝒇(𝒙)=𝒇(𝒙+𝟐)𝒇(𝒙𝟐)称上述𝒇(𝒙)、𝒇(𝒙)、𝒇(𝒙)为函数在x点处以h(h0)为步长的向前差分、向后差分、中心差分，当步长h足够小时，有：𝒇(𝒙)𝒇(𝒙)𝒇(𝒙)𝒇(𝒙)𝒇(𝒙)𝒇(𝒙)𝒇(𝒙)、𝒇(𝒙)、𝒇(𝒙)也分别被称为函数在x点处以h(h0)为步长的向前差商、向后差商、中心差商.当h足够小时，函数f(x)在x点处的导数接近于在该点的任意一种差商，微分接近于在该点的任意一种差分.2.函数导数的求法62.1用多项式或样条函数g(x)对函数f(x)进行逼近(插值或拟合)，然后用逼近函数g(x)在点x处的导数作为f(x)在该点处的导数.2.2用f(x)在点x处的差商作为其导数.3.数值微分的实现方法Matlab中，只有计算向前差分的函数diff，其调用格式为：·DX=diff(X):计算向量X的向前差分，DX(i)=X(i+1)-X(i),i=1,2,…,n-1·DX=diff(X,n):计算向量X的n阶向前差分，例如diff(X,2)=diff(diff(X))·DX=diff(A,n,dim):计算矩阵A的n阶向前差分，dim=1(默认值)按列计算差分，dim=2按行计算差分.---------------------------------------------------------------------例如：生成6阶范德蒙德矩阵，然后分别按行、按列计算二阶向前差分A=vander(1:6)A=111111321684212438127931102425664164131256251252551777612962163661D2A1=diff(A,2,1)D2A1=180501220057011018200132019424200255030230200D2A2=diff(A,2,2)D2A2=000084211083612457614436920004008016540090015025---------------------------------------------------------------------7例如：设𝒇(𝒙)=√𝒙+𝟐𝒙𝟐𝒙+𝟏𝟐+√(𝒙+𝟓)𝟓𝟔+𝟓𝒙+𝟐求函数f(x)的数值导数，并在同一坐标系中作出f’(x)的图像.已知函数f(x)的导函数如下：𝒇(𝒙)=𝒙𝟐+𝒙𝟏𝟐√𝒙+𝟐𝒙𝟐𝒙+𝟏𝟐+𝟏𝟔√(𝒙+𝟓)𝟓𝟔+𝟓编辑函数文件fun7.m和fun8.mfunctionf=fun7(x)f=sqrt(x.^3+2*x.^2-x+12)+(x+5).^(1/6)+5*x+2;functionf=fun8(x)f=(3*x.^2+4*x-1)/2./sqrt(x.^3+2*x.^2-x+12)+1/6./(x+5).^(5/6)+5;x=-3:0.01:3;p=polyfit(x,fun7(x),5);用5次多项式拟合曲线dp=polyder(p);对拟合多项式进行求导dpx=polyval(dp,x);对dp在假设点的求函数值dx=diff(fun7([x,3.01]))/0.01;直接对dx求数值导数gx=fun8(x);求函数f的函数在假设点的导数plot(x,dpx,x,dx,'.',x,gx,'-')可以发现，最后得到的三条曲线基本重合.------------------------------------------------------------------