matlab diff函数用法_MATLAB数值求导与积分

4.4 数值求导与积分

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

4.4.1  导数与梯度

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

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)              %  生成随机数列

A =

     8     2    7     4     3    6     9     5    1

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

B =

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

>> C = pascal(6)

C =

     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     5   15    35    70  126

     1     6   21    56   126  252

>> D = diff(C)              % 对矩阵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

2．梯度

梯度也经常在数值计算中使用。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)

Q =

   -0.4605

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)

Q =

   -0.4605

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 =

    1.9998

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

Z =

    1.9998

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)

ans =

        0         0         0

   1.5000    2.5000    3.5000

>> cumtrapz(Y,2)

ans =

        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)

Q =

   -9.8696

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);

然后在命令窗口输入：

>> q=triplequad(@integrnd3,0,pi,0,1,-1,1)

q =

    2.0000