跳转至

MATLAB 微分

微分是描述函数变化率的核心数学工具，在科学与工程计算中无处不在。MATLAB 作为一个强大的计算平台，提供了两大类截然不同的方法来进行微分运算：数值微分与符号微分。理解它们各自的原理、适用场景和局限性，是准确高效地解决问题的关键。

  • 数值微分 (Numerical Differentiation):

    • 适用场景: 当您只有一组离散的数据点（例如，实验测量值），而没有函数的解析表达式时。
    • 核心思想: 基于有限差分法，用离散点的差商来近似导数。
    • 特点: 结果是近似值，对数据噪声非常敏感。

  • 符号微分 (Symbolic Differentiation):

    • 适用场景: 当你知道函数的精确数学表达式时。
    • 核心思想: 应用微积分的求导法则，得到导函数的精确解析表达式。
    • 特点: 结果是精确的符号表达式，完全无误差。需要 Symbolic Math Toolbox。

数值微分 (Numerical Differentiation)

原理分析：有限差分法

数值微分的理论基础源于导数的极限定义：

\[ f'(x) = \lim_{\dt{x} \to 0} \frac{f(x+\dt{x}) - f(x)}{h} \]

有限差分法通过选取一个很小但非零的步长 \(h\,(h>0)\) 来近似这个极限。根据采样点的不同，主要分为三种形式：

  • 前向差分 (Forward Difference): 使用当前点和前方一个点来近似导数。其公式为：

    \[ f'(x) \approx \frac{f(x+h) - f(x)}{h} \]

    它是一阶精度近似，误差与步长 \(h\) 成正比，记为 \(O(h)\)

  • 后向差分 (Backward Difference): 使用当前点和后方一个点来近似导数。其公式为：

    \[ f'(x) \approx \frac{f(x) - f(x-h)}{h} \]

    它同样是一阶精度近似，\(O(h)\)

  • 中心差分 (Central Difference): 使用当前点前后各一个点来近似导数。其公式为：

    \[ f'(x) \approx \frac{f(x+h) - f(x-h)}{2h} \]

    中心差分通常比前向或后向差分更精确。它是二阶精度近似，误差与 \(h^2\) 成正比，记为 \(O(h^2)\)

函数展示

diff - 基础差分

  • 功能 计算向量或矩阵中相邻元素之间的差值。它本身是差分运算，除以步长后可近似导数。

  • 语法 

    Y = diff(X)
Y = diff(X, n)
Y = diff(X, n, dim)

  • 语法说明

    • diff(X): 计算向量 X 的一阶前向差分。如果 X 是矩阵，则对每一列进行操作。返回结果的长度会减 1。
    • diff(X, n): 计算 n 阶差分。例如，diff(X, 2) 等价于 diff(diff(X))。
    • diff(X, n, dim): 沿着指定的维度 dim 计算差分。

  • 示例 

    x = 0:0.1:pi;
y = sin(x);

% 使用前向差分近似导数
dy_dx_approx = diff(y) ./ diff(x);

% 注意：结果的长度比原始向量少1，且导数值对应于原x点的中点
x_mid = x(1:end-1) + 0.05;

figure;
plot(x, cos(x), 'r--', 'DisplayName', '解析导数 cos(x)');
hold on;
plot(x_mid, dy_dx_approx, 'bo', 'DisplayName', '数值导数 (diff)');
legend;
title('diff - 基础差分求导');

gradient - 推荐的数值梯度函数

  • 功能 计算一维、二维或多维数组的数值梯度。它对内部点使用更精确的中心差分，对端点使用单边差分，并返回与输入数组同等大小的结果。

  • 语法 

    FX = gradient(F)
FX = gradient(F, h)
[FX, FY] = gradient(F, hx, hy)
[FX, FY, FZ, ...] = gradient(F, hx, hy, hz, ...)

  • 语法说明

    • gradient(F): 假设采样点间距为 1，计算 F 的梯度。
    • gradient(F, h): 为一维向量 F 指定均匀的采样间距 h。
    • [FX, FY] = gradient(F, hx, hy): 为二维矩阵 F 指定 x 和 y 方向的均匀采样间距 hx 和 hy。FX 是 F 对 x 的偏导数 (\(\partial F / \partial x\))，FY 是 F 对 y 的偏导数 (\(\partial F / \partial y\))。

  • 示例 (二维) 

    [X, Y] = meshgrid(-2:0.2:2);
Z = X .* exp(-X.^2 - Y.^2);

% 计算 Z 的梯度，[PX, PY] 分别是 dZ/dX 和 dZ/dY
% 步长 hx 和 hy 均为 0.2
[PX, PY] = gradient(Z, 0.2, 0.2);

figure;
contour(X, Y, Z); % 绘制原始函数等高线
hold on;
quiver(X, Y, PX, PY); % 叠加梯度向量场
hold off;
title('gradient - 二维求导（梯度场）');
axis equal;

del2 - 离散拉普拉斯算子

  • 功能 计算数组的离散拉普拉斯算子。对于二维函数 \(u(x, y)\)，其定义为：

    \[ \nabla^2 u = \frac{\partial^2 u}{\partial x^2} + \frac{\partial^2 u}{\partial y^2} \]

    常用于图像的边缘检测和物理学中的偏微分方程。

  • 语法 

    L = del2(U)
L = del2(U, h)
L = del2(U, hx, hy, ...)

  • 语法说明

    • del2(U): 假设所有方向的步长均为 1。
    • del2(U, h): 指定一个均匀的步长 h。
    • del2(U, hx, hy, ...): 为每个维度指定不同的步长。

  • 示例 

    [X, Y] = meshgrid(-2:0.1:2);
Z = peaks(X, Y);

% 计算 Z 的拉普拉斯，步长为 0.1
L = del2(Z, 0.1);

figure;
surf(X, Y, L);
title('del2 - peaks 函数的离散拉普拉斯');

符号微分 (Symbolic Differentiation)

原理分析：应用求导法则

符号微分不依赖于离散的数据点，而是直接对数学表达式应用微积分的求导法则（如链式法则、乘积法则等）。这需要 Symbolic Math Toolbox。其核心优势在于精确性，结果是导函数的精确数学表达式，而非数值近似。

函数展示

diff (符号)

  • 功能 对符号表达式求导。可以指定求导变量和阶数。

  • 语法 

    df = diff(f)
df = diff(f, var)
df = diff(f, n)
df = diff(f, var, n)

  • 语法说明

    • diff(f): 对默认符号变量（通常是 x）求一阶导数。
    • diff(f, var): 对指定的符号变量 var 求一阶导数。
    • diff(f, n): 对默认符号变量求 n 阶导数。
    • diff(f, var, n): 对指定的符号变量 var 求 n 阶导数。

  • 示例 (常微分与偏微分) 

    syms x y
f = sin(x^2) * exp(-y);

% 对 x 求一阶偏导
df_dx = diff(f, x);

% 对 y 求二阶偏导
d2f_dy2 = diff(f, y, 2);

% 先对 x 再对 y 求混合偏导
d2f_dxdy = diff(diff(f, x), y);

disp('∂f/∂x:'); disp(df_dx);
disp('∂²f/∂y²:'); disp(d2f_dy2);
disp('∂²f/∂y∂x:'); disp(d2f_dxdy);

gradient (符号)

  • 功能 计算一个标量符号函数相对于一个向量变量的梯度。结果是一个符号向量。

    \[ \nabla f = \left[ \frac{\partial f}{\partial x_1}, \frac{\partial f}{\partial x_2}, \dots, \frac{\partial f}{\partial x_n} \right] \]

  • 语法 

    g = gradient(f)
g = gradient(f, vars)

  • 语法说明

    • gradient(f): 计算 f 相对于其所有符号变量的梯度。
    • gradient(f, vars): 计算 f 相对于 vars 中指定的符号变量向量的梯度。

  • 示例 

    syms x y z
f = x^2 * sin(y) + z^3;

% 计算梯度 ∇f
grad_f = gradient(f, [x, y, z]);

disp('函数 f 的梯度:');
disp(grad_f); % 结果: [ 2*x*sin(y), x^2*cos(y), 3*z^2 ]

jacobian (符号)

  • 功能 计算一个向量函数的雅可比矩阵，即所有一阶偏导数构成的矩阵。

  • 语法 

    J = jacobian(f, v)

  • 语法说明

    • jacobian(f, v): 计算向量函数 f 相对于向量 v 的雅可比矩阵。

  • 示例 

    syms r theta

% 极坐标到笛卡尔坐标的变换
x = r * cos(theta);
y = r * sin(theta);

F = [x; y];
vars = [r; theta];

% 计算雅可比矩阵
J = jacobian(F, vars);

disp('极坐标变换的雅可比矩阵:');
disp(J);

hessian (符号)

  • 功能 计算一个标量函数的二阶偏导数（海森）矩阵。它在多元函数的极值判断和优化问题中至关重要。

  • 语法 

    H = hessian(f)
H = hessian(f, vars)

  • 语法说明

    • hessian(f, vars): 计算函数 f 相对于 vars 中指定的变量向量的海森矩阵。

  • 示例 

    syms x y
f = x^3 + y^3 + 3*x*y;

% 计算海森矩阵
H = hessian(f, [x, y]);

disp('函数 f 的海森矩阵:');
disp(H);

高级与工具箱特定函数

differentiate (Curve Fitting Toolbox)

  • 功能 对拟合曲线 (cfit 对象) 或原始数据 (x, y) 求导。它能有效地处理含噪数据，因为它是在平滑的拟合曲线上进行微分，而不是在噪声数据上。

  • 语法 

    [df, d2f] = differentiate(fitresult, xd)

  • 语法说明

    • differentiate(fitresult, xd): 计算 fitresult (一个 cfit 对象) 在 xd 点处的一阶导数。[df, d2f] 则同时返回一阶和二阶导数。

  • 示例 

    x = linspace(0, 10, 101);
y = cos(x) + 0.1*randn(size(x)); % 带噪声的余弦

% 创建一个平滑样条拟合
f_fit = fit(x', y', 'smoothingspline', 'SmoothingParam', 0.1);

% 对拟合曲线求一阶导数
df = differentiate(f_fit, x);

figure;
plot(f_fit, x, y); % 绘制原始数据和拟合曲线
hold on;
plot(x, -sin(x), 'k--', 'LineWidth', 2, 'DisplayName', '真实导数');
plot(x, df, 'r-', 'LineWidth', 2, 'DisplayName', '拟合后求导');
legend('Location', 'northwest');
title('differentiate - 对含噪数据求导');

dlgradient (Deep Learning Toolbox)

  • 功能 实现深度学习中的核心功能——自动微分。它计算一个标量损失函数相对于一个或多个 dlarray 变量的梯度，是训练神经网络的基础。

  • 语法 

    [g1, g2, ...] = dlgradient(f, Y1, Y2, ...)

  • 语法说明

    • dlgradient(f, Y1, ...): f 是一个返回标量 dlarray 的函数句柄，Y1, ... 是 f 的输入参数，也是要求梯度的 dlarray 变量。g1, ... 是 f 的输出相对于 Y1, ... 的梯度。

  • 示例 

    % 定义一个简单的函数: f(w, b) = (w*x + b - y)^2
loss_fun = @(w, b) sum((w*2 + b - 5).^2);

% 初始化深度学习数组
w = dlarray(1);
b = dlarray(0);

% 计算梯度
[gW, gB] = dlgradient(loss_fun, w, b);

disp('损失函数对 w 的梯度:');
disp(extractdata(gW)); % 结果是 -12
disp('损失函数对 b 的梯度:');
disp(extractdata(gB)); % 结果是 -6