数値微分(高校数学)

同義語:差分近似(高校数学)numerical differentiation (high school mathematics)

概要

数値微分(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$ を小さくしすぎると誤差はかえって大きくなりうる。

$$\newcommand{C}[0]{\mathbb{C}} \newcommand{div}[0]{\mathbin{÷}} \newcommand{N}[0]{\mathbb{N}} \newcommand{Q}[0]{\mathbb{Q}} \newcommand{R}[0]{\mathbb{R}} \newcommand{Z}[0]{\mathbb{Z}} $$

前提知識: 1次近似と誤差の伝わり方(高校数学), 微分係数と導関数の定義(高校数学), 関数の近似とテイラー展開

高校での出発点:極限をとる前で止める

微分係数は、平均変化率の極限
$$ f'(a)=\lim_{h\to0}\frac{f(a+h)-f(a)}{h} $$
として定義される(微分係数と導関数の定義(高校数学))。関数の式が分からず、いくつかの点での値しか分からないとき(実験のデータや、計算機のプログラムが出す値)は、極限をとることができない。そこで、$h$ を小さな数に止めて平均変化率を計算し、それを $f'(a)$ の近似値とする。これを 数値微分 という。
同じ $h$ でも、平均変化率のとり方で精度が大きく変わる。手で計算できる例で比べる。

$x^2$ と $x^3$ の $x=1$ での微分係数を見積もる

$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$
  1. の誤差は $h$ に比例して減り、(2) の誤差は $h^2$ に比例して減る。$x^2$ では (2) がちょうど正しい。

この記事で答える問いは次の 3 つである。

  1. 右側だけの差と左右対称の差で、誤差はどれだけ違うか。その上界を保証できるか。→ thm-ndf-truncation
  2. $h$ はどこまでも小さくすればよいか。→ thm-ndf-rounding
  3. 差をとる方法がうまくいかないのはどんなときか。→ ex-ndf-cx-abs、ex-ndf-cx-osc
    高校の計算この記事の言葉大学の言葉
    平均変化率 $\dfrac{f(a+h)-f(a)}{h}$前進差分差分商
    極限をとる前の $h$刻み幅刻み幅
    平均変化率と $f'(a)$ の差打ち切り誤差Taylor の定理の剰余項
    計算機の値の小さな誤差丸め誤差浮動小数点数の誤差、桁落ち

3 つの差分

以下、$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)$ の値を使わない。

$e^x$ の $x=0$ での微分係数

$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 を結ぶ)を描き、中心差分の直線の傾きが接線の傾きに近いことを見る図 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))$ を通らないが、接線とほぼ平行である。

主定理 1:前進差分は $h$、中心差分は $h^2$ に比例する誤差

差分と $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 $$
である。

