数値解析（03/XX）

はじめに

　プログラミング云々ばかり言っていても、それがどのように計算されているかを知らないと仕様が無い。ということで数値解析を学んでいく。
　基本的には

数値解析入門 (サイエンスライブラリ現代数学への入門)

を参照しつつ、他書で補足する。

前回の復習

• 反復法の基本的な構成手順を説明せよ。　適当な初期値$x_0$から出発し、求めたい解$\alpha$に収束するような数列$\{x_{\nu}\}$を構成し、$x_{\nu}$が$\alpha$に充分近づいたら計算を打ち切ってその値を近似解とする。　方程式$f(x)=0$を変形して
  $\begin{aligned}x=g(x)\end{aligned}$
  を得る。　次に求めたい解$\alpha$を含むような適当な閉区間$I$において$\varphi(x)\neq0$になるように選び、
  $\begin{aligned}g(x)=x-\varphi(x)f(x)\end{aligned}$
  とおく。次に$x_0$を適当に与えて
  $\begin{aligned}x_{\nu+1}=g(x_{\nu}),\ \nu=0,1,2,\cdots\end{aligned}$
  とおく。$\varphi(x)$を適当に与えることで様々な反復法を得る。
• 関数$f(x):I\rightarrow\mathbb{R},I\subset\mathbb{R}$が縮小写像であることの定義を述べよ。　関数$f(x)$が完備であるときに　(1) $x\in I\Rightarrow g(x)\in I$　(2) $\mathrm{Lipschitz}$条件　
  $\begin{aligned}x,x^{\prime}\in I\Rightarrow{}^{\exists}\lambda\in\left[0,1\right)\ \mathrm{s.t.}\ \left|g(x)-g(x^{\prime})\right|\leq\lambda\left|x-x^{\prime}\right|\end{aligned}$
  を満たす場合、$f(x)$は縮小写像であると呼ぶ。

前回

power-of-awareness.com

2.　非線形方程式の解

2.2　具体的な反復法

2.2.1　2分法

　反復法の一種ではあるがより簡単な方法である2分法をまずは扱う。
　2分法は中間値の定理を数学的な基礎とする。すなわち$f$が$[\min\{a,b\},\max\{a,b\}]$で連続であるとして、もし$f(a)\lt0,f(b)\gt0$となるような$a,b$が存在すれば、$f(\alpha)=0$となるような$\alpha$が$a,b$の間に少なくとも1つ存在する。これを基に以下のようなアルゴリズムを基に近似解を求める。

1. $f(a)\lt0,f(b)\gt0$を満たすようなとある$a,b$を$x$の値として取る。また収束判定用に$\varepsilon\gt0$を設定する。
2. 新たに
   $\begin{aligned}c=\displaystyle{\frac{a+b}{2}}\end{aligned}$
   を定義する。
3. もし$\displaystyle{\frac{\left|a-b\right|}{2}}\lt\varepsilon$ならば$5.$に進む。そうでなければ$4.$に進む。
4. もし$f(c)\gt0$ならば$b=c,$$f(c)\lt0$ならば$a=c$として$2.$に進む。$f(c)=0$ならば$5.$に進む。
5. $\alpha=c$とする。

　収束判定条件に着目すると、初めて収束判定条件を満たすような反復回数を$n$とおくと

$\begin{aligned} \displaystyle{\frac{\left|a-b\right|}{2^{n+1}}}\lt\varepsilon \end{aligned}$

が成り立つ。したがって$n\gt\log_2\left(\displaystyle{\frac{|a-b|}{\varepsilon}}\right)-1$を満たすような最小の$n$がこのアルゴリズムの計算量である。

2.2.2　Newton法

　二分法は安全確実ではあるものの計算に時間がかかる。より高速な方法が$\mathrm{Newton}$法である。

　微分可能な関数$y=f(x)$のグラフの解$\alpha$の付近の領域を考える。$x$の適当な初期値$x_0$を取る。次に点$\left(x_0,f\left(x_0\right)\right)$における$f(x)$の接線を引き、この接線と$x$軸の交点を$x_1$とおく。更に点$\left(x_1,f\left(x_1\right)\right)$における$f(x)$の接線を引き、この接線と$x$軸の交点を同様に$x_2$とおく。この操作を$n$回続けることで得た交点を$x_n$とすれば、点$\left(x_n,f\left(x_n\right)\right)$における$f(x)$の接線は

$\begin{aligned} y=f^{\prime}(x_n)(x-x_n)+f(x_n) \end{aligned}$

で与えられる。これと$y=0$との交点$x_{n+1}$を求めれば、

$\begin{aligned} x_{n+1}=x_n-\displaystyle{\frac{f(x_n)}{f^{\prime}(x_n)}} \end{aligned}$

