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

はじめに

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

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

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

Amazon

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

power-of-awareness.com

はじめに
2.　非線形方程式の解
- 2.2　具体的な反復法
  - 2.2.4　代数方程式の解法
- 2.3　連立非線形方程式
参考文献

2.　非線形方程式の解

2.2　具体的な反復法

2.2.4　代数方程式の解法

　実係数を持つ $n$ 次方程式

$\begin{aligned} P(z)=z^n+a_{1}z^{n-1}+\cdots+a_{n}=0,\ a_n\neq0 \end{aligned}$

を解く方法はいくつか知られている。原理的には何らかの方法、たとえば $\mathrm{Newton}$ 法で1つの解 $\alpha$ を求めた上で $z-\alpha$ もしくは $(z-\alpha)(z-\bar{\alpha})$ で上の方程式を割り、方程式の次数を下げる。この操作を繰り返すことで解ける。問題は初期値を如何にして見つけるかである。また割り算を繰り返すことで、丸め誤差により数値解が崩れることがある。
　ここではこうした難点の少ない $\mathrm{Durand}$ - $\mathrm{Kerner}$ - $\mathrm{Aberth}$ の方法を考える。

　まず $z=w-\displaystyle{\frac{a_1}{n}}$ とおき、

$\begin{aligned} P\left(w-\displaystyle{\frac{a_1}{n}}\right)&=\left(w-\displaystyle{\frac{a_1}{n}}\right)^n+a_{1}\left(w-\displaystyle{\frac{a_1}{n}}\right)^{n-1}+\cdots+a_{n}\\ &=w^n+c_2w^{n-2}+\cdots+c_n=0 \end{aligned}$

と変形する。このとき $(c_2,\cdots,c_n)\neq(0,\cdots,0),n\geq2$ ならば

$\begin{aligned} S(w)=w^n-|c_2|w^{n-2}-\cdots-|c_n|=0 \end{aligned}$

は1つの正解 $r$ を持ち、元の方程式 $P(z)=0$ の解 $\alpha,\cdots,\alpha_n$ はすべて閉円板

$\begin{aligned} \left|z+\displaystyle{\frac{a_1}{n}}\right|\leq r \end{aligned}$

に含まれる。そこで $r\leq r_0$ を満たすような $r_0$ を選び、

$\begin{aligned} z_k^{(0)}&=-\displaystyle{\frac{a_1}{n}}+r_0\exp\left[\sqrt{-1}\left(\displaystyle{\frac{2(k-1)\pi}{n}}+\displaystyle{\frac{\pi}{2n}}\right)\right],\ k=1,2,\cdots,n\\ z_k^{(\nu+1)}&=z_k^{(\nu)}-\displaystyle{\frac{P\left(z_k^{(\nu)}\right)}{\displaystyle{\prod_{j\neq k}z_k^{(\nu)}-z_j^{(\nu)}}}},\ 1\leq k\leq n,\nu=0,1,2,\cdots \end{aligned}$

とおく。もしある $j$ において $z_{k}^{(\nu)}=z_{j}^{(\nu)}$ となれば $z_{k}^{(\nu)}$ は収束したものとみなし、以後 $z_{k}^{(\nu+1)}=z_{k}^{(\nu)}$ とする。このようにして

$\begin{aligned} \left|z_{k}^{(\nu+1)}-z_{k}^{(\nu)}\right|\leq\varepsilon\left|z_{k}^{(\nu+1)}\right|,\ 1\leq k\leq n \end{aligned}$

となるまで反復する。

$\begin{aligned} z_k^{(1)}+\displaystyle{\frac{a_1}{n}}=\left(1-\displaystyle{\frac{1}{n}}\right)\left(z_k^{(0)}+\displaystyle{\frac{a_1}{n}}\right)+O\left(\left(z_k^{(0)}+\displaystyle{\frac{a_1}{n}}\right)^{-1}\right) \end{aligned}$

が知られており、したがって $r_0$ が充分に大きければ

$\begin{aligned} \left|z_k^{(1)}+\displaystyle{\frac{a_1}{n}}\right|\approx\left(1-\displaystyle{\frac{1}{n}}\right)r_0　 \end{aligned}$