1 次近似の誤差と、2 次の項の打ち消し
  1. 1次近似と誤差の伝わり方(高校数学) の定理「1次近似の誤差の上界」により
    $$ \bigl\lvert f(a+h)-f(a)-f'(a)h\bigr\rvert\le\frac{M_2}2h^2 $$
    である。左辺の絶対値の中身を $h$ で割ると
    $$ \frac{f(a+h)-f(a)-f'(a)h}{h}=\frac{f(a+h)-f(a)}h-f'(a)=D_+(h)-f'(a) $$
    なので、両辺を $h>0$ で割って $\bigl\lvert D_+(h)-f'(a)\bigr\rvert\le\dfrac{M_2}2h$ を得る。後退差分は、同じ定理を $h$ の代わりに $-h$ で使うと $\bigl\lvert f(a-h)-f(a)+f'(a)h\bigr\rvert\le\dfrac{M_2}2h^2$ で、$\dfrac{f(a-h)-f(a)+f'(a)h}{-h}=D_-(h)-f'(a)$ だから、同様に示される。
  2. 方針:$f(a+h)$ と $f(a-h)$ を 2 次の式と剰余項で書き、引き算で $f(a)$ と 2 次の項が消えることを使う。
    段 1(2 次の近似式)。Taylor の定理の Lagrange の形(関数の近似とテイラー展開 の定理「Lagrange の形の剰余項」を $n=2$ で使う。Leb26a Theorem 4.3.2 も同じ)により、$a$ と $a+h$ の間の点 $c_1$ と、$a-h$ と $a$ の間の点 $c_2$ で
    $$ f(a+h)=f(a)+f'(a)h+\frac{f''(a)}2h^2+\frac{f'''(c_1)}6h^3,\qquad f(a-h)=f(a)-f'(a)h+\frac{f''(a)}2h^2-\frac{f'''(c_2)}6h^3 $$
    となるものがある(2 つ目は $x=a-h$ とした式で、$(-h)^2=h^2$、$(-h)^3=-h^3$ を使った)。
    段 2(引き算)。1 つ目から 2 つ目を引くと、$f(a)$ と $\dfrac{f''(a)}2h^2$ が打ち消し合う。
    $$ f(a+h)-f(a-h)=\cancel{f(a)-f(a)}+2f'(a)h+\cancel{\frac{f''(a)}2h^2-\frac{f''(a)}2h^2}+\frac{f'''(c_1)+f'''(c_2)}6h^3 $$
    段 3($2h$ で割る)。両辺を $2h$ で割ると
    $$ D_0(h)-f'(a)=\frac{f'''(c_1)+f'''(c_2)}{12}h^2 $$
    である。$c_1$、$c_2$ はどちらも $a-h$ と $a+h$ の間にあるので $\lvert f'''(c_1)\rvert,\lvert f'''(c_2)\rvert\le M_3$ で、右辺の大きさは $\dfrac{2M_3}{12}h^2=\dfrac{M_3}6h^2$ 以下である。$\square$

前進差分では、$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$ と等しい。

$\sin x$ の $x=1$ での打ち切り誤差

$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$ で微分係数を求めるときは、左の点がないので中心差分は使えず、前進差分を使う。

主定理 2:$h$ を小さくしすぎると悪くなる

thm-ndf-truncation だけを見ると、$h$ はいくらでも小さくすればよさそうである。ところが、計算機で実際に計算すると、そうはならない。

計算機で $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$ は小さすぎても大きすぎてもいけない」。

最適な $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 乗あたりから丸め誤差の線に沿って増えることを見る図 横軸に h(右へ行くほど小さい)、縦軸に誤差の大きさをとった両対数のグラフで、前進差分と中心差分の誤差が h を小さくすると一度減り、h が 10 の −8 乗や −5 乗あたりから丸め誤差の線に沿って増えることを見る図
図 2 は、ex-ndf-double の誤差を細かい $h$ で描いたものである。右へ行くほど $h$ が小さい。左側では誤差は打ち切りの上界(点線)に沿って減り、右側では丸めの上界(破線)の程度に増える。2 つの線の交わるあたりで誤差が最も小さい。丸め誤差による部分はでこぼこしているが、これは個々の値の丸め方が $h$ ごとにばらばらだからである。
同じ考え方で、2 階の微分係数も差分で見積もれる。

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$ 以下である。


$f(x)=\sin x$、$a=1$ では $f''(1)=-\sin1$、$M_4=1$ で、$h=0.1$ の誤差は $7.01\times10^{-4}$(上界 $8.33\times10^{-4}$)、$h=0.01$ では $7.01\times10^{-6}$(上界 $8.33\times10^{-6}$)である。ただし、$h^2$ で割るので丸め誤差の影響はもっと強い。値の誤差を $\delta$ とすると、丸めの部分の上界は $\dfrac{4\delta}{h^2}$ になる。倍精度で計算すると、誤差の大きさは $h=10^{-4}$ で $3.0\times10^{-9}$、$h=10^{-5}$ で $7.6\times10^{-7}$、$h=10^{-6}$ で $7.8\times10^{-5}$ と、$h=10^{-4}$ より先では $h$ を小さくするほど悪くなる。

例と反例

差分による見積もりには、次の前提がある。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=|x| の原点で、前進差分の直線は傾き 1、後退差分の直線は傾き −1、中心差分の直線は水平で傾き 0 になることを見る図
y=sin 50x の原点で、接線の傾きは 50 だが、h=0.1 の中心差分の直線は右下がりで傾きが約 −9.59 になることを見る図 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)。前進差分と後退差分が食い違うことが、微分できないことの手がかりになる。

反例:$h$ に比べて速く変わる関数

$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$ は、関数が変化する細かさよりずっと小さくとる必要がある。

大学数学で見る:Richardson の補外

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$ の項を打ち消している。


計算機の丸め誤差の一般的な扱い(浮動小数点数のモデル、桁落ち、条件数)は 誤差(数値解析) で扱う。同記事の節「打ち切り誤差」にあるとおり、刻み幅を小さくすると打ち切り誤差は減るが丸め誤差は増えうるので、全体の釣り合いを見る必要がある。thm-ndf-rounding は、その最も簡単な場合である。

さらに先へ

  • 補間多項式を微分する:3 点 $(a-h,f(a-h))$、$(a,f(a))$、$(a+h,f(a+h))$ を通る 2 次以下の多項式を $x=a$ で微分すると、ちょうど中心差分 $D_0(h)$ になり、2 回微分すると $D_2(h)$ になる。点の値を通る多項式とその誤差は、次の記事 Lagrange補間(高校数学) で扱う。
  • 差分の記号:数列の差 $a_{n+1}-a_n$ を記号で扱う方法は 数列の和と差分 にある。$h=1$ の前進差分は数列の階差である。
  • 微分方程式への応用:微分方程式 $y'=F(x,y)$ で、微分係数 $y'$ を前進差分で置き換えると Euler法(高校数学) になる。
  • 関数の式が分かっているとき:微分の公式で導関数を求めれば、丸め誤差の問題は起こらない。数値微分は、値しか分からないときの道具である。

関連項目

参考文献

Mathpediaは寄付と、参考文献の書籍リンク(Amazonアソシエイト)の紹介料で運営されています。 支援について / 寄付する