偏导与差分

贡献者： addis

本文处于草稿阶段。

预备知识　多元泰勒展开，导数与差分

　　对于可微的二元函数 $f(x,y)$，可以证明

\begin{equation} \begin{aligned} \frac{\partial^2 }{\partial x \partial y} f(x_0,y_0) = \lim_{a\to 0}\lim_{b\to 0}\frac{1}{4ab}[-& f(x_0-a,\ y_0+b) + f(x_0+a,\ y_0+b)\\ & f(x_0-a,\ y_0-b) - f(x_0+a,\ y_0-b)]~, \end{aligned} \end{equation}

对每个 $f$ 使用多元泰勒展开即可证明。

Matlab 代码

　　计算过程中，由于出现了大数加小数，所以实际上 x(i) + h 增加的值可能比 h 有较大不同，可以用 (x(i) + h) - x(i) 来计算实际的增量，这样会使结果更精确。

　　注意该代码的有效数字估计只包括计算减法的截断误差，不包括用差分代替微分的误差。

代码 1：D2_ij.m

% 数值二阶偏导
% digi 是结果的有效数字
% D^2 f / D(x_i)D(x_j)
function [ret, digi] = D2_ij(f, i, j, x, h)
if i == j
    h = (x(i) + h) - x(i);
    if h == 0
        error('h 太小');
    end
    f2 = f(x);
    x(i) = x(i) + h;
    f3 = f(x);
    x(i) = x(i) - 2*h;
    f1 = f(x);
    val = f1 + f3;
    dif = val - 2*f2;
    ret = dif/h^2;
else % i ~= j
    x(i) = x(i) - 0.5*h; x(j) = x(j) - 0.5*h;
    hi = (x(i) + h) - x(i);
    hj = (x(j) + h) - x(j);
    if hi == 0 || hj == 0
        error('h 太小');
    end
    f1 = f(x);
    x(i) = x(i) + hi;
    f2 = f(x);
    x(j) = x(j) + hj;
    f3 = f(x);
    x(i) = x(i) - hi;
    f4 = f(x);
    val = f1+f3;
    dif = val-(f2+f4);
    ret = dif/(hi*hj);
end
if nargout == 2
    digi = max(0, 15.6536 + log10(abs(dif/val)));
end
end

　　另外给出一个变精度计算版本，可以使截断误差精确到最后一位有效数字。显然函数 f 需要支持变精度 sym 类型的变量。

代码 2：D2_ij_vpa.m

% 数值二阶偏导 (变精度)
% D^2 f / D(x_i)D(x_j)
% 精确到最后一位有效数字
function ret = D2_ij_vpa(f, i, j, x, h)
[ret, digi] = D2_ij(f, i, j, x, h);
if digi >= 15
    return;
end
old_digi = digits;
digits(2*15.6536 - digi - 8);
ret = D2_ij(f, i, j, arrayfun(@num2vpa, x), h);
ret = double(ret);
digits(old_digi);
end

致读者：小时百科一直以来坚持所有内容免费无广告，这导致我们处于严重的亏损状态。长此以往很可能会最终导致我们不得不选择大量广告以及内容付费等。因此，我们请求广大读者热心打赏 ，使网站得以健康发展。如果看到这条信息的每位读者能慷慨打赏 20 元，我们一周就能脱离亏损，并在接下来的一年里向所有读者继续免费提供优质内容。但遗憾的是只有不到 1% 的读者愿意捐款，他们的付出帮助了 99% 的读者免费获取知识，我们在此表示感谢。