【数値情報処理a】第6回 数値微分

前回の課題

【数値情報処理a】第5回の答え

微分

関数 $f(x)$ の微分とは導関数 $f^{\prime}(x)$ を求めることである。

導関数 $f^{\prime}(x)$ の定義は以下になる。 \[ f^{\prime}(x)=\lim_{h\rightarrow0}f_{1}(x, h) \] ここで \[ f_{1}(x, h)=\frac{f(x+h)-f(x)}{h} \] である。 $f_{1}(x,h)$ を前進差分という。

$f_{1}(x,h)$ は以下の図の点Pと点Qを直線で結んだ傾きである。

前進差分

例えば、$f(x)=x^{2}$ の前進差分は \begin{eqnarray*} f_{1}(x,h) &=&\frac{(x+h)^{2}-x^{2}}{h}\\ &=&\frac{x^{2}+2xh+h^{2}-x^{2}}{h}\\ &=&2x+h \end{eqnarray*} であることから $f(x)=x^{2}$ の導関数は \[ f^{\prime}(x)=\lim_{h\rightarrow 0}(2x+h)=2x \] となる。

これを一般化すると $f(x)=x^{n}$ の前進差分は \begin{eqnarray*} f_{1}(x,h) &=&\frac{(x+h)^{n}-x^{n}}{h}\\ &=&nx^{n-1}+\frac{n(n-1)}{2}x^{n-2}h+\cdots \end{eqnarray*} である。 ここで二項定理 \begin{eqnarray*} (x+h)^{n}&=&_{n}C_{0}x^{n}h^{0}+_{n}C_{1}x^{n-1}h^{1}+_{n}C_{2}x^{n-2}h^{2}+\cdots+_{n}C_{n}x^{0}h^{n}\\ &=&x^{n}+nx^{n-1}h+\frac{n(n-1)}{2}x^{n-2}h^{2}+\cdots+h^{n} \end{eqnarray*} を使った。 したがって $f(x)=x^{n}$ の導関数は \[ f^{\prime}(x)=\lim_{h\rightarrow 0}\left\{nx^{n-1}+\frac{n(n-1)}{2}x^{n-2}h+\cdots\right\}=nx^{n-1} \] となる。

小テスト①

Moodle Server（非公式）で第6回の小テスト①を受験しなさい。

前進差分

前進差分 $f_{1}(x,h)$ の定義では分母に $h$ があるため、分子を変形して分母の $h$ を打ち消さないと $h\rightarrow 0$ の極限をとれない。 分子を変形して分母の $h$ を打ち消せればいいが、簡単にはできない場合もある。 そのような場合は $h$ をできるだけ小さい有限の値にとって $f_{1}(x, h)$ を計算することで $f^{\prime}(x)$ を近似的に求める方法がある。

例として次の関数を考えよう。 \[ f(x)=x^{5}-4x^{4} \] この導関数は \[ f^{\prime}(x)=5x^{4}-16x^{3} \] であることはすぐに求められる。 $h$ を有限にした $f_{1}(x,h)$ と $f^{\prime}(x)$ を比較してみる。

次のプログラム diff1.R では $f_{1}(x, h)$ のグラフを描く。

diff1.R
  2. # 数値微分（前進差分）
  3.  
  4. # 関数
  5. f <- function(x) {
  6.   return(x ^ 5 - 4 * x ^ 4)
  7. }
  8.  
  9. # 前進差分
  10. f1 <- function(x, h) {
  11.   return((f(x + h) - f(x)) / h)
  12. }
  13.  
  14. # 前進差分のグラフ
  15. h <- 1e-1
  16. curve(f1(x, h), xlim = c(-1, 3), ylab = "y")

続けて、次のコマンドで導関数 $f^{\prime}(x)=5x^{4}-16x^{3}$ を赤色の線で描いて前進差分の結果と比較する。 

> curve(5 * x ^ 4 - 16 * x ^ 3, add = TRUE, col = "red")

