贡献者： addis
1. 一阶导数
我们在导数的定义中已经知道1
\begin{equation}
f'(x) = \lim_{h\to 0}\frac{f(x + h/2) - f(x - h/2)}{h}~.
\end{equation}
在一些应用（如数值计算）中，我们只能把 $h$ 取一个很小的数值（如 $10^{-10}$）而并非无穷小，这就需要我们估计用上式右边的
差分来代替 $f'(x)$ 有多精确。为了估算误差，我们可以将 $f(x \pm h/2)$ 展开为关于 $h$ 的泰勒级数
\begin{equation}
f(x \pm h/2) = f(x) \pm f'(x)\frac h2 + \frac12 f''(x) \left(\frac h2 \right) ^2 + \mathcal{O}\left(h^3 \right) ~.
\end{equation}
代入
式 1 得
\begin{equation}
\lim_{h\to 0} \frac{f'(x)h + \mathcal{O}\left(h^3 \right) }{h} = f'(x) + \mathcal{O}\left(h^2 \right) ~,
\end{equation}
所以用差分代替一阶导数可以精确到 $h$ 的二阶无穷小 $ \mathcal{O}\left(h^2 \right) $。
2. 二阶导数
能否用类似的方法来表示二阶导数呢？根据二阶导数的定义，我们需要用双重极限来表示
\begin{equation} \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}{lh} [f(x + l/2 + h/2) - f(x + l/2 - h/2)\\
&\qquad\qquad - f(x - l/2 + h/2) + f(x - l/2 - h/2)]~,
\end{aligned} \end{equation}
但我们希望只用一个极限来表示二阶导数。然而我们不确定 $h$ 是否需要是 $l$ 的高阶无穷小。我们不妨来试试令 $l = h$，即
\begin{equation}
f''(x) = \lim_{h\to 0} \frac{f(x + h) - 2f(x) +f(x-h)}{h^2}~.
\end{equation}
要验证该式成立与否，将 $f(x \pm h)$ 关于 $h$ 做泰勒展开得
\begin{equation}
f(x \pm h) = f(x) \pm f'(x) h + \frac12 f''(x) h^2 \pm \frac16 f'''(x) h^3 + \mathcal{O}\left(h^4 \right) ~.
\end{equation}
代入
式 5 右边得
\begin{equation}
\lim_{h\to 0} \frac{f''(x)h^2 + \mathcal{O}\left(h^4 \right) }{h^2} = f''(x) + \mathcal{O}\left(h^2 \right) ~.
\end{equation}
这就验证了
式 5 的正确性。另外我们得知用差分来近似二阶导数 $f''(x)$ 同样是精确到二阶无穷小 $ \mathcal{O}\left(h^2 \right) $。
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)$ 在某区间内处处可导。