4.4 数值求导与积分

在数学计算中,积分和求导是最常见的运算。

4.4.1  导数与梯度

导数的数值计算是数值计算的基本操作之一。如牛顿法求根、微分方程求解、泰勒级数展开等,都离不开导数。

在MATLAB中,diff函数被用来计算数值差分或者符号导数。本小节只介绍diff函数如何用来计算差分,符号导数的计算将在下一章介绍。

diff函数的调用语法如下。

(1)Y = diff(X):求X相邻行元素之间的一阶差分。

(2)Y = diff(X,n):求X相邻行元素之间的n阶差分。

(3)Y = diff(X,n,dim):在dim指定维上求X相邻行元素之间的n阶差分。

【例4-31】  diff函数应用示例。

>> rand('state',0)            %  设置随机种子

>> A=randperm(9)              %  生成随机数列

8     2    7     4     3    6     9     5    1

>> B = diff(A)                %  求数列的差分

-6     5   -3    -1     3    3    -4    -4

>> C = pascal(6)

1     1    1     1    1     1

1     2    3     4     5    6

1     3    6    10    15   21

1     4   10    20    35   56

1     6   21    56   126  252

>> D = diff(C)              % 对矩阵C列方向各元素进行差分计算

0    1     2     3    4     5

0     1    3     6    10   15

0     1    4    10    20   35

0     1    5    15    35   70

0     1    6    21    56  126

梯度也经常在数值计算中使用。MATLAB提供了gradient函数来进行梯度计算。gradient函数的调用语法如下。

(1)FX = gradient(F):返回F的一维数值梯度,F是一个向量。

(2)[FX,FY] = gradient(F):返回二维数值梯度的x和y部分,F是一个矩阵。

(3)[FX,FY,FZ,...] = gradient(F):求高维矩阵F的数值梯度。

(4)[...] = gradient(F,h):h是一个标量,用于指定各个方向上点之间的间距。

(5)[...] = gradient(F,h1,h2,...):指定各个方向上的间距。

【例4-32】  梯度求解示例。

>> 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 off

%  绘制等高线和梯度方向

运行以上命令,可以得到如图4-3所示的结果。

图4-3 梯度计算示例

4.4.2 一元函数的数值积分

1.quad函数

函数quad采用自适应Simpson方法计算积分,特点是精度较高,较为常用。quad函数的调用语法如下。

(1)q=quad(fun,a,b):计算函数fun在a到b区间内的数值积分,其中fun是一个函数句柄,默认误差为10-6。若给fun输入向量x,应返回y,即fun是一个单值函数。

(2)q=quad(fun,a,b,tol):用指定的绝对误差tol代替默认误差。tol越大,函数计算的次数越少,速度越快,但相应计算结果的精度会越低。

(3)[q,fcnt]=quad(fun,a,b):计算函数fun的数值积分,同时返回函数计算的次数fcnt。

【例4-33】  使用quad函数求的数值积分。

为了求函数的积分,首先要在工作目录下建立函数文件,不妨命名为myfun.m,该函数文件的内容如下:

function y = myfun(x)

y = 1./(x.^3-2*x-5);

函数文件的创建与使用,将在第6章介绍。在建立好函数文件之后,需要传递相应的函数句柄@myfun到quad函数,同时需要指定上下限。

>> Q = quad(@myfun,0,2)

2.quadl函数

函数quadl采用自适应Lobatto方法计算积分,特点是精度较高,最为常用。quadl函数的主要调用语法如下。

(1)q=quadl(fun,a,b):计算函数fun在a到b区间内的数值积分,其中fun是一个函数句柄,误差为10-6。若给fun输入向量x,应返回y,即fun是一个单值函数。

(2)q=quadl(fun,a,b,tol):用指定的绝对误差tol代替默认误差。tol越大,函数计算的次数越少,速度越快,但相应计算结果的精度会越低。

【例4-34】  使用quadl函数求的数值积分。

本例在前例建立的myfun文件基础上进行计算。

>> Q = quadl(@myfun,0,2)

3.trapz函数