$h$ を変えて実行すると以下のようになる。

$h=10^{-1}$ の場合

h <- 1e-1

h=1e-1 の場合

$h=10^{-10}$ の場合

h <- 1e-10

h=1e-10 の場合

$h=10^{-15}$ の場合

h <- 1e-15

h=1e-15 の場合

$h$ が小さいほど $f_{1}(x,h)$ は $f^{\prime}(x)$ に近づくように思われるが、$h=10^{-15}$ になると逆に合わなくなっていることが分かる。 これはどうしてだろうか。

前進差分の誤差

$h$ を小さくしたときに前進差分の $f'(x)$ に対する誤差はどうなるだろうか。 $h$ を $10^{0}$ から $10^{-20}$ まで 0.1 倍ずつ小さくしていったときの誤差を見てみよう。

次のプログラム differr1.R では、$x=0.1$ として誤差の絶対値 $|f^{\prime}(x)-f_{1}(x,h)|$ を計算する。

differr1.R
  2. # 数値微分（前進差分）の誤差
  3.  
  4. # 関数
  5. f <- function(x) {
  6.   return(x ^ 5 - 4 * x ^ 4)
  7. }
  8.  
  9. # 導関数（誤差確認用）
  10. fp <- function(x) {
  11.   return(5 * x ^ 4 - 16 * x ^ 3)
  12. }
  13.  
  14. # 前進差分
  15. f1 <- function(x, h) {
  16.   return((f(x + h) - f(x)) / h)
  17. }
  18.  
  19. x <- 0.1
  20. h <- 10^(0 : -20)
  21. error <- abs(fp(x) - f1(x, h))
  22. range = c(1e-16, 1)
  23. plot(h, error, ylim = range, log = "xy", type = "o", col = "blue")

前進差分の誤差

グラフは縦軸と横軸ともに対数表示で、例えば 1e-05 は $10^{-5}$ のことである。 縦軸は誤差（error）で、小さくなるほど（グラフでは下になるほど）よい。

$h$ が $10^{0}$ から $10^{-9}$ まで小さくなるにしたがって誤差は小さくなる傾向にある（グラフの右上から中央まで）。 これは打ち切り誤差が小さくなることによる。

打ち切り誤差はテイラー展開によって評価できる。 テイラー展開とは、$f(x)$ が無限回微分可能であるとき、 $f(x+h)$ を $x$ の付近の $h$ の無限の多項式で展開することである。 \[ f(x+h)=f(x)+f^{\prime}(x)h+\frac{1}{2!}f^{\prime\prime}(x)h^{2}+\frac{1}{3!}f^{\prime\prime\prime}(x)h^{3}+\cdots \] ここで $f^{\prime}(x)$ は $f(x)$ の一階微分、$f^{\prime\prime}(x)$ は $f(x)$ の二階微分、$f^{\prime\prime\prime}(x)$ は $f(x)$ の三階微分などである。 これを使うと \[ f_{1}(x,h)=f^{\prime}(x)+\frac{1}{2!}f^{\prime\prime}(x)h+\frac{1}{3!}f^{\prime\prime\prime}(x)h^{2}+\cdots \] より \[ |f^{\prime}(x)-f_{1}(x,h)|=\left|\frac{1}{2!}f^{\prime\prime}(x)h+\frac{1}{3!}f^{\prime\prime\prime}(x)h^{2}+\cdots\right| \] となることから、$f_{1}(x,h)$ の打ち切り誤差は $h$ のオーダーで（$h$ に比例して）小さくなる。

ところが $h$ が $10^{-10}$ より小さくなると逆に誤差は大きくなっていく（グラフの中央から左上まで）。 これは $f_{1}(x,h)$ における $f(x+h)-f(x)$ の計算で桁落ちが発生するためである。

$h$ がどの程度で桁落ちが発生するかは丸め誤差による。 倍精度浮動小数点数の前進差分では $h$ は小さくても $h=10^{-10}$ 程度にするのが良さそうである。

