导数与差分

                     

贡献者: 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)$ 在某区间内处处可导。

                     

© 小时科技 保留一切权利