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

===== 前回の課題 =====

[[ans5]]


===== 微分 =====

関数 $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を直線で結んだ傾きである。

{{ diff_forward.png?nolink |前進差分}}

例えば、$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}
\]
となる。

===== 小テスト① =====

[[KMS>|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)$ のグラフを描く。

<file rsplus diff1.R[enable_line_numbers="true"]>
# 数値微分（前進差分）

# 関数
f <- function(x) {
  return(x ^ 5 - 4 * x ^ 4)
}

# 前進差分
f1 <- function(x, h) {
  return((f(x + h) - f(x)) / h)
}

# 前進差分のグラフ
h <- 1e-1
curve(f1(x, h), xlim = c(-1, 3), ylab = "y")
</file>

続けて、次のコマンドで導関数 $f^{\prime}(x)=5x^{4}-16x^{3}$ を赤色の線で描いて前進差分の結果と比較する。
<code rsplus>
> curve(5 * x ^ 4 - 16 * x ^ 3, add = TRUE, col = "red")
</code>

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

=== $h=10^{-1}$ の場合 ===
<code rsplus>
h <- 1e-1
</code>

{{ diff1_h1e-1.png?nolink |h=1e-1 の場合}}

=== $h=10^{-10}$ の場合 ===
<code rsplus>
h <- 1e-10
</code>

{{ diff1_h1e-10.png?nolink |h=1e-10 の場合}}

=== $h=10^{-15}$ の場合 ===
<code rsplus>
h <- 1e-15
</code>

{{ diff1_h1e-15.png?nolink |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)|$ を計算する。

<file rsplus differr1.R[enable_line_numbers="true"]>
# 数値微分（前進差分）の誤差

# 関数
f <- function(x) {
  return(x ^ 5 - 4 * x ^ 4)
}

# 導関数（誤差確認用）
fp <- function(x) {
  return(5 * x ^ 4 - 16 * x ^ 3)
}

# 前進差分
f1 <- function(x, h) {
  return((f(x + h) - f(x)) / h)
}

x <- 0.1
h <- 10^(0 : -20)
error <- abs(fp(x) - f1(x, h))
range = c(1e-16, 1)
plot(h, error, ylim = range, log = "xy", type = "o", col = "blue")
</file>

{{ differr1.png?nolink | 前進差分の誤差}}

グラフは縦軸と横軸ともに対数表示で、例えば 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|
\]
となることから、**<fc #ff0000>$f_{1}(x,h)$ の打ち切り誤差は $h$ のオーダーで（$h$ に比例して）小さくなる</fc>**。

ところが $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)$ を__**中心差分**__という。

{{ diff_center.png?nolink |中心差分}}

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

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

<file rsplus diff2.R[enable_line_numbers="true"]>
# 数値微分（中心差分）

# 関数
f <- function(x) {
  return(x ^ 5 - 4 * x ^ 4)
}

# 中心差分
f2 <- function(x, h) {
  return((f(x + h) - f(x - h)) / (2 * h))
}

# 中心差分のグラフ
h <- 1e-1
curve(f2(x, h), xlim = c(-1, 3), ylab = "y")
</file>

続けて、次のコマンドで導関数 $f^{\prime}(x)=5x^{4}-16x^{3}$ を赤色の線で描いて中心差分の結果と比較する。
<code rsplus>
> curve(5 * x ^ 4 - 16 * x ^ 3, add = TRUE, col = "red")
</code>

=== $h=10^{-1}$ の場合 ===
{{ diff2_h1e-1.png?nolink |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|
\]
となることから、**<fc #ff0000>$f_{2}(x,h)$ の打ち切り誤差は $h^{2}$ のオーダーで（$h^{2}$ に比例して）小さくなる</fc>**。

===== 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)$ のグラフを描く。

<file rsplus diffr.R[enable_line_numbers="true"]>
# 数値微分（Romberg1段公式）

# 関数
f <- function(x) {
  return(x ^ 5 - 4 * x ^ 4)
}

# 中心差分
f2 <- function(x, h) {
  return((f(x + h) - f(x - h)) / (2 * h))
}

# Romberg1段公式
fr <- function(x, h) {
  return((4 * f2(x, h) - f2(x, 2 * h)) / 3)
}

# Romberg1段公式のグラフ
h <- 1e-1
curve(fr(x, h), xlim = c(-1, 3), ylab = "y")
</file>


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|
\]
となることから、**<fc #ff0000>$f_{r}(x,h)$ の打ち切り誤差は $h^{4}$ のオーダーで（$h^{4}$ に比例して）小さくなる</fc>**。

===== 誤差の比較 =====

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

<file rsplus differr2.R[enable_line_numbers="true"]>
# 数値微分（前進差分,中央差分,Romberg1段公式）の誤差

# 関数
f <- function(x) {
  return(x ^ 5 - 4 * x ^ 4)
}

# 導関数（誤差確認用）
fp <- function(x) {
  return(5 * x ^ 4 - 16 * x ^ 3)
}

# 前進差分
f1 <- function(x, h) {
  return((f(x + h) - f(x)) / h)
}

# 中心差分
f2 <- function(x, h) {
  return((f(x + h) - f(x - h)) / (2 * h))
}

# Romberg1段公式
fr <- function(x, h) {
  return((4 * f2(x, h) - f2(x, 2 * h)) / 3)
}

x <- 0.1
h <- 10^(0 : -20)
e1 <- abs(fp(x) - f1(x, h))
e2 <- abs(fp(x) - f2(x, h))
er <- abs(fp(x) - fr(x, h))
range <- c(1e-16, 1)
plot(h, e1, ylab = "error", ylim = range, log = "xy", type = "o", col = "blue")
par(new = TRUE)
plot(h, e2, ylim = range, log = "xy", type = "o", col = "red",
     ann = FALSE, axes = FALSE)
par(new = TRUE)
plot(h, er, ylim = range, log = "xy", type = "o", col = "green",
     ann = FALSE, axes = FALSE)
</file>
{{ differr2.png?nolink | 前進差分の誤差（青色）と中心差分の誤差（赤色）とRomberg1段公式の誤差（緑色）}}

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

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

===== 小テスト② =====

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

===== 定期試験 =====

  * 日時：2025年5月29日（木）1限
  * 場所：205教室
  * Moodle で行う
  * 持ち込み可

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

===== 授業アンケート =====

  * 科目名：数値情報処理a
  * 授業コード：1494

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