中心差分

導関数 $f^{\prime}(x)$ は他に以下の定義もできる。 \[ f^{\prime}(x)=\lim_{h\rightarrow0}f_{2}(x, h) \] ここで \[ f_{2}(x, h)=\frac{f(x+h)-f(x-h)}{2h}=\frac{1}{2}f_{1}(x,h)+\frac{1}{2}f_{1}(x-h,h) \] である。 $f_{2}(x, h)$ を中心差分という。

中心差分

前進差分のときと同じように $h$ をできるだけ小さい有限の値にとって $f_{2}(x, h)$ を計算すれば $f^{\prime}(x)$ を近似できる。

次のプログラム diff2.R では $f_{2}(x, h)$ のグラフを描く。

diff2.R
  2. # 数値微分（中心差分）
  3.  
  4. # 関数
  5. f <- function(x) {
  6.   return(x ^ 5 - 4 * x ^ 4)
  7. }
  8.  
  9. # 中心差分
  10. f2 <- function(x, h) {
  11.   return((f(x + h) - f(x - h)) / (2 * h))
  12. }
  13.  
  14. # 中心差分のグラフ
  15. h <- 1e-1
  16. curve(f2(x, h), xlim = c(-1, 3), ylab = "y")

続けて、次のコマンドで導関数 $f^{\prime}(x)=5x^{4}-16x^{3}$ を赤色の線で描いて中心差分の結果と比較する。 

> curve(5 * x ^ 4 - 16 * x ^ 3, add = TRUE, col = "red")

$h=10^{-1}$ の場合

h=1e-1 の場合

中心差分では $h=10^{-1}$ の時点で $h_{2}(x,h)$ は十分 $f^{\prime}(x)$ に近いことが分かる。

中心差分の打ち切り誤差をテイラー展開によって評価する。 $f(x+h)$ と $f(x-h)$ のテイラー展開はそれぞれ \[ f(x+h)=f(x)+f^{\prime}(x)h+\frac{1}{2!}f^{\prime\prime}(x)h^{2}+\frac{1}{3!}f^{\prime\prime\prime}(x)h^{3}+\cdots \] \[ f(x-h)=f(x)-f^{\prime}(x)h+\frac{1}{2!}f^{\prime\prime}(x)h^{2}-\frac{1}{3!}f^{\prime\prime\prime}(x)h^{3}+\cdots \] となる。 これを使うと \[ f_{2}(x,h)=f^{\prime}(x)+\frac{1}{3!}f^{\prime\prime\prime}(x)h^{2}+\cdots \] より \[ |f^{\prime}(x)-f_{2}(x,h)|=\left|\frac{1}{3!}f^{\prime\prime\prime}(x)h^{2}+\cdots\right| \] となることから、$f_{2}(x,h)$ の打ち切り誤差は $h^{2}$ のオーダーで（$h^{2}$ に比例して）小さくなる。

Romberg1段公式

中心差分は次の公式を使って改良できる。 \[ f_{r}(x, h)=\frac{f_{2}(x,h)-\frac{1}{4}f_{2}(x,2h)}{1-\frac{1}{4}}=\frac{4f_{2}(x,h)-f_{2}(x,2h)}{3} \] これをRomberg1段公式という。

次のプログラム diffr.R では $f_{r}(x, h)$ のグラフを描く。

diffr.R
  2. # 数値微分（Romberg1段公式）
  3.  
  4. # 関数
  5. f <- function(x) {
  6.   return(x ^ 5 - 4 * x ^ 4)
  7. }
  8.  
  9. # 中心差分
  10. f2 <- function(x, h) {
  11.   return((f(x + h) - f(x - h)) / (2 * h))
  12. }
  13.  
  14. # Romberg1段公式
  15. fr <- function(x, h) {
  16.   return((4 * f2(x, h) - f2(x, 2 * h)) / 3)
  17. }
  18.  
  19. # Romberg1段公式のグラフ
  20. h <- 1e-1
  21. curve(fr(x, h), xlim = c(-1, 3), ylab = "y")

