第七章 数值积分、微分

1、Matlaband Engineering Calculation第7章MATLAB数值积分与微分7.1 数值积分7.2 数值微分XiamenUniversityMatlaband Engineering Calculation7.1 数值积分7.1.1 数值积分基本原理求解定积分的数值方法多种多样，如简单的梯形法、辛普生(Simpson) 法、牛顿柯特斯(Newton-Cotes)法等都是经常采用的方法。它们的基本思想都是将整个积分区间a,b分成n个子区间xi,xi+1，i=1,2,n，其中x1=a，xn+1=b。这样求定积分问题就分解为求和问题。XiamenUniversityMatlab

2、and Engineering Calculation7.1.2 数值积分的实现方法(Quadrature)1自适应辛普生法基于变步长辛普生法，MATLAB给出了quad函数来求定积分。该函数的调用格式为：I,n=quad(fun,a,b,tol,trace)fun是被积函数名、句柄或内联函数对象;a和b分别是定积分的下限和上限;tol用来控制积分精度，缺省时取tol=10-6;trace控制是否展现积分过程，若取非0则展现积分过程，取0则不展现，缺省时取trace=0I即定积分值，n为被积函数的调用次数。XiamenUniversityMatlaband Engineering Calcul

3、ation例7-1 求定积分。(1) 建立被积函数文件fesin.m。function f=fesin(x)f=exp(-0.5*x).*sin(x+pi/6);(2) 调用数值积分函数quad求定积分。S,n=quad('fesin',0,3*pi)S =0.9008n =77XiamenUniversityMatlaband Engineering Calculation2牛顿柯特斯法基于牛顿柯特斯法，MATLAB给出了quad8函数来求定积分。该函数的调用格式为：I,n=quad8(fun,a,b,tol,trace)其中参数的含义和quad函数相似。 该函数可以更精确地

4、求出定积分的值，且一般情况下函数调用的步数明显小于quad函数，从而保证能以更高的效率求出所需的定积分值。XiamenUniversityMatlaband Engineering Calculation例7-2 求定积分。(1) 被积函数文件fx.m。f=inline('x.*sin(x)./(1+cos(x).*cos(x) ');I=quad8(f,0,pi)I =2.4674XiamenUniversityMatlaband Engineering Calculation2高阶方法MATLAB给出了quadl函数来替代quad8求定积分I,n=quadl(fun,a,b

5、,tol,trace)其中参数的含义和quad函数相似。 该函数可以更精确地求出定积分的值，且一般情况下函数调用的步数明显小于quad函数，从而保证能以更高的效率求出所需的定积分值。XiamenUniversityMatlaband Engineering Calculation例7-3 分别用quad函数、quad8和quadl函数求定积分的近似值，并在相同的积分精度下，比较函数的调用次数。调用函数quad求定积分：format long;fx=inline('exp(-x)');I,n=quad(fx,1,2.5,1e-10)I =0.28579444254766n =65

6、调用函数quad8求定积分：I,n=quad8(fx,1,2.5,1e-10)I =0.28579444254754n =33调用函数quadl求定积分：I,n=quadl(fx,1,2.5,1e-10)I =0.28579444254754n =33XiamenUniversityMatlaband Engineering Calculation7.1.3 二重定积分的数值求解二重积分可由单重积分函数来构建。1、自己利用单重积分函数编程进行二重积分计算；2、使用MATLAB提供的dblquad函数求二重定积分的数值解。I=dblquad(f,a,b,c,d,tol,method)该函数求f(

7、x,y)在a,b×c,d区域上的二重定积分。参数tol的用法与函数quad完全相同。method 用于指定用哪个数值积分函数，如quad8,quadl或用户自己定义的积分函数XiamenUniversityMatlaband Engineering Calculation例7-5 计算二重定积分(1) 建立一个函数文件fxy.m：function f=fxy(x,y)global ki;ki=ki+1; %ki用于统计被积函数的调用次数f=exp(-x.2/2).*sin(x.2+y);(2) 调用dblquad函数求解。global ki;ki=0;I=dblquad('f

8、xy',-2,2,-1,1)kiI =1.57449318974494ki=1038XiamenUniversityMatlaband Engineering Calculation？已知二自变量函数在矩形网格下采样点的值，如何计算其积分？x=-1:0.1:1;y=-2:0.1:2;X,Y=meshgrid(x,y);Z=sqrt(X.2+Y.2);mesh(X,Y,Z)I1=trapz(y,trapz(x,Z,2)I2=trapz(x,trapz(y,Z),2)XiamenUniversityMatlaband Engineering Calculation7.1.4 三重定积分的数

9、值求解三重积分可由单重积分函数来构建。1、自己利用单重和二重积分函数编程进行三重积分计算；2、使用MATLAB提供的triplequad函数求三重定积分的数值解。triplequad(fun,xmin,xmax,ymin,ymax,zmin,zmax,tol,method)e.g.Q = triplequad(inline('y*sin(x)+z*cos(x)'),0,pi,0,1,-1,1)或Q = triplequad(integrnd,0,pi,0,1,-1,1)where integrnd.mis the M-filefunction f = integrnd(x,y,

10、z)f = y*sin(x)+z*cos(x);一重积分计算平面图形面积，二重积分计算体积，三重积分的物理意义是什么？如何使用triplequad计算一个非立方体的质量?XiamenUniversityMatlaband Engineering Calculation？已知三自变量函数在空间立方网格下采样点的值，如何计算其积分？x=-1:0.1:1;y=-2:0.1:2;z=1:0.1:2;X,Y,Z=meshgrid(x,y,z);F=sqrt(X.2+Y.2+Z.2);slice(X,Y,Z,F, -0.5 0 0.5,-1 0 1, 1.5)shading interpdaspect(1

11、 1 1)camlight;I1=trapz(z,trapz(y,trapz(x,Z,2),3)I2=trapz(z,trapz(x,trapz(y,Z),2),3)I3=trapz(x,trapz(y,trapz(z,Z,3),2)XiamenUniversityMatlaband Engineering Calculation例6-6 生成以向量V=1,2,3,4,5,6为基础的范得蒙矩阵，按列进行差分运算。命令如下：V=vander(1:6)DV=diff(V) %计算V的一阶差分XiamenUniversityMatlaband Engineering Calculation例6-7

12、用不同的方法求函数f(x)的数值导数，并在同一个坐标系中做出f'(x)的图像。程序如下：f=inline('sqrt(x.3+2*x.2-x+12)+(x+5).(1/6)+5*x+2');g=inline('(3*x.2+4*x-1)./sqrt(x.3+2*x.2-x+12)/2+1/6./(x+5).(5/6)+5');x=-3:0.01:3;p=polyfit(x,f(x),5); %用5次多项式p拟合f(x)dp=polyder(p); %对拟合多项式p求导数dpdpx=polyval(dp,x); %求dp在假设点的函数值dx=diff(f(

13、x,3.01)/0.01; %直接对f(x)求数值导数gx=g(x); %求函数f的导函数g在假设点的导数plot(x,dpx,x,dx,'.',x,gx, 'r-'); %作图XiamenUniversityMatlaband Engineering Calculation标量场梯度计算(gradient)Fx,Fy,Fz,. = gradient(F,h1,h2,.)Fx,Fy,Fz,. = gradient(F)F是一个矩阵，F的维数表示自变量的个数n，返回n个大小相同的矩阵v = -2:0.2:2;x,y = meshgrid(v);z = x .* exp(-x.2 -y.2);px,py = gradient(z,.2,.2);contour(v,v,z), hold on, quiver(v,v,px,py), hold