が得られる。こうして定義できる数列$\left\{x_n\right\}$を用いることで近似解を求めることができる。
　いま求めたい方程式の解$\alpha$を含む適当な区間$\tilde{I}$において$f(x)$は$C^2$級かつ$f^{\prime}(x)\neq0$とすれば、

$\begin{aligned} g^{\prime}(x)&=1-\displaystyle{\frac{\left(f^{\prime}(x)\right)^2-f(x)f^{\prime\prime}(x)}{\left(f^{\prime}(x)\right)^2}}\\ &=\displaystyle{\frac{f(x)f^{\prime\prime}(x)}{\left(f^{\prime}(x)\right)^2}},\ g^{\prime}(\alpha)=0 \end{aligned}$

を得る。したがって$\alpha$を含む適当な区間$I\subset\tilde{I}$を取れば、系2の仮定が満たされるから、$\mathrm{Newton}$法は$f(x)=0$の解$\alpha$に収束する。

　$\mathrm{Newton}$法では現在の$x_n$の値が真の解に近づいているかを知る確実な方策が無い。そのため、ある正数$\varepsilon\gt0$を与え、

$\begin{aligned} \left|\displaystyle{\frac{x_{n+1}-x_{n}}{x_{n+1}}}\right|\lt\varepsilon \end{aligned}$

となった時点で計算を終了させ、$\alpha\approx x_{n+1}$とする。

　以上をまとめると、$\mathrm{Newton}$法は以下のようなアルゴリズムで定義できる：

1. 初期値$x$および$\varepsilon\gt0$を設定する。
2. $x_{\mathrm{new}}=x-\displaystyle{\frac{f(x)}{f^{\prime}(x)}}$とおく。
3. $\left|x_{\mathrm{new}}-x\right|\lt\varepsilon\left|x_{\mathrm{new}}\right|$ならば$5.$に移り、そうでなければ$4.$に移る。
4. $x=x_{\mathrm{new}}$とおき、$2.$に戻る。
5. $\alpha=x_{\mathrm{new}}$として計算を終了する。

　$\mathrm{Newton}$法の計算量は、2分法よりも少なく高速である。求めたい方程式の解を$\alpha$とし、$\varepsilon_n=x_n-\alpha$とおく。このとき$\varepsilon_{n}$と$\varepsilon_{n+1}$との関係を考えると、

$\begin{aligned} \varepsilon_{n+1}=x_{n+1}-\alpha=\left(x_n-\displaystyle{\frac{f(x_n)}{f^{\prime}(x_n)}}\right)-\alpha=\varepsilon_n-\displaystyle{\frac{f(x_n)}{f^{\prime}(x_n)}} \end{aligned}$

である。$x_n=\alpha+(x_n-\alpha)=\alpha+\varepsilon_n$と書き換えることができることから、$\mathrm{Taylor}$の定理より

$\begin{aligned} \begin{cases} f(x_n)&=f(\alpha)+\varepsilon_n f^{\prime}(\alpha)+\displaystyle{\frac{\varepsilon_n^2}{2!}}f^{\prime\prime}(\alpha)+\displaystyle{\frac{\varepsilon_n^3}{3!}}f^{\prime\prime\prime}(\xi_1)\\ f^{\prime}(x_n)&=f^{\prime}(\alpha)+\varepsilon_n f^{\prime\prime}(\alpha)+\displaystyle{\frac{\varepsilon_n^2}{2!}}f^{\prime\prime\prime}(\xi_2) \end{cases} \end{aligned}$

を満たすような$\xi_1,\xi_2\in\left[x_n,\alpha\right]$が存在する。これらを代入し、$f(\alpha)=0$を考慮すれば、

$\begin{aligned} \varepsilon_{n+1}&=\varepsilon_n-\displaystyle{\frac{f(x_n)}{f^{\prime}(x_n)}}\\ &=\varepsilon_n-\displaystyle{\frac{f(\alpha)+\varepsilon_n f^{\prime}(\alpha)+\displaystyle{\frac{\varepsilon_n^2}{2!}}f^{\prime\prime}(\alpha)+\displaystyle{\frac{\varepsilon_n^3}{3!}}f^{\prime\prime\prime}(\xi_1)}{f^{\prime}(\alpha)+\varepsilon_n f^{\prime\prime}(\alpha)+\displaystyle{\frac{\varepsilon_n^2}{2!}}f^{\prime\prime\prime}(\xi_2)}}\\ &=\varepsilon_n-\displaystyle{\frac{\varepsilon_n f^{\prime}(\alpha)+\displaystyle{\frac{\varepsilon_n^2}{2!}}f^{\prime\prime}(\alpha)+\displaystyle{\frac{\varepsilon_n^3}{3!}}f^{\prime\prime\prime}(\xi_1)}{f^{\prime}(\alpha)+\varepsilon_n f^{\prime\prime}(\alpha)+\displaystyle{\frac{\varepsilon_n^2}{2!}}f^{\prime\prime\prime}(\xi_2)}}\\ &=\displaystyle{\frac{\varepsilon_n^2}{2}}\cdot\displaystyle{\frac{f^{\prime\prime}(\alpha)-\displaystyle{\frac{\varepsilon_n}{3}}f^{\prime\prime\prime}(\xi_1)}{f^{\prime}(\alpha)+\varepsilon_n f^{\prime\prime}(\alpha)+\displaystyle{\frac{\varepsilon_n^2}{2!}}f^{\prime\prime\prime}(\xi_2)}}\\ &\approx\varepsilon_n^2\displaystyle{\frac{f^{\prime\prime}(\alpha)}{2f^{\prime}(\alpha)}} \end{aligned}$