となり、 $z_k^{(\nu)}$ は当初円の中心 $\left(-\displaystyle{\frac{a_1}{n}},0\right)$ に向かって1次収束するが、単解の場合、それらの付近で2次収束する。

定理2.9　 $\mathrm{Durand}$ - $\mathrm{Kerner}$ - $\mathrm{Aberth}$ の方法における解の精度　 $z_1,\cdots,z_n$ を相違する複素数とすれば、実係数を持つ方程式

$\begin{aligned} P(z)=z^n+a_{1}z^{n-1}+\cdots+a_{n}=0,\ a_n\neq0 \end{aligned}$

のすべての解は閉円板

$\begin{aligned} \mathit{\Gamma}_k:\left|z-z_k\right|\leq n\left| \displaystyle{\frac{P(z_k)}{\displaystyle{\prod_{j\neq k}(z_k-z_j)}}}\right|,\ 1\leq k\leq n \end{aligned}$

の合併に含まれる。その連結成分の1つが $m$ 個の閉円板から成り立つならば、その中に丁度 $m$ 個の解を持つ。

　適用対象において $z_k^{(\nu)}$ がすべて相違すれば $z_k=z_k^{(\nu)}$ として上の結果を適用すればよい。もしある $k,j$ において $z_k^{(\nu)}=z_j^{(\nu)}$ となるならば、充分に小さい $\delta\gt0$ を与えて $z_j=z_k^{(\nu)}+\delta$ などとして適用する。もし解の一部に重解や近接解があっても、それらと離れた解がかなり存在すれば、 $z_1,\cdots,z_n$ が良い近似のとき、 $\mathit{\Gamma}_k$ の半径は小さくなり、十分実用に耐え得る。もし $\mathit{\Gamma}_k$ が互いに素であれば、各 $\mathit{\Gamma}_k$ の半径が各近似解 $z_k$ の誤差限界を与える。

2.3　連立非線形方程式

　反復法は多変数の方程式にも有効である。基本的には単変数のときの考え方を援用すればよい。

2.3.1　例：2変数関数の場合

　まずは援用の仕方を学ぶべく、詳細な議論を避けて2変数の連立方程式

$\begin{aligned} \begin{cases} f(x,y)=0,\\ g(x,y)=0 \end{cases} \end{aligned}$

における $\mathrm{Newton}$ 法を考える。
　この数値解を考えるのに、解 $(a,b)$ に収束するような点列 $(x_1,y_1),\cdots,(x_n,y_n),\cdots$ を構成したい。反復が $(x_k,y_k)$ まで進んでいるとする。このとき $z=f(x,y)$ の $(x_k,y_k,f(x_k,y_k))$ における接平面および $z=g(x,y)$ の $z=f(x,y)$ の $(x_k,y_k,g(x_k,y_k))$ における接平面は、それぞれ

$\begin{aligned} z-f(x_k,y_k)&=(x-x_k)f_x(x_k,y_k)+(y-y_k)f_y(x_k,y_k),\\ z-g(x_k,y_k)&=(x-x_k)g_x(x_k,y_k)+(y-y_k)g_y(x_k,y_k) \end{aligned}$

で与えられる。ここで

$\begin{aligned} f_x&=\displaystyle{\frac{\partial f}{\partial x}},\\ f_y&=\displaystyle{\frac{\partial f}{\partial y}},\\ g_x&=\displaystyle{\frac{\partial g}{\partial x}},\\ g_y&=\displaystyle{\frac{\partial g}{\partial y}} \end{aligned}$

とおいた。
　これら接平面の式において $z=0$ とおいた

$\begin{aligned} -f(x_k,y_k)&=(x-x_k)f_x(x_k,y_k)+(y-y_k)f_y(x_k,y_k),\\ -g(x_k,y_k)&=(x-x_k)g_x(x_k,y_k)+(y-y_k)g_y(x_k,y_k) \end{aligned}$

は $xy$ 平面上の2本の直線を表し、それらの交点を $(x_{k+1},y_{k+1})$ とおく。このとき

