数値解析（03/XX） - 「大人の教養・知識・気付き」を伸ばすブログ

はじめに

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

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

作者:山本哲朗
サイエンス社

Amazon

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

はじめに
前回の復習
前回
2.　非線形方程式の解
- 2.2　具体的な反復法
  - 2.2.1　2分法
  - 2.2.2　Newton法
今回の復習
次回
参考文献

前回の復習

反復法の基本的な構成手順を説明せよ。
　適当な初期値 $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つ存在する。これを基に以下のようなアルゴリズムを基に近似解を求める。

$f(a)\lt0,f(b)\gt0$ を満たすようなとある $a,b$ を $x$ の値として取る。また収束判定用に $\varepsilon\gt0$ を設定する。
新たに
$\begin{aligned}c=\displaystyle{\frac{a+b}{2}}\end{aligned}$
を定義する。
もし $\displaystyle{\frac{\left|a-b\right|}{2}}\lt\varepsilon$ ならば $5.$ に進む。そうでなければ $4.$ に進む。
もし $f(c)\gt0$ ならば $b=c,$ $f(c)\lt0$ ならば $a=c$ として $2.$ に進む。 $f(c)=0$ ならば $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}$ 法は以下のようなアルゴリズムで定義できる：

初期値 $x$ および $\varepsilon\gt0$ を設定する。
$x_{\mathrm{new}}=x-\displaystyle{\frac{f(x)}{f^{\prime}(x)}}$ とおく。
$\left|x_{\mathrm{new}}-x\right|\lt\varepsilon\left|x_{\mathrm{new}}\right|$ ならば $5.$ に移り、そうでなければ $4.$ に移る。
$x=x_{\mathrm{new}}$ とおき、 $2.$ に戻る。
$\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)

はじめに

前回の復習

前回

2. 非線形方程式の解

2.2 具体的な反復法

2.2.1 2分法

2.2.2 Newton法