数値微分(numerical differentiation)とは、微分係数 $f'(a)$ を前進差分 $\frac{f(a+h)-f(a)}{h}$ や中心差分 $\frac{f(a+h)-f(a-h)}{2h}$($h>0$)で近似することである。使う点の間で $f''$ が連続で $\lvert f''\rvert\le M_2$ なら前進差分の誤差は $\frac{M_2}{2}h$ 以下、$f'''$ が連続で $\lvert f'''\rvert\le M_3$ なら中心差分の誤差は $\frac{M_3}{6}h^2$ 以下である。関数の値に $\delta$ 以下の誤差があると中心差分の誤差の上界は $\frac{M_3}{6}h^2+\frac{\delta}{h}$ で、$h$ を小さくしすぎると誤差はかえって大きくなりうる。
微分係数は、平均変化率の極限
$$
f'(a)=\lim_{h\to0}\frac{f(a+h)-f(a)}{h}
$$
として定義される(微分係数と導関数の定義(高校数学))。関数の式が分からず、いくつかの点での値しか分からないとき(実験のデータや、計算機のプログラムが出す値)は、極限をとることができない。そこで、$h$ を小さな数に止めて平均変化率を計算し、それを $f'(a)$ の近似値とする。これを 数値微分 という。
同じ $h$ でも、平均変化率のとり方で精度が大きく変わる。手で計算できる例で比べる。
$f'(1)$ の真の値は、$f(x)=x^2$ で $2$、$f(x)=x^3$ で $3$ である。
(1) 右側の点だけを使う:$\dfrac{f(1+h)-f(1)}{h}$。$x^2$ では $\dfrac{(1+h)^2-1}{h}=\dfrac{2h+h^2}{h}=2+h$、$x^3$ では $\dfrac{(1+h)^3-1}{h}=\dfrac{3h+3h^2+h^3}{h}=3+3h+h^2$ である。
(2) 左右の点を対称に使う:$\dfrac{f(1+h)-f(1-h)}{2h}$。$x^2$ では $\dfrac{(1+h)^2-(1-h)^2}{2h}=\dfrac{4h}{2h}=2$、$x^3$ では $\dfrac{(1+h)^3-(1-h)^3}{2h}=\dfrac{6h+2h^3}{2h}=3+h^2$ である。
| $h=0.1$ | $h=0.01$ | 真の値 | |
|---|---|---|---|
| $x^2$、(1) 右側だけ | $2.1$ | $2.01$ | $2$ |
| $x^2$、(2) 左右対称 | $2$ | $2$ | $2$ |
| $x^3$、(1) 右側だけ | $3.31$ | $3.0301$ | $3$ |
| $x^3$、(2) 左右対称 | $3.01$ | $3.0001$ | $3$ |
この記事で答える問いは次の 3 つである。
| 高校の計算 | この記事の言葉 | 大学の言葉 |
|---|---|---|
| 平均変化率 $\dfrac{f(a+h)-f(a)}{h}$ | 前進差分 | 差分商 |
| 極限をとる前の $h$ | 刻み幅 | 刻み幅 |
| 平均変化率と $f'(a)$ の差 | 打ち切り誤差 | Taylor の定理の剰余項 |
| 計算機の値の小さな誤差 | 丸め誤差 | 浮動小数点数の誤差、桁落ち |
以下、$h>0$ とする。
$f'(a)$ の近似値として使う次の 3 つの式を、それぞれ 前進差分、後退差分、中心差分 という。
$$
D_+(h)=\frac{f(a+h)-f(a)}{h},\qquad D_-(h)=\frac{f(a)-f(a-h)}{h},\qquad D_0(h)=\frac{f(a+h)-f(a-h)}{2h}
$$
前進差分は点 $a$ と右の点 $a+h$、後退差分は左の点 $a-h$ と点 $a$ を結ぶ直線の傾き、中心差分は左右の点 $a-h$、$a+h$ を結ぶ直線の傾きである。
中心差分は、前進差分と後退差分の平均である。実際
$$
\frac{D_+(h)+D_-(h)}2=\frac{\bigl(f(a+h)-f(a)\bigr)+\bigl(f(a)-f(a-h)\bigr)}{2h}=\frac{f(a+h)-f(a-h)}{2h}=D_0(h)
$$
である。中心差分は $f(a)$ の値を使わない。
$f(x)=e^x$、$a=0$、$h=0.1$ とする。$f'(0)=1$ である。
$$
D_+(0.1)=\frac{e^{0.1}-1}{0.1}=1.05171\ldots,\qquad D_-(0.1)=\frac{1-e^{-0.1}}{0.1}=0.95163\ldots,\qquad D_0(0.1)=\frac{e^{0.1}-e^{-0.1}}{0.2}=1.00167\ldots
$$
である。前進差分は $0.05$ ほど大きく、後退差分は $0.05$ ほど小さく、平均の中心差分では誤差が $0.0017$ まで下がる。$e^x$ は下に凸なので、右の点を結ぶ直線は接線より急で、左の点を結ぶ直線はゆるやかになる。2 つの誤差がほぼ打ち消し合うのが中心差分である。
y=sin x の x=a での接線と、前進差分の直線(点 a と a+h を結ぶ)、中心差分の直線(点 a−h と a+h を結ぶ)を描き、中心差分の直線の傾きが接線の傾きに近いことを見る図
図 1 は $y=\sin x$、$a=1$、$h=0.5$ で 3 本の直線を比べたものである。接線の傾きは $\cos1=0.5403\ldots$ で、前進差分の直線の傾きは $0.3120\ldots$、中心差分の直線の傾きは $0.5181\ldots$ である。中心差分の直線は、点 $(a,f(a))$ を通らないが、接線とほぼ平行である。
差分と $f'(a)$ の差を 打ち切り誤差 という(極限をとらずに途中で打ち切ったことによる誤差)。
$h>0$ とする。
(1) $f$ が $a$ と $a+h$ を含む区間で 2 回微分可能で $f''$ が連続とし、$a\le t\le a+h$ で $\lvert f''(t)\rvert\le M_2$ とする。このとき
$$
\bigl\lvert D_+(h)-f'(a)\bigr\rvert\le\frac{M_2}2h
$$
である。後退差分についても、$a-h\le t\le a$ で $\lvert f''(t)\rvert\le M_2$ なら $\bigl\lvert D_-(h)-f'(a)\bigr\rvert\le\dfrac{M_2}2h$ である。
(2) $f$ が $a-h$ と $a+h$ を含む区間で 3 回微分可能で $f'''$ が連続とし、$a-h\le t\le a+h$ で $\lvert f'''(t)\rvert\le M_3$ とする。このとき
$$
\bigl\lvert D_0(h)-f'(a)\bigr\rvert\le\frac{M_3}6h^2
$$
である。
前進差分では、$f(a+h)$ の 2 次の項 $\dfrac{f''(a)}2h^2$ が残り、$h$ で割ると $\dfrac{f''(a)}2h$ の誤差になる。中心差分では、左右対称にとったので 2 次の項が打ち消し合い、次の 3 次の項が誤差になる。これが問い 1 の答えである。
ex-ndf-start で確かめる。$x^2$ では $f'''=0$ なので $M_3=0$ で、中心差分の誤差は $0$、つまりちょうど正しい。$x^3$ では $f'''=6$ なので上界は $\dfrac66h^2=h^2$ で、実際の誤差 $h^2$ と等しい。前進差分では、$x^2$ で $f''=2$ なので上界は $\dfrac22h=h$ で、実際の誤差 $h$ と等しい。
$f(x)=\sin x$、$a=1$ とする。$f'(1)=\cos1=0.5403023\ldots$ である。$\lvert f'''(t)\rvert=\lvert\cos t\rvert\le1$ なので $M_3=1$ とでき、中心差分の上界は $\dfrac{h^2}6$ である。前進差分では $f''(t)=-\sin t$ で、$1\le t\le1+h$($h\le0.1$)では $\lvert f''(t)\rvert\le\sin(1+h)$ なので、$M_2=\sin(1+h)$ とする。誤差(差分 $-$ 真の値)は次のとおりである(打ち切り誤差だけを見るため、40 桁で計算した)。
| $h$ | 前進差分の誤差 | 上界 $\dfrac{\sin(1+h)}2h$ | 中心差分の誤差 | 上界 $\dfrac{h^2}6$ |
|---|---|---|---|---|
| $10^{-1}$ | $-4.29\times10^{-2}$ | $4.46\times10^{-2}$ | $-9.001\times10^{-4}$ | $1.67\times10^{-3}$ |
| $10^{-2}$ | $-4.22\times10^{-3}$ | $4.23\times10^{-3}$ | $-9.005\times10^{-6}$ | $1.67\times10^{-5}$ |
| $10^{-4}$ | $-4.21\times10^{-5}$ | $4.21\times10^{-5}$ | $-9.005\times10^{-10}$ | $1.67\times10^{-9}$ |
$h$ を $\dfrac1{10}$ にすると、前進差分の誤差は約 $\dfrac1{10}$、中心差分の誤差は約 $\dfrac1{100}$ になる。前進差分の誤差を $h$ で割ると $-\dfrac{\sin1}2=-0.4207\ldots$ に、中心差分の誤差を $h^2$ で割ると $-\dfrac{\cos1}6=-0.09005\ldots$ に近づく。どちらも上界の中にある。
3 つの差分を比べる。
| 前進差分 $D_+$ | 後退差分 $D_-$ | 中心差分 $D_0$ | |
|---|---|---|---|
| 使う点 | $a$、$a+h$ | $a-h$、$a$ | $a-h$、$a+h$ |
| 誤差の上界 | $\dfrac{M_2}2h$ | $\dfrac{M_2}2h$ | $\dfrac{M_3}6h^2$ |
| $h$ を $\dfrac1{10}$ にすると誤差は | 約 $\dfrac1{10}$ | 約 $\dfrac1{10}$ | 約 $\dfrac1{100}$ |
| ちょうど正しい関数 | 1 次以下の多項式 | 1 次以下の多項式 | 2 次以下の多項式 |
| 必要な滑らかさ | $f''$ が連続 | $f''$ が連続 | $f'''$ が連続 |
区間の端 $a$ で微分係数を求めるときは、左の点がないので中心差分は使えず、前進差分を使う。
thm-ndf-truncation だけを見ると、$h$ はいくらでも小さくすればよさそうである。ところが、計算機で実際に計算すると、そうはならない。
$f(x)=\sin x$、$a=1$ の前進差分と中心差分を、倍精度の浮動小数点数(有効数字およそ 16 桁で計算する、計算機の標準的な数の表し方)で計算すると、誤差(差分 $-\cos1$)は次のようになる。
| $h$ | 前進差分の誤差 | 中心差分の誤差 |
|---|---|---|
| $10^{-1}$ | $-4.29\times10^{-2}$ | $-9.00\times10^{-4}$ |
| $10^{-2}$ | $-4.22\times10^{-3}$ | $-9.00\times10^{-6}$ |
| $10^{-4}$ | $-4.21\times10^{-5}$ | $-9.00\times10^{-10}$ |
| $10^{-5}$ | $-4.21\times10^{-6}$ | $-1.11\times10^{-11}$ |
| $10^{-6}$ | $-4.21\times10^{-7}$ | $2.77\times10^{-11}$ |
| $10^{-8}$ | $-2.97\times10^{-9}$ | $2.58\times10^{-9}$ |
| $10^{-12}$ | $4.32\times10^{-5}$ | $-1.23\times10^{-5}$ |
中心差分は $h=10^{-4}$ までは ex-ndf-sin と同じく $h^2$ に比例して減るが、$h=10^{-5}$ あたりで最も小さくなり、その先は $h$ を小さくするほど悪くなる。$h=10^{-12}$ の誤差は、$h=10^{-2}$ のときより大きい。
原因は、計算機が出す $\sin(1+h)$ などの値そのものに、ごく小さな誤差(丸め誤差)があることである。$h$ が小さいと $f(a+h)$ と $f(a-h)$ は近い 2 数になり、その差をとると、値の誤差が相対的に大きく残る。1次近似と誤差の伝わり方(高校数学) の例「反例:近い 2 数の差」で見た桁落ちである。さらに、その差を小さな $2h$ で割るので、誤差は $\dfrac1{2h}$ 倍に拡大される。
この現象を式で押さえる。計算機が出す値を $F(x)$ と書き、真の値 $f(x)$ との差が $\delta$ 以下であるとする。ここでは関数の値の誤差だけを考え、引き算・割り算や $a\pm h$ の計算での丸めは考えない。
thm-ndf-truncation (2) の仮定のもとで、2 つの数 $F(a+h)$、$F(a-h)$ が $\bigl\lvert F(a\pm h)-f(a\pm h)\bigr\rvert\le\delta$ を満たすとする($\delta>0$)。$\tilde D_0(h)=\dfrac{F(a+h)-F(a-h)}{2h}$ とおくと
$$
\bigl\lvert\tilde D_0(h)-f'(a)\bigr\rvert\le\underbrace{\frac{M_3}6h^2}_{\text{打ち切り誤差}}+\underbrace{\frac\delta h}_{\text{丸め誤差}}
$$
が成り立つ。さらに $M_3>0$ のとき、右辺を $h>0$ の関数と見ると、$h=h_*=\left(\dfrac{3\delta}{M_3}\right)^{\frac13}$ で最小になり、最小値は $\dfrac{M_3}2h_*^{\,2}$ である。
段 1(2 つに分ける)。$D_0(h)$ を真の値での中心差分とすると
$$
\tilde D_0(h)-f'(a)=\bigl(\tilde D_0(h)-D_0(h)\bigr)+\bigl(D_0(h)-f'(a)\bigr)
$$
で、三角不等式から $\bigl\lvert\tilde D_0(h)-f'(a)\bigr\rvert\le\bigl\lvert\tilde D_0(h)-D_0(h)\bigr\rvert+\bigl\lvert D_0(h)-f'(a)\bigr\rvert$ である。
段 2(丸めの部分)。
$$
\tilde D_0(h)-D_0(h)=\frac{\bigl(F(a+h)-f(a+h)\bigr)-\bigl(F(a-h)-f(a-h)\bigr)}{2h}
$$
で、分子の大きさは $\delta+\delta=2\delta$ 以下なので、$\bigl\lvert\tilde D_0(h)-D_0(h)\bigr\rvert\le\dfrac{2\delta}{2h}=\dfrac\delta h$ である。
段 3(打ち切りの部分)。thm-ndf-truncation (2) により $\bigl\lvert D_0(h)-f'(a)\bigr\rvert\le\dfrac{M_3}6h^2$ である。段 1〜3 で不等式が示された。
段 4(最小値)。$g(h)=\dfrac{M_3}6h^2+\dfrac\delta h$($h>0$)とおくと
$$
g'(h)=\frac{M_3}3h-\frac\delta{h^2}=\frac{M_3h^3-3\delta}{3h^2}
$$
である。分子 $M_3h^3-3\delta$ は $h$ について増加し、$h^3=\dfrac{3\delta}{M_3}$ すなわち $h=h_*$ で $0$ になる。よって $0< h< h_*$ で $g'(h)<0$、$h>h_*$ で $g'(h)>0$ で、$g$ は $h=h_*$ で最小になる。$h_*^{\,3}=\dfrac{3\delta}{M_3}$ より $\dfrac\delta{h_*}=\dfrac{M_3h_*^{\,3}}{3h_*}=\dfrac{M_3}3h_*^{\,2}$ なので、最小値は $g(h_*)=\dfrac{M_3}6h_*^{\,2}+\dfrac{M_3}3h_*^{\,2}=\dfrac{M_3}2h_*^{\,2}$ である。$\square$
打ち切り誤差は $h$ を小さくすると減り、丸め誤差の上界 $\dfrac\delta h$ は $h$ を小さくすると増える。両者の釣り合う所で全体の上界が最小になる。これが問い 2 の答えで、「$h$ は小さすぎても大きすぎてもいけない」。
倍精度では、$1$ 程度の大きさの数は、相対的におよそ $1.1\times10^{-16}$ の精度で表される。$\sin(1\pm h)$ の値が正しく丸められているとして、$\delta=1.1\times10^{-16}$ とする。$\sin$ では $M_3=1$ なので
$$
h_*=\left(3\times1.1\times10^{-16}\right)^{\frac13}=6.9\times10^{-6},\qquad\frac{M_3}2h_*^{\,2}=2.4\times10^{-11}
$$
である。ex-ndf-double で、中心差分の誤差が最も小さかったのは $h=10^{-5}$ の $1.1\times10^{-11}$ で、見積もりと合っている。どんなに $h$ を選んでも、中心差分で得られる $\cos1$ の精度は、上界で言えば $10^{-11}$ 程度(有効数字 10 桁ほど)までである。
前進差分では、同じように $\bigl\lvert\tilde D_+(h)-f'(a)\bigr\rvert\le\dfrac{M_2}2h+\dfrac{2\delta}h$ となる(丸めの部分は $\dfrac{\delta+\delta}h$)。右辺は $h=2\sqrt{\dfrac\delta{M_2}}$ で最小値 $2\sqrt{M_2\delta}$ をとる(相加平均と相乗平均の関係 $\dfrac{M_2}2h+\dfrac{2\delta}h\ge2\sqrt{M_2\delta}$ からも分かる)。$M_2=1$、$\delta=1.1\times10^{-16}$ では $h=2.1\times10^{-8}$、最小値 $2.1\times10^{-8}$ で、ex-ndf-double の $h=10^{-8}$ の誤差の大きさ $3.0\times10^{-9}$ は、$h=10^{-8}$ での上界 $\dfrac{10^{-8}}2+\dfrac{2.2\times10^{-16}}{10^{-8}}=2.7\times10^{-8}$ 以下である。最も良い $h$ での上界で比べると、中心差分($2.4\times10^{-11}$)は前進差分($2.1\times10^{-8}$)より 900 倍ほど小さい。
横軸に h(右へ行くほど小さい)、縦軸に誤差の大きさをとった両対数のグラフで、前進差分と中心差分の誤差が h を小さくすると一度減り、h が 10 の −8 乗や −5 乗あたりから丸め誤差の線に沿って増えることを見る図
図 2 は、ex-ndf-double の誤差を細かい $h$ で描いたものである。右へ行くほど $h$ が小さい。左側では誤差は打ち切りの上界(点線)に沿って減り、右側では丸めの上界(破線)の程度に増える。2 つの線の交わるあたりで誤差が最も小さい。丸め誤差による部分はでこぼこしているが、これは個々の値の丸め方が $h$ ごとにばらばらだからである。
同じ考え方で、2 階の微分係数も差分で見積もれる。
$f$ が $a-h$ と $a+h$ を含む区間で 4 回微分可能で $f^{(4)}$ が連続とし、$a-h\le t\le a+h$ で $\lvert f^{(4)}(t)\rvert\le M_4$ とする。このとき
$$
D_2(h)=\frac{f(a+h)-2f(a)+f(a-h)}{h^2},\qquad\bigl\lvert D_2(h)-f''(a)\bigr\rvert\le\frac{M_4}{12}h^2
$$
である。$D_2(h)$ は、前進差分と後退差分の差を $h$ で割ったもの $\dfrac{D_+(h)-D_-(h)}h$ に等しい。
prf-ndf-truncation の段 1 と同じく、Taylor の定理の Lagrange の形を $n=3$ で使うと、$a$ と $a+h$ の間の $c_1$、$a-h$ と $a$ の間の $c_2$ で
$$f(a\pm h)=f(a)\pm f'(a)h+\frac{f''(a)}2h^2\pm\frac{f'''(a)}6h^3+\frac{f^{(4)}(c_{1,2})}{24}h^4$$
(複号同順、$+$ のとき $c_1$、$-$ のとき $c_2$)となる。2 つを足すと $f'(a)h$ と $\dfrac{f'''(a)}6h^3$ の項が打ち消し合い、$f(a+h)+f(a-h)=2f(a)+f''(a)h^2+\dfrac{f^{(4)}(c_1)+f^{(4)}(c_2)}{24}h^4$ である。$2f(a)$ を移して $h^2$ で割ると $D_2(h)-f''(a)=\dfrac{f^{(4)}(c_1)+f^{(4)}(c_2)}{24}h^2$ で、大きさは $\dfrac{2M_4}{24}h^2=\dfrac{M_4}{12}h^2$ 以下である。
差分による見積もりには、次の前提がある。1 つを外すと崩れる。
| 外す条件 | 反例 | 成り立たなくなること |
|---|---|---|
| $f$ が $x=a$ で微分可能 | $\lvert x\rvert$ の $x=0$ | 中心差分の値が微分係数の近似になる |
| $\dfrac{M_3}6h^2$ が小さい($h$ に比べて $f$ がゆっくり変わる) | $\sin50x$ の $x=0$、$h=0.1$ | 中心差分が $f'(a)$ に近い |
| $h$ が小さすぎない(丸め誤差が小さい) | $\sin x$ の $x=1$、$h=10^{-12}$ | $h$ を小さくすると誤差が減る(ex-ndf-double) |
y=|x| の原点で、前進差分の直線は傾き 1、後退差分の直線は傾き −1、中心差分の直線は水平で傾き 0 になることを見る図
y=sin 50x の原点で、接線の傾きは 50 だが、h=0.1 の中心差分の直線は右下がりで傾きが約 −9.59 になることを見る図
$f(x)=\lvert x\rvert$、$a=0$ とする。どんな $h>0$ でも
$$
D_+(h)=\frac{\lvert h\rvert-0}{h}=1,\qquad D_-(h)=\frac{0-\lvert -h\rvert}{h}=-1,\qquad D_0(h)=\frac{\lvert h\rvert-\lvert -h\rvert}{2h}=0
$$
である。中心差分は $h$ によらず $0$ という落ち着いた値を出すが、$\lvert x\rvert$ は $x=0$ で微分できない(右微分係数 $1$ と左微分係数 $-1$ が違う。微分係数と導関数の定義(高校数学) の反例「左右で傾きが違う(絶対値)」)。中心差分の値が $h$ を変えても一定だからといって、それが微分係数だとは限らない(図 3)。前進差分と後退差分が食い違うことが、微分できないことの手がかりになる。
$f(x)=\sin50x$、$a=0$ とする。$f'(x)=50\cos50x$ なので $f'(0)=50$ である。中心差分は
$$
D_0(h)=\frac{\sin50h-\sin(-50h)}{2h}=\frac{\sin50h}{h}
$$
で、$h=0.1$ では $\dfrac{\sin5}{0.1}=-9.589\ldots$ となり、符号まで違う(図 4)。$h=0.01$ では $47.94\ldots$(誤差 $-2.06$)、$h=0.001$ では $49.979\ldots$(誤差 $-0.0208$)である。
thm-ndf-truncation (2) は成り立っている。$f'''(x)=-50^3\cos50x$ なので $M_3=125000$ で、$h=0.1$ の上界は $\dfrac{125000}6\times0.01=208.3\ldots$ である。実際の誤差 $59.6$ はこれ以下だが、上界は何の役にも立たない。使う 2 点 $x=-0.1$ と $x=0.1$ の間で、$\sin50x$ は $\dfrac{50\times0.2}{2\pi}=1.6$ 周期ほども振動するので、2 点だけで傾きを読むのは無理である。$h$ は、関数が変化する細かさよりずっと小さくとる必要がある。
thm-ndf-truncation (2) の証明をもう 1 次先まで進めると、$f$ が十分なめらかなら、中心差分の誤差は $D_0(h)-f'(a)=Ch^2+(h^4\text{ の程度})$($C=\dfrac{f'''(a)}6$)と書ける。そこで、$h$ と $\dfrac h2$ の 2 つの中心差分を
$$
\frac{4D_0\left(\frac h2\right)-D_0(h)}3
$$
と組み合わせると、$4\cdot C\left(\dfrac h2\right)^2-Ch^2=0$ なので $h^2$ の項が消え、誤差は $h^4$ の程度になる。この方法を Richardson の補外 という。
$\sin x$ の $x=1$ で、$D_0(0.1)$ の誤差は $-9.0005\times10^{-4}$、$D_0(0.05)$ の誤差は $-2.2510\times10^{-4}$(約 $\dfrac14$)である。$\dfrac{4D_0(0.05)-D_0(0.1)}3$ の誤差は $-1.125\times10^{-7}$ で、同じ 4 つの値から、$D_0(0.05)$ よりおよそ 2000 倍正確な値が得られる(打ち切り誤差だけを見るため、40 桁で計算した)。
この組み合わせ方は、台形公式とSimpsonの公式(高校数学) の定理「Simpson の公式の 2 つの性質」の $S_{2m}=\dfrac{4T_{2m}-T_m}3$ と同じ形である。誤差が $h^2$ に比例する近似を 2 つ組み合わせて、$h^2$ の項を打ち消している。
Mathpediaは寄付と、参考文献の書籍リンク(Amazonアソシエイト)の紹介料で運営されています。 支援について / 寄付する