$\begin{aligned} -\begin{bmatrix}f(x_k,y_k)\\g(x_k,y_k)\end{bmatrix}=\begin{bmatrix}f_x(x_k,y_k)&f_y(x_k,y_k)\\g_x(x_k,y_k)&g_y(x_k,y_k)\end{bmatrix}\begin{bmatrix}x_{k+1}-x_k\\y_{k+1}-y_k\end{bmatrix} \end{aligned}$

であるから、

$\begin{aligned} \begin{bmatrix}x_{k+1}\\y_{k+1}\end{bmatrix}=\begin{bmatrix}x_{k}\\y_{k}\end{bmatrix}-\begin{bmatrix}f_x(x_k,y_k)&f_y(x_k,y_k)\\g_x(x_k,y_k)&g_y(x_k,y_k)\end{bmatrix}^{-1}\begin{bmatrix}f(x_k,y_k)\\g(x_k,y_k)\end{bmatrix} \end{aligned}$

を得る。初期値 $(x_0,y_0)$ を適当に定め、この式を用いて反復列 $\left\{(x_k,y_k)\right\}_{k=1,2,\cdots}$ を生成する方法が2変数 $\mathrm{Newton}$ 法である。
　ここで現れた

$\begin{aligned} J_{f,g}(x,y)=\begin{bmatrix}f_x(x,y)&f_y(x,y)\\g_x(x,y)&g_y(x,y)\end{bmatrix} \end{aligned}$

を $f,g$ のヤコビ行列(ヤコビアン)という。
　ベクトル表示すると、 $\mathrm{Newton}$ 法は

$\begin{aligned} \boldsymbol{x}_{k+1}=\boldsymbol{x}_{k}-J_{\boldsymbol{f}}(\boldsymbol{x}_{k})^{-1}\boldsymbol{f}(\boldsymbol{x}_k) \end{aligned}$

で書ける。同様に2変数の簡易 $\mathrm{Newton}$ 法および緩和反復法は、緩和反復に適当な行列を $B$ としてそれぞれ

$\begin{aligned} \boldsymbol{x}_{k+1}&=\boldsymbol{x}_{k}-J_{\boldsymbol{f}}(\boldsymbol{x}_{0})^{-1}\boldsymbol{f}(\boldsymbol{x}_k),\\ \boldsymbol{x}_{k+1}&=\boldsymbol{x}_{k}-B\boldsymbol{f}(\boldsymbol{x}_k) \end{aligned}$

で与えられる。
　なお2変数の $\mathrm{Newton}$ 法は、複素数および複素関数を用いて表現することができる。 $z=x+iy,\ i=\sqrt{-1}$ とおき、正則な複素関数

$\begin{aligned} u(z)=f(x,y) + ig(x,y)=0 \end{aligned}$

を導入する。そして反復列 $\{z_{k}\}_{k=1,2,\cdots}$ を $z_k=x_k+iy_k$ として

$\begin{aligned} z_{k+1}=z_{k}-\displaystyle{\frac{u(z_k)}{u^{\prime}(z_k)}} \end{aligned}$

で定める。これを特に複素 $\mathrm{Newton}$ 法という。
　なお $u(z)$ が正則関数であるためには全微分可能かつ $\mathrm{Cauchy}$ - $\mathrm{Riemann}$ の微分方程式

$\begin{aligned} f_x=g_y,\ f_y=-g_x \end{aligned}$

を満たさなければならないため、複素 $\mathrm{Newton}$ 法の反復法が定義できなくとも、2変数の $\mathrm{Newton}$ 法が適用できる場合が存在する。

2.3.2　縮小写像の原理

　方程式

$\begin{aligned} \boldsymbol{x}&=\boldsymbol{g}(\boldsymbol{x}),\\ \boldsymbol{g}(\boldsymbol{x})&=\begin{bmatrix}g_1(\boldsymbol{x})\\\vdots\\g_n(\boldsymbol{x})\end{bmatrix} \end{aligned}$

を解く反復

$\begin{aligned} \boldsymbol{x}_{\nu+1}&=\boldsymbol{g}(\boldsymbol{x}_{\nu}) \end{aligned}$

を考える。このとき、1変数のときの縮小写像の原理を以下のように多変数に拡張することができる。証明は、1変数の場合のものにおいて、変数をベクトルに絶対値をノルムに置き換えて全く同じ議論をすればよいため省略する。