を得る。これは$\mathrm{Newton}$法が急速に誤差を収束させていることを意味する（またこれは反復列$\{x_{\nu}\}$が2次収束することを意味する。）。
　他方で2分法に比べ、$\mathrm{Newton}$法は計算に失敗する確率が高い。これは与えられた$f(x)$と初期値$x_0$に由来する。

　前回示した反復法の収束条件を基に$\mathrm{Newton}$法の収束性を示す。

定理2.6　$\mathrm{Newton}$法の収束　区間$I\subset\mathbb{R}$で定義された$C^2$級関数$f(x)$を考える。区間$I$内に方程式$f(x)=0$には唯一の解$\alpha\in I$が存在し、$f^{\prime}(a)\neq0$が成り立つと仮定する。このとき$\mathrm{Newton}$法は初期値を$a$の近くから取る限り収束する。

($\because$　$\varphi(x)=x-\displaystyle{\frac{f(x)}{f^{\prime}(x)}}$とおいて$\left|\varphi^{\prime}(a)\right|\lt1$を示す。$f^{\prime}$の連続性および$f^{\prime}(a)\neq0$であるから、$f^{\prime}(x)\neq0,x\in J\subset I$を満たすような有界閉区間$J$が存在する。$\varphi$は$J$上では$C^{1}$級である。このとき

$\begin{aligned} \varphi^{\prime}(x)&=1-\displaystyle{\frac{\left(f^{\prime}(x)\right)^2-f(x)f^{\prime\prime}(x)}{\left(f^{\prime}(x)\right)^2}}\\ &=\displaystyle{\frac{f(x)f^{\prime\prime}(x)}{\left(f^{\prime}(x)\right)^2}},\ x\in J \end{aligned}$

が成り立つから、$\varphi^{\prime}(a)=0$が成り立ち、特に$\left|\varphi^{\prime}(a)\right|\lt1$　である。$\blacksquare$)

　なお$\mathrm{Newton}$法以外にも反復列を得る方法はいくつか存在する。
　反復する度に$f^{x_k}$を計算するのが難しければ、簡易$\mathrm{Newton}$法

$\begin{aligned} x_{k+1}=x_{k}-\displaystyle{\frac{f(x_k)}{f^{\prime}(x_0)}} \end{aligned}$

を用いる方法がある。また微分値を求めること自体が難しい場合は、

$\begin{aligned} f^{\prime}(x)=\displaystyle{\lim_{x\rightarrow c}\frac{f(x)-f(c)}{x-c}} \end{aligned}$

に注意して、セカント法

$\begin{aligned} x_{k+1}=x_{k}-\displaystyle{\frac{x_{k}-x_{k-1}}{f(x_{k})-f(x_{k-1})}}f(x_{k}) \end{aligned}$

が利用できる。1次収束も2次収束もしないものの、その中間くらいの収束性を示すことがある。ただしこの方法には$x_0$に加え$x_1$も初期値として用いる必要がある。

今回の復習

• 2分法の基礎となる定理を述べよ。
• $\mathrm{Newton}$法のアルゴリズムを説明せよ。

次回

power-of-awareness.com

参考文献

• 伊理正夫・藤野和建（1985）「数値計算の常識」（共立出版）
• 菊地文雄・齊藤宣一（2016）「数値解析の原理　現象の解明をめざして」（岩波書店）
• 齊藤宣一（2012）「数値解析入門」（東京大学出版）
• 齊藤宣一（2017）「数値解析」（共立出版）
• 高橋大輔（1996）「理工系の基礎数学8　数値計算」（岩波書店）
• 山本哲朗（2003）「数値解析入門[増訂版]」（サイエンス社）
• 日本応用数理学会 監修・ 櫻井 鉄也 編・ 松尾 宇泰 編・ 片桐 孝洋 編（2018）「数値線形代数の数理とHPC」（共立出版）
• Ridgway Scott, Larkin (2011), "Numerical Analysis", (Princeton University Press)