Romberg1段公式の打ち切り誤差をテイラー展開で評価する。 \[ f_{2}(x,h)=f^{\prime}(x)+\frac{1}{3!}f^{\prime\prime\prime}(x)h^{2}+\frac{1}{5!}f^{\prime\prime\prime\prime\prime}(x)h^{4}+\cdots \] \[ f_{2}(x,2h)=f^{\prime}(x)+\frac{2^{2}}{3!}f^{\prime\prime\prime}(x)h^{2}+\frac{2^{4}}{5!}f^{\prime\prime\prime\prime\prime}(x)h^{4}+\cdots \] これを使うと \[ f_{r}(x,h)=f^{\prime}(x)-\frac{4}{5!}f^{\prime\prime\prime\prime\prime}(x)h^{4}+\cdots \] より \[ |f^{\prime}(x)-f_{r}(x,h)|=\left|\frac{4}{5!}f^{\prime\prime\prime\prime\prime}(x)h^{4}+\cdots\right| \] となることから、$f_{r}(x,h)$ の打ち切り誤差は $h^{4}$ のオーダーで（$h^{4}$ に比例して）小さくなる。

誤差の比較

次のプログラム differr2.R は $x=0.1$ で前進差分、中心差分、Romberg1段公式のそれぞれの誤差を計算して比較する。

differr2.R
  2. # 数値微分（前進差分,中央差分,Romberg1段公式）の誤差
  3.  
  4. # 関数
  5. f <- function(x) {
  6.   return(x ^ 5 - 4 * x ^ 4)
  7. }
  8.  
  9. # 導関数（誤差確認用）
  10. fp <- function(x) {
  11.   return(5 * x ^ 4 - 16 * x ^ 3)
  12. }
  13.  
  14. # 前進差分
  15. f1 <- function(x, h) {
  16.   return((f(x + h) - f(x)) / h)
  17. }
  18.  
  19. # 中心差分
  20. f2 <- function(x, h) {
  21.   return((f(x + h) - f(x - h)) / (2 * h))
  22. }
  23.  
  24. # Romberg1段公式
  25. fr <- function(x, h) {
  26.   return((4 * f2(x, h) - f2(x, 2 * h)) / 3)
  27. }
  28.  
  29. x <- 0.1
  30. h <- 10^(0 : -20)
  31. e1 <- abs(fp(x) - f1(x, h))
  32. e2 <- abs(fp(x) - f2(x, h))
  33. er <- abs(fp(x) - fr(x, h))
  34. range <- c(1e-16, 1)
  35. plot(h, e1, ylab = "error", ylim = range, log = "xy", type = "o", col = "blue")
  36. par(new = TRUE)
  37. plot(h, e2, ylim = range, log = "xy", type = "o", col = "red",
  38.      ann = FALSE, axes = FALSE)
  39. par(new = TRUE)
  40. plot(h, er, ylim = range, log = "xy", type = "o", col = "green",
  41.      ann = FALSE, axes = FALSE)

前進差分の誤差（青色）と中心差分の誤差（赤色）とRomberg1段公式の誤差（緑色）

グラフの青色は前進差分の誤差、赤色は中心差分の誤差、緑色は Romberg1段公式の誤差を示す。 中心差分は前進差分より速く誤差が小さくなり、Romberg1段公式はさらに速く誤差が小さくなっていることが分かる。

三つとも $h$ が小さくなると打ち切り誤差が見られるようになるが、$h$ を適切に選ぶことで桁落ちの影響が出てくる前に打ち切り誤差を小さくできる。

小テスト②

Moodle Server（非公式）で第6回の小テスト②を受験しなさい。

定期試験

今までの小テストと課題に似た問題を出すので復習しておくように。

授業アンケート

「er-学籍番号8桁@k.koeki-u.ac.jp」に空メールを送信すると返事のメールが届くので、後はその指示に従うこと。