专题六数值微积分与方程求解6.1 数值微分与数值积分☐数值微分☐数值积分1．数值微分（1）数值差分与差商微积分中，任意函数f(x)在x 0点的导数是通过极限定义的：hx f h x f x f h )()(lim)('00-+=→hh x f x f x f h )()(lim)('0000--=→hh x f h x f x f h )2/()2/(lim)('0--+=→)()()(000x f h x f xf -+=∆)()()(0h x f x f x f --=∇)2/()2/()(0h x f h x f x f --+=δ如果去掉极限定义中h 趋向于0的极限过程，得到函数在x 0点处以h （h>0）为步长的向前差分、向后差分和中心差分公式： 向前差分： 向后差分： 中心差分：函数f(x)在点x 0的微分接近于函数在该点的差分，而f 在点x 的导数接近于函数在该点的差商。

hx f h x f xf )()(≈)('000-+hh x f x f x f )()(≈)('000--hh x f h x f x f )2/()2/(≈)('0--+向前差商： 向后差商： 中心差商：当步长h 充分小时，得到函数在x 0点处以h （h>0）为步长的向前差商、向后差商和中心差商公式：（2）数值微分的实现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，按行计算差分。

注意：diff函数计算的是向量元素间的差分，故差分向量元素的个数比原向量少了一个。

同样，对于矩阵来说，差分后的矩阵比原矩阵少了一行或一列。

另外，计算差分之后，可以用f(x)在某点处的差商作为其导数的近似值。

例1 设f(x)=sin x，在[0，2π]范围内随机采样，计算f'(x)的近似值，并与理论值f'(x)=cos x进行比较。

>> x=[0,sort(2*pi*rand(1,5000)),2*pi];>> y=sin(x);>> f1=diff(y)./diff(x);>> f2=cos(x(1:end-1));>> plot(x(1:end-1),f1,x(1:end-1),f2);>> d=norm(f1-f2)d =0.04332．数值积分（1）数值积分基本原理在有些情况下，应用牛顿—莱布尼兹公式有困难，例如，当被积函数的原函数无法用初等函数表示，或被积函数是用离散的表格形式给出的。

这时就需要用数值解法来求定积分的近似值。

在高等数学中，计算定积分依靠微积分基本定理，只要找到被积函数f(x)的原函数大F(x)，则可用牛顿—莱布尼兹（Newton-Leibniz ）公式：⎰-=ba a Fb F xx f)()(d )(求定积分的数值方法多种多样，如梯形法、辛普森（Simpson ）法、高斯求积公式等。

它们的基本思想都是将积分区间[a ，b]分成n 个子区间[x i ，x i+1]，i=1，2，…，n ，其中x 1=a ，x n+1=b ，这样求定积分问题就分解为下面的求和问题：在每一个小的子区间上定积分的值可以近似求得,从而避免了牛顿—莱布尼兹公式需要寻求原函数的困难。

⎰∑⎰=+==ba nix x i ixx f xx f S11d )(d )(（2）数值积分的实现☐基于自适应辛普森方法[I,n]=quad(filename,a,b,tol,trace)☐基于自适应Gauss-Lobatto方法[I,n]=quadl(filename,a,b,tol,trace)其中，filename是被积函数名；a和b分别是定积分的下限和上限，积分限[a，b]必须是有限的，不能为无穷大（Inf）；tol用来控制积分精度，默认时取tol=10-6；trace控制是否展现积分过程，若取非0则展现积分过程，取0则不展现，默认时取trace=0；返回参数I即定积分的值，n为被积函数的调用次数。

例2 分别用quad 函数和quadl 函数求定积分的近似值，并在相同的积分精度下，比较被积函数的调用次数。

x d x1412⎰+>> format long>> f=@(x) 4./(1+x.^2); >> [I,n]=quad(f,0,1,1e-8) I =3.141592653733437 n = 61>> [I,n]=quadl(f,0,1,1e-8) I =3.141592653589806 n = 48>> (atan(1)-atan(0))*4 ans =3.141592653589793 >> format short基于全局自适应积分方法I=integral(filename,a,b)其中，I是计算得到的积分；filename是被积函数；a和b分别是定积分的下限和上限，积分限可以为无穷大。

dxx x I e ⎰-=12ln 11例3 求定积分。

被积函数文件fe.m ：function f=fe(x)f=1./(x.*sqrt(1-log(x).^2));>> I=integral(@fe,1,exp(1))I =1.5708基于自适应高斯-克朗罗德方法[I,err]=quadgk(filename,a,b)其中，err返回近似误差范围，其他参数的含义和用法与quad函数相同。

积分上下限可以是无穷大（−Inf或Inf），也可以是复数。

如果积分上下限是复数，则quadgk函数在复平面上求积分。

例4 求定积分。

⎰∞+ 2 2d x1sin 1πx x 被积函数文件fsx.m ：function f=fsx(x)f=sin(1./x)./x.^2;>> I=quadgk(@fsx,2/pi,+Inf)I =1.0000基于梯形积分法已知（xi ，yi）（i=1，2，…，n），且a=x1<x2<…<xn=b，求近似值。

⎰=b a xfI(x)dI=trapz(x,y)其中，向量x、y定义函数关系y=f(x)。

trapz 函数采用梯形积分法则，积分的近似值为：其中， 。

sum(diff(x).*(y(1:end-1)+y(2:end))/2)I = ℎi n−1i=1y i+1+y i ℎi =x i+1−x i 可用以下语句实现：例5 设x=1:6，y=[6,8,11,7,5,2],用trapz函数计算定积分。

>> x=1:6;>> y=[6,8,11,7,5,2];>> plot(x,y,'-ko');>> grid on>> axis([1,6,0,11]);>> I1=trapz(x,y)I1 =35>> I2=sum(diff(x).*(y(1:end-1)+y(2:end))/2)I2 =35（3）多重定积分的数值求解I=integral3(filename,a,b,c,d,e,f)I=triplequad(filename,a,b,c,d,e,f,tol)(,)d d d b c af x y x y ⎰⎰(,,)d d d f d b e c a f x y z x y z⎰⎰⎰☐求二重积分的数值解： I=integral2(filename,a,b,c,d) I=quad2d(filename,a,b,c,d)I=dblquad(filename,a,b,c,d,tol) ☐求三重积分的数值解：>> f1=@(x,y) exp(-x.^2/2).*sin(x.^2+y); >> I1=quad2d(f1,-2,2,-1,1)I1 =1.574>> f2=@(x,y,z) 4*x.*z.*exp(-z.*z.*y-x.*x); >> I2=integral3(f2,0,pi,0,pi,0,1) I2 =1.7328 ⎰⎰---+11 2 2 22/d d )sin(e 2y x y x x 例6 分别求二重积分和三重积分。

2210004e d d d z y x xz x y zππ--⎰⎰⎰。