定理2.10　縮小写像の原理(多変数版)　 $n$ 次元 $\mathrm{Euclid}$ 空間 $\mathbb{R}^n$ 上のノルム付けされた閉集合 $\mathscr{D}$ で定義されたベクトル値関数 $\boldsymbol{g}(\boldsymbol{x})$ が条件

$\begin{aligned} \mathrm{(1)}\ \ \ \ &\boldsymbol{x}\in\mathscr{D}\Longrightarrow \boldsymbol{g}(\boldsymbol{x})\in\mathscr{D},\\ \mathrm{(2)}\ \ \ \ &\left\|\boldsymbol{g}(\boldsymbol{x})-\boldsymbol{g}(\boldsymbol{x}^{\prime})\right\|\leq\lambda\left\|\boldsymbol{x}-\boldsymbol{x}^{\prime}\right\|,\ \boldsymbol{x},\boldsymbol{x}^{\prime}\in\mathscr{D},\\ \mathrm{(3)}\ \ \ \ &\lambda\in \mathbb{R},\ 0\leq\lambda\lt1 \end{aligned}$

を満たすならば、方程式 $\boldsymbol{x}=\boldsymbol{g}(\boldsymbol{x})$ は一意の解 $\boldsymbol{\alpha}$ を持ち、それは反復

$\begin{aligned} \boldsymbol{x}_{\nu+1}&=\boldsymbol{g}(\boldsymbol{x}_{\nu}) \end{aligned}$

の極限として得られる、すなわち

$\begin{aligned} \boldsymbol{\alpha}=\displaystyle{\lim_{\nu\rightarrow\infty}\boldsymbol{x}_{\nu}} \end{aligned}$

が成り立つ。

　なお条件(2),(3)のみを満たすような $\boldsymbol{g}(\boldsymbol{x})$ を $\mathrm{Lipschitz}$ 連続であるという。

定理2.11　 $\mathrm{Lipschitz}$ 連続　閉集合 $\mathit{\Omega}\subset\mathbb{R}^n$ で定義された $\boldsymbol{g}(\boldsymbol{x})$ が、 $\lambda\gt0$ に対して、ノルム $\|\cdot\|$ が定義されていて、

$\begin{aligned} \left\|\boldsymbol{g}\left(\boldsymbol{x}\right)-\boldsymbol{g}\left(\boldsymbol{x}^{\prime}\right)\right\|\leq\lambda\left\|\boldsymbol{x}-\boldsymbol{x}^{\prime}\right\|,\ \boldsymbol{x},\boldsymbol{x}^{\prime}\in\mathit{\Omega} \end{aligned}$

が成り立つとき、 $\boldsymbol{g}(\boldsymbol{x})$ は（ $\|\cdot\|$ に関して） $\mathrm{Lipschitz}$ 連続であるという。特に $0\lt\lambda\lt1$ の場合を縮小写像という。

　縮小写像の原理から以下の系が得られる。

系2.1　縮小写像の原理　方程式 $\boldsymbol{x}=\boldsymbol{g}(\boldsymbol{x})$ の解 $\boldsymbol{\alpha}$ は、

$\begin{aligned} \mathscr{D}=\left\{\boldsymbol{x}\in\mathbb{R}^n\left|\right.\left\|\boldsymbol{x}-\boldsymbol{\alpha}\right\|\leq d\right\},\ d\gt0 \end{aligned}$

において、 $\boldsymbol{g}$ が縮小写像の原理における条件

$\begin{aligned} \mathrm{(2)}\ \ \ \ &\left\|\boldsymbol{g}(\boldsymbol{x})-\boldsymbol{g}(\boldsymbol{x}^{\prime})\right\|\leq\lambda\left\|\boldsymbol{x}-\boldsymbol{x}^{\prime}\right\|,\ \boldsymbol{x},\boldsymbol{x}^{\prime}\in\mathscr{D},\\ \mathrm{(3)}\ \ \ \ &\lambda\in \mathbb{R},\ 0\leq\lambda\lt1 \end{aligned}$

を満たすならば、

