「大人の教養・知識・気付き」を伸ばすブログ

一流の大人（ビジネスマン、政治家、リーダー…）として知っておきたい、教養・社会動向を意外なところから取り上げ学ぶことで“気付く力”を伸ばすブログです。データ分析・語学に力点を置いています。　→諸事情により、2023年4月～8月くらいまで更新を不定期にします(8/23追記：10月末まで更新不定期期間を延期します。)。

MENU

数値解析（05/XX）

数値解析

はじめに

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

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

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

Amazon

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

power-of-awareness.com

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}について、

\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.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}

で与えられる。

※今日の復習はなし

参考文献

プライバシーポリシー お問い合わせ