导数与差分

贡献者： addis

预备知识　泰勒展开

1. 一阶导数

　　我们在导数的定义中已经知道¹

\begin{matrix} (1) & f^{'} (x) = lim_{h \to 0} \frac{f (x + h / 2) - f (x - h / 2)}{h} . \end{matrix}

在一些应用（如数值计算）中，我们只能把

h

取一个很小的数值（如

10^{- 10}

）而并非无穷小，这就需要我们估计用上式右边的差分来代替

f^{'} (x)

有多精确。为了估算误差，我们可以将

f (x \pm h / 2)

展开为关于

h

的泰勒级数

\begin{matrix} (2) & f (x \pm h / 2) = f (x) \pm f^{'} (x) \frac{h}{2} + \frac{1}{2} f^{″} (x) {(\frac{h}{2})}^{2} + O (h^{3}) . \end{matrix}

代入式 1 得

\begin{matrix} (3) & lim_{h \to 0} \frac{f^{'} (x) h + O (h^{3})}{h} = f^{'} (x) + O (h^{2}), \end{matrix}

所以用差分代替一阶导数可以精确到

h

的二阶无穷小

O (h^{2})

。

2. 二阶导数

　　能否用类似的方法来表示二阶导数呢？根据二阶导数的定义，我们需要用双重极限来表示

\begin{matrix} (4) & \begin{aligned} f^{″} (x) & = lim_{l \to 0} \frac{f^{'} (x + l / 2) - f^{'} (x - l / 2)}{l} \\ = lim_{l \to 0} lim_{h \to 0} \frac{1}{l h} [f (x + l / 2 + h / 2) - f (x + l / 2 - h / 2) \\ - f (x - l / 2 + h / 2) + f (x - l / 2 - h / 2)], \end{aligned} \end{matrix}

但我们希望只用一个极限来表示二阶导数。然而我们不确定

h

是否需要是

l

的高阶无穷小。我们不妨来试试令

l = h

，即

\begin{matrix} (5) & f^{″} (x) = lim_{h \to 0} \frac{f (x + h) - 2 f (x) + f (x - h)}{h^{2}} . \end{matrix}

要验证该式成立与否，将

f (x \pm h)

关于

h

做泰勒展开得

\begin{matrix} (6) & f (x \pm h) = f (x) \pm f^{'} (x) h + \frac{1}{2} f^{″} (x) h^{2} \pm \frac{1}{6} f^{‴} (x) h^{3} + O (h^{4}) . \end{matrix}

代入式 5 右边得

\begin{matrix} (7) & lim_{h \to 0} \frac{f^{″} (x) h^{2} + O (h^{4})}{h^{2}} = f^{″} (x) + O (h^{2}) . \end{matrix}

这就验证了式 5 的正确性。另外我们得知用差分来近似二阶导数

f^{″} (x)

同样是精确到二阶无穷小

O (h^{2})

。

Matlab 代码

代码 1：D_i.m

% 数值偏偏导
% digi 是结果的有效数字
function [ret, digi] = D_i(f, i, x, h)
x(i) = x(i) - 0.5*h;
f1 = f(x);
h = (x(i) + h) - x(i);
if h == 0
    error('h 太小');
end
x(i) = x(i) + h;
f2 = f(x);
dif = f2 - f1;
ret = dif/h;
if nargout == 2
    digi = max(0, 15.6536 + log10(abs(dif/f1)));
end
end

代码 2：D_i_vpa.m

% 数值偏导 (变精度)
% f 需要支持 sym 类型的变精度参数
% 精确到最后一位有效数字
function ret = D_i_vpa(f, i, x, h)
[ret, digi] = D_i(f, i, x, h);
if digi >= 15
    return;
end
old_digi = digits;
digits(2*15.6536 - digi - 8);
ret = D_i(f, i, arrayfun(@num2vpa,x), h);
ret = double(ret);
digits(old_digi);
end

1. ^ 以下假设 $f (x)$ 在某区间内处处可导。

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