Matlab 8数值积分与数值微分.doc

Matlab 数值积分与数值微分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.010-6，trace控制是否展开积分过程，若为0则不展开，非0则展开，默认不展开.返回值I为积分数值；n为调用函数的次数.-例如：求03e0.5xsinx+6dx的值.先建立函数文件fesin.mfunction f=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函数求积分03e0.5xsinx+6dx的近似值，比较函数调用的次数.先建立函数文件fesin.mfunction f=fesin(x)f=exp(-0.5*x).*sin(x+(pi/6);format longI,n=quadl(fesin,0,3*pi,1e-10)I = 0.900840787818886n = 198I,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，也可以是复数，若为复数则在复平面上求积分.-例如：求积分0xsinx1+cos2xdx的数值.先编写被积函数的m文件fsx.mfunction f=fsx(x)f=x.*sin(x)./(1+cos(x).2);再调用quadgk函数I=quadgk(fsx,0,pi)I =2.4674- 例如：求积分-+xsinx1+cos2xdx的值.先编写被积函数的m文件fsx.mfunction f=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)(x1x20，引进记号：fx=fx+h-fxfx= fx-fx-hfx= fx+h2-fx-h2称上述fx、fx、fx为函数在x点处以h(h0)为步长的向前差分、向后差分、中心差分，当步长h足够小时，有：f(x)fxhfxfxhfxfxhfxh、fxh、fxh也分别被称为函数在x点处以h(h0)为步长的向前差商、向后差商、中心差商.当h足够小时，函数f(x)在x点处的导数接近于在该点的任意一种差商，微分接近于在该点的任意一种差分.2. 函数导数的求法2.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-1DX=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 = 1 1 1 1 1 1 32 16 8 4 2 1 243 81 27 9 3 1 1024 256 64 16 4 1 3125 625 125 25 5 1 7776 1296 216 36 6 1D2A1=diff(A,2,1)D2A1 = 180 50 12 2 0 0 570 110 18 2 0 0 1320 194 24 2 0 0 2550 302 30 2 0 0D2A2=diff(A,2,2)D2A2 = 0 0 0 0 8 4 2 1 108 36 12 4 576 144 36 9 2000 400 80 16 5400 900 150 25-例如：设fx=x3+2x2-x+12+6x+55+5x+2求函数f(x)的数值导数，并在同一坐标系中作出f(x)的图像.已知函数f(x)的导函数如下：fx=3x2+4x-12x3+2x2-x+12+166x+55+5编辑函数文件fun7.m和fun8.mfunction f=fun7(x)f=sqrt(x.3+2*x.2-x+12)+(x+5).(1/6)+5*x+2;function f=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,-)可以发现，最后得到的三条曲线基本重合.-练习：A. 用高斯-克朗罗德法求积分-+dx1+x2的值并讨论计算方法的精确度.（该积分值为）function f=fun9(x)f=1./(1+x.2);format longI,err=quadgk(fun9,-Inf,Inf)I = 3.141592653589793err = 5.329070518200750e-015B. 设函数fx=sinxx+cos2x 用不同的办法求该函数的数值导数，并在同一坐标系中作出f(x)的图像. 已知fx=xcosx+cosxcos2x-sinx+2sinxsin2xx+cos2x2function f=fun10(x)f=sin(x)./(x+cos(2*x);function f=fun11(x)f=(x.*cos(x)+cos(x).*cos(2*x)-sin(x)-2*sin(x).*sin(2*