用Matlab计算数值微分

导论

        函数微分的概念这里不再赘述。阅读此文前希望读者能够知晓前向差分、后向差分与中心差分的概念。

        本文以函数f(x)为例进行说明如何利用Matlab计算数值微分

diff函数计算微分

        如果输入一个长度为n的一维向量，diff函数返回一个长度为n-1的向量，其中包含原向量相邻元素的差。

clear
clc
f = @(x) 0.2 + 25*x - 200*x.^2 + 675*x.^3 - 900*x.^4 + 400* x.^5;
x = 0:0.1:0.8;
y = f(x);
d = diff(y) ./ diff(x);
n = length(x);
xm = (x(1:n-1) + x(2:n)) / 2;

xa = 0:0.01:0.8;
ya = 25 - 400*xa + 3*675*xa.^2 - 4*900*xa.^3 + 5*400*xa.^4;
figure(1)
plot(xm, d, '*', xa, ya)

结果如下：

        请注意：diff(y)/diff(x)求得的是每个小区间中点处的导数（且为中心差分格式），在代码中画图部分使用xm作为横坐标也正是其体现。

gradient函数计算微分

        输入一个长度为n的向量，返回一个长度为n的向量，返回向量的第一个值为输入向量前两个值之差除以数据点间距，返回向量的最后一个值为输入向量的最后两个值之差除以数据点间距（这里可以认为在端点处使用向前和向后差分）；对于中间的值返回的是与其相邻的两个值的中心型差分（注意，这里已经做了差分运算，除以2*间距）

clear
clc
f = @(x) 0.2 + 25*x - 200*x.^2 + 675*x.^3 - 900*x.^4 + 400* x.^5;
x = 0:0.05:0.8;
n = length(x);
y = f(x);
dy = gradient(y, 0.05);
x = [(x(1)+x(2))/2, x(2:end-1), (x(end-1)+x(end))/2];
figure(2)
plot(x, dy, '*', xa, ya)

• fx = gradient(f, h);其中h为数据点之间的间距，若缺省则默认为1
• 返回的前后两点导数值也可以认为是第一个和最后一个区间中点处的中心差分，程序正是按照这个思路进行编写的。

向量场的可视化

上边提到的gradient函数还特别适用于计算矩阵的偏导数，其用法如下：

[fx, fy] = gradient(f, h);

fx对应于x方向(列)的差值，fy对应于y方向（行）的差值，h为数据点之间的间距。

clear
clc

f = @(x, y) y - x - 2*x.^2 - 2*x.*y - y.^2;
[x, y] = meshgrid(-2:0.25:0, 1:0.25:3);
z = f(x, y);

[fx, fy] = gradient(z, 0.25);
cs = contour(x, y, z);clabel(cs);hold on

quiver(x, y, -fx, -fy); hold off
set(gca, 'Fontname', 'Times')