$\begin{aligned} \mathrm{(1)}\ \ \ \ &\boldsymbol{\alpha}\mathrm{のみが}\mathscr{D}\mathrm{内における解である},\\ \mathrm{(2)}\ \ \ \ &\boldsymbol{x}_{0}\in\mathscr{D},\boldsymbol{x}_{\nu+1}=\boldsymbol{g}(\boldsymbol{x}_{\nu})\Longrightarrow \boldsymbol{x}_{\nu}\in\mathscr{D}\land \boldsymbol{x}_{\nu}\rightarrow\boldsymbol{\alpha}(\nu\rightarrow\infty) \end{aligned}$

が成り立つ。

系2　縮小写像の原理　 $\boldsymbol{g}(\boldsymbol{x})$ は $\mathscr{D}=\left\{\boldsymbol{x}\left|\right.\left\|\boldsymbol{x}-\boldsymbol{\alpha}\right\|_{\infty}\leq d\right\},$ $d\gt0$ において $C^1$ 級で

$\begin{aligned} G(\boldsymbol{x})=\left(\displaystyle{\frac{\partial g_i(\boldsymbol{x})}{\partial x_i}}\right),\ \left\|G\left(\boldsymbol{x}\right)\right\|_{\infty}\leq\lambda\lt1,\ \boldsymbol{x}\in\mathscr{D} \end{aligned}$

であるとする。このとき方程式 $\boldsymbol{x}=\boldsymbol{g}(\boldsymbol{x})$ の解 $\boldsymbol{\alpha}$ について、

が成り立つ。

2.3.3　多変数の反復法

　次に一般の多変数方程式

$\begin{aligned} f_1(x_1,\cdots,x_n)&=0,\\ f_2(x_1,\cdots,x_n)&=0,\\ \vdots&\\ f_n(x_1,\cdots,x_n)&=0 \end{aligned}$

を考える。この場合も、

$\begin{aligned} \boldsymbol{x}&=\begin{bmatrix}x_1\\x_2\\\vdots\\x_n\end{bmatrix},\\ \boldsymbol{f}(\boldsymbol{x})&=\begin{bmatrix}f_1(\boldsymbol{x})\\f_2(\boldsymbol{x})\\\vdots\\f_n(\boldsymbol{x})\end{bmatrix},\\ J_{\boldsymbol{f}}(\boldsymbol{x})&=\begin{bmatrix} f_{1,1}(\boldsymbol{x})&\cdots&f_{1,n}(\boldsymbol{x})\\ \vdots&\ddots&\vdots\\ f_{n,1}(\boldsymbol{x})&\cdots&f_{n,n}(\boldsymbol{x}) \end{bmatrix},\\ f_{i,j}(\boldsymbol{x})&=\displaystyle{\frac{\partial f_i}{\partial x_j}}(x_1,\cdots,x_n) \end{aligned}$

とおけば、 $\mathrm{Newton}$ 法、簡易 $\mathrm{Newton}$ 法および緩和反復法は

$\begin{aligned} \boldsymbol{x}_{k+1}&=\boldsymbol{x}_{k}-J_{\boldsymbol{f}}(\boldsymbol{x}_{k})^{-1}\boldsymbol{f}(\boldsymbol{x}_k),\\ \boldsymbol{x}_{k+1}&=\boldsymbol{x}_{k}-J_{\boldsymbol{f}}(\boldsymbol{x}_{0})^{-1}\boldsymbol{f}(\boldsymbol{x}_k),\\ \boldsymbol{x}_{k+1}&=\boldsymbol{x}_{k}-B\boldsymbol{f}(\boldsymbol{x}_k) \end{aligned}$

で与えられる。

※今日の復習はなし

参考文献

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

はじめに

2. 非線形方程式の解

2.2 具体的な反復法

2.2.4 代数方程式の解法

2.3 連立非線形方程式

2.3.1 例：2変数関数の場合

2.3.2 縮小写像の原理

2.3.3 多変数の反復法

参考文献

2.　非線形方程式の解

2.2　具体的な反復法

2.2.4　代数方程式の解法

2.3　連立非線形方程式

2.3.1　例：2変数関数の場合

2.3.2　縮小写像の原理

2.3.3　多変数の反復法