trapz函数使用梯形法进行积分,特点是速度快,但精度差。trapz函数的调用语法如下。

(1)T=trapz(y):用等距梯形法近似计算y的积分。若y是一个向量,则trapz(y)为y的积分;若y是一个矩阵,则trapz(y)为y的每一列的积分。

(2)trapz(x,y):用梯形法计算y在x点的积分。

(3)T=trapz(…,dim):沿着dim指定的维对y进行积分。若参量当中包含x,则应有length(x)=size(y,dim)。

【例4-35】  trapz函数使用示例。计算。

通过数学推导可知,的精确值为2。现在我们以trapz函数进行计算以对比结果。首先需要划分梯形法使用的均匀区间:

>> X = 0:pi/100:pi;

>> Y = sin(X);        %  被积函数

>> Z = trapz(X,Y)

>> Z = pi/100*trapz(Y)

4.cumtrapz函数

cumtrapz函数用于求累积的梯形数值积分,其调用语法如下。

(1)Z=cumtrapz(Y):通过梯形积分法计算单位步长时Y的累积积分值,若步长不是一个单位量,则算出的Z值还应该乘以步长。当Y为向量时,返回该向量的累积积分;若Y为矩阵,则返回该矩阵每列元素的累积积分。

(2)Z=cumtrapz(X,Y):采用梯形积分求Y对X的积分,注意X与Y的长度必须相等。

(3)Z=cumtrapz(X,Y,dim)或Z=cumtrapz(Y,dim):对Y的dim维求积分,X的长度必须等于size(Y,dim)。

【例4-36】  cumtrapz函数应用示例。

>> Y = [0 1 2; 3 4 5];

>> cumtrapz(Y,1)

0         0         0

1.5000    2.5000    3.5000

>> cumtrapz(Y,2)

0    0.5000    2.0000

0    3.5000    8.0000

4.4.3 二重积分的数值计算

二重积分的数值计算可由函数dblquad实现。dblquad函数有以下几种语法格式。

(1)q=dblquad(fun,xmin,xmax,ymin,ymax):调用函数quad在区域xmin≤x≤xmax, ymin≤y≤ymax上计算二元函数z=fun(x,y)的二重积分。函数fun需要满足条件:输入向量x、标量y,则f(x,y)必须返回一个用于积分的向量。

(2)q=dblquad(fun,xmin,xmax,ymin,ymax,tol):用指定的精度tol代替默认精度10-6,再进行计算。

(3)q=dblquad(fun,xmin,xmax,ymin,ymax,tol,method):用指定的计算方法method代替默认算法quad。method的取值有@quadl或用户指定的与命令quad和quadl有相同调用次序的函数句柄。

【例4-37】  应用dblquad函数求重积分示例。

首先在工作目录下建立一个M文件integrnd.m,该文件内容为:

function z = integrnd(x, y)

z = y*sin(x)+x*cos(y);

然后调用dblquad函数求重积分:

>> Q = dblquad(@integrnd,pi,2*pi,0,pi)

4.4.4 三重积分的数值计算

MATLAB提供了triplequad函数来对三重积分进行数值计算。该函数的调用语法如下。

(1)q=triplequad(fun,xmin,xmax,ymin,ymax,zmin,zmax):调用函数quad在区域[xmin,xmax, ymin,ymax,zmin,zmax]上进行三重积分运算。fun参数是一个M文件函数或匿名函数的句柄。

(2)q=triplequad(fun,xmin,xmax,ymin,ymax,zmin,zmax,tol):用指定的精度tol代替默认精度10-6进行计算。

(3)q=triplequad(fun,xmin,xmax,ymin,ymax,zmin,zmax,tol,method):用指定的计算方法method代替默认算法quad。method的取值有@quadl或用户指定的与命令quad与quadl有相同调用次序的函数句柄。

【例4-38】  用triplequad函数求三重积分。

首先需要在工作目录下建立函数M文件integrnd3.m,该文件内容为:

function z=integrnd3(x,y,z)

z=y*sin(x)+z*cos(x);

然后在命令窗口输入: