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

一流の大人（ビジネスマン、政治家、リーダー…）として知っておきたい、教養・社会動向を意外なところから取り上げ学ぶことで“気付く力”を伸ばすブログです。現在、コンサルタントの雛になるべく、少しずつ勉強中です(※2024年12月10日改訂)。

MENU

数値解析（04/XX）

数値解析

はじめに

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

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

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

Amazon

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

前回の復習

  • 2分法の基礎となる定理を述べよ。

    　中間値の定理
  • \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}}として計算を終了する。

前回

power-of-awareness.com

2.　非線形方程式の解

2.2　具体的な反復法

2.2.3　Aitkenの加速法とSteffensenの反復法

　一般に\alphaに収束する反復列\{x_{\nu}\}が

\begin{aligned} x_{\nu+1}-\alpha=(\kappa-\varepsilon_{\nu})(x_{\nu}-\alpha),\ 0\lt|\kappa|\lt1,\ \displaystyle{\lim_{\nu\rightarrow\infty}\varepsilon_{\nu}}=0 \end{aligned}

を満たすとき、線形収束(1次収束)するという。
　反復列\{x_{\nu}\}が\alphaに線形収束するとすれば、0\lt|\kappa|\lt1を満たすような定数\kappaおよび0に収束する数列\{\varepsilon_{\nu}\}が存在し、

\begin{aligned} \begin{cases} x_{\nu+1}-\alpha&=(\kappa+\varepsilon_{\nu})(x_{\nu}-\alpha)\approx\kappa(x_{\nu}-\alpha)\\ x_{\nu+2}-\alpha&=(\kappa+\varepsilon_{\nu+1})(x_{\nu+1}-\alpha)\approx\kappa(x_{\nu+1}-\alpha) \end{cases} \end{aligned}

が成り立つ。これらから\kappaを消去すると、

\begin{aligned} &(x_{\nu+2}-\alpha)(x_{\nu}-\alpha)\approx (x_{\nu+1}-\alpha)^2\\ \Leftrightarrow&x_{\nu+2}x_{\nu}-(x_{\nu+2}+x_{\nu})\alpha\approx x_{\nu+1}^2-2x_{\nu+1}\alpha\\ \Leftrightarrow&(x_{\nu+2}+x_{\nu}-2x_{\nu+1})\alpha\approx x_{\nu+2}x_{\nu}-x_{\nu+1}^2\\ \Leftrightarrow&\alpha\approx\displaystyle{\frac{x_{\nu+2}x_{\nu}-x_{\nu+1}^2}{x_{\nu+2}+x_{\nu}-2x_{\nu+1}}}\\ \Leftrightarrow&\alpha\approx x_{\nu}-\displaystyle{\frac{(x_{\nu+1}-x_{\nu})^2}{x_{\nu+2}+x_{\nu}-2x_{\nu+1}}} \end{aligned}

を得る。そこで数列\{y_{\nu}\}を

\begin{aligned} y_{\nu}=x_{\nu}-\displaystyle{\frac{(x_{\nu+1}-x_{\nu})^2}{x_{\nu+2}-2x_{\nu+1}+x_{\nu}}}=x_{\nu+2}-\displaystyle{\frac{(x_{\nu+2}-x_{\nu+1})^2}{x_{\nu+2}-2x_{\nu+1}+x_{\nu}}} \end{aligned}

と定義すれば、\{y_{\nu}\}は\{x_{\nu}\}よりも速く\alphaに近づくと予想される。この\{y_{\nu}\}を用いた反復法を\mathrm{Aitken}の加速法という。この方法は、\{x_{\nu}\}の収束が遅い場合に用いる。一般に収束列からより収束の早い収束列をつくる方法を加速法という。



定理2.7　\mathrm{Aitken}の加速法の収束速度　反復列\{x_{\nu}\}が線形収束ならば、\mathrm{Aitken}加速列\{y_{\nu}\}

\begin{aligned} y_{\nu}=x_{\nu}-\displaystyle{\frac{(x_{\nu+1}-x_{\nu})^2}{x_{\nu+2}-2x_{\nu+1}+x_{\nu}}}=x_{\nu+2}-\displaystyle{\frac{(x_{\nu+2}-x_{\nu+1})^2}{x_{\nu+2}-2x_{\nu+1}+x_{\nu}}} \end{aligned}

は\{x_{\nu}\}を加速する。

(\because　e_{\nu}=x_{\nu}-\alphaとおけば、

\begin{aligned} x_{\nu+1}-\alpha&=(\kappa+\varepsilon_{\nu})(x_{\nu}-\alpha)\approx\kappa(x_{\nu}-\alpha),\\ x_{\nu+2}-\alpha&=(\kappa+\varepsilon_{\nu+1})(x_{\nu+1}-\alpha)\approx\kappa(x_{\nu+1}-\alpha) \end{aligned}

より、

\begin{aligned} e_{\nu+1}&=(\kappa+\varepsilon_{\nu})e_{\nu},\\ e_{\nu+2}&=(\kappa+\varepsilon_{\nu+1})(\kappa+\varepsilon_{\nu})e_{\nu} \end{aligned}

が成り立つ。\{y_{\nu}\}の定義式から、

\begin{aligned} y_{\nu}-\alpha&=e_{\nu}-\displaystyle{\frac{(x_{\nu+1}-\alpha-(x_{\nu}-\alpha))^2}{x_{\nu+2}-\alpha-2(x_{\nu+1}-\alpha)+(x_{\nu}-\alpha)}}\\ &=e_{\nu}-\displaystyle{\frac{(e_{\nu+1}-e_{\nu})^2}{e_{\nu+2}-2e_{\nu+1}+e_{\nu}}}\\ &=\displaystyle{\frac{e_{\nu+2}e_{\nu}-e_{\nu+1}^2}{e_{\nu+2}-2e_{\nu+1}+e_{\nu}}}\\ \displaystyle{\frac{y_{\nu}-\alpha}{e_{\nu+2}}}&=\displaystyle{\frac{e_{\nu}-\frac{e_{\nu+1}^2}{e_{\nu+2}}}{e_{\nu+2}-2e_{\nu+1}+e_{\nu} }}\\ &=\displaystyle{\frac{e_{\nu}-\frac{1}{\kappa+\varepsilon_{\nu+1}}e_{\nu+1}}{e_{\nu+2}-2e_{\nu+1}+e_{\nu} }}\\ &=\displaystyle{\frac{e_{\nu}-\frac{\kappa+\varepsilon_{\nu}}{\kappa+\varepsilon_{\nu+1}}e_{\nu}}{(\kappa+\varepsilon_{\nu+1})(\kappa+\varepsilon_{\nu})e_{\nu}-2(\kappa+\varepsilon_{\nu})e_{\nu}+e_{\nu}}}\\ &=\displaystyle{\frac{1-\frac{\kappa+\varepsilon_{\nu+1}}{\kappa+\varepsilon_{\nu}}}{(\kappa+\varepsilon_{\nu+1})(\kappa+\varepsilon_{\nu})-2(\kappa+\varepsilon_{\nu})+1}}\\ &\rightarrow0(\nu\rightarrow0)\ (\because\ \kappa\neq1) \end{aligned}

が得られる。　\blacksquare)

　他方で\mathrm{Steffensen}は別の反復

\begin{aligned} z_{\nu+1}&=\displaystyle{\frac{z_{\nu}g(g(z_{\nu}))-\left(g(z_{\nu})\right)^2}{g(g(z_{\nu}))-2g(z_{\nu})+z_{\nu}}}\\ &=g(g(z_{\nu}))- \displaystyle{\frac{\left(g(g(z_{\nu}))-g(z_{\nu})\right)^2}{g(g(z_{\nu}))-2g(z_{\nu})+z_{\nu}}} \end{aligned}

を提案し、この方法は\left|g^{\prime}\right|\gt1の場合であっても\alphaに急速に収束する。
　この\{z_{\nu+1}\}は

\begin{aligned} x_{\nu}^{*}=z_{\nu},\ x_{\nu+1}^{*}=g\left(x_{\nu}\right),\ x_{\nu+2}^{*}=g\left(x_{\nu}\right) \end{aligned}

に対して\mathrm{Aitken}加速を施したものに他ならない。しかしこの\mathrm{Steffensen}反復は単なる反復法ではなく、反復列x_{\nu+1}=g(x_{\nu})が発散する場合であっても、初期値z_{0}を\alphaの充分近くに選べば、\{z_{\nu}\}は1次以上の収束をする。すなわち



定理2.8　\mathrm{Steffensen}反復の性質　C^1級関数g(x)がx=g(x)の解\alphaについて、g^{\prime}(\alpha)\neq1ならば、\alphaに充分に近いz_0を初期値とした場合に

\begin{aligned} z_{\nu}-\alpha=o\left(\left|z_{\nu}-\alpha\right|\right) \end{aligned}

であり、特にg(x)がC^2級ならば、

\begin{aligned} z_{\nu}-\alpha=O\left(\left|z_{\nu}-\alpha\right|^2\right) \end{aligned}

が成り立つ。

(\because　関数G(z),g^{*}(z),G^{*}(z)をそれぞれ

\begin{aligned} G(z)&=\begin{cases} \displaystyle{\frac{zg(g(z))-\left(g(z)\right)^2}{g(g(z))-2g(z)+z}},z\neq\alpha\\ \alpha,z=\alpha, \end{cases}\\ g^{*}(z)&=g(z+\alpha)-\alpha,\\ G^{*}(z)&=\begin{cases} \displaystyle{\frac{zg^{*}(g^{*}(z))-\left(g^{*}(z)\right)^2}{g^{*}(g^{*}(z))-2g^{*}(z)+z}},z\neq0\\ 0,z=0 \end{cases} \end{aligned}

とおけば、g^{*}(0)=g(\alpha)-\alpha=\alpha-\alpha=0であり、

\begin{aligned} g^{*}(g^{*}(z))=g(g^{*}(z)+\alpha)-\alpha=g(g(z+\alpha))-\alpha \end{aligned}

である。したがって

\begin{aligned} G(z)-\alpha=G^{*}(z-\alpha) \end{aligned}

が成立する。
　ここでg^{*\prime}=c(=g^{\prime}(\alpha))とおけば、

\begin{aligned} g^{*}(z)&=cz+o(z)\\ g^{*}(g^{*}(z))&=c(cz+o(z))+o(cz+o(z))=c^2z+o(z)\\ g^{*}(g^{*}(z))&-2g^{*}(z)+z=(c-1)^2z+o(z)\\ g^{*}(z)^2&=c^2z^2+o(z^2)\\ zg^{*}(g^{*}(z))&-g^{*}(z)^2=o(z^2) \end{aligned}

が成り立つ。したがってc\neq1ならば、z\rightarrow0のときG^{*}(z)=o(z)が成り立つから、反復の定義式

\begin{aligned} z_{\nu+1}&=\displaystyle{\frac{z_{\nu}g(g(z_{\nu}))-\left(g(z_{\nu})\right)^2}{g(g(z_{\nu}))-2g(z_{\nu})+z_{\nu}}}\\ &=g(g(z_{\nu}))- \displaystyle{\frac{\left(g(g(z_{\nu}))-g(z_{\nu})\right)^2}{g(g(z_{\nu}))-2g(z_{\nu})+z_{\nu}}} \end{aligned}

から、

\begin{aligned} z_{\nu+1}-\alpha=G(z_{\nu})-\alpha=o(z_{\nu}-\alpha) \end{aligned}

が得られる。特にgがC^2級ならば、g^{*}(z)=cz+O(z^2)が成り立つから、同様の議論によりG^{*}(z)=O(z^2),z_{\nu+1}-\alpha=O\left(\left(z_{\nu}-\alpha\right)^2\right)を得る。　\blacksquare)

2.2.4　Pythonによるシミュレーション

　方程式x-e^{-x}=0を数値的に解く。手法として、

  • 二分法
  • 反復法(\varphi(x)=I：恒等写像)
  • \mathrm{Aitken}加速
  • \mathrm{Steffensen}の反復
  • \mathrm{Newton}法

を試してみる。

############################
### 非線形方程式の数値解 ###
############################

# x - e^{-x} = 0を解く

import math

# 初期値関係
a_0 = 0.0 # 2分法の初期値①
b_0 = 2.0 # 2分法の初期値②

x_0 = 1.0 # Newton法等の初期値
eps = 1e-8 # 誤差評価の値

def f(x):
    return x - math.exp(-x)

def f_prime(x):
    return 1 + math.exp(-x)

def g(x):
    return math.exp(-x)

def g_prime(x):
    return -math.exp(-x)

#############
### 2分法 ###
#############
val1 = []

error = float('inf')

a = a_0
b = b_0

c = (a + b) * 0.5
error = abs(a-b) * 0.5

val1.append(c)

while error > eps:
    f_val = f(c)
    
    if f_val > 0:
        b = c
    elif f_val < 0:
        a = c
    else:
        break
    c = (a + b) * 0.5
    val1.append(c)
    error = abs(a-b) * 0.5

##############
### 反復法 ###
##############
val2 = []

x_old = 99999999999.0
x_new = x_0

val2.append(x_new)

while abs(x_new - x_old) > eps:
    x_old = x_new
    x_new = g(x_old)
    
    val2.append(x_new)


######################
### Aitkenの加速法 ###
######################
val3_1 = []
val3_2 = []

x_old = 99999999999.0
x_new = x_0

val3_1.append(x_new)

# 最初にx_1, x_2をつくっておく
for i in range(2):
    x_old = x_new
    x_new = g(x_old)
    val3_1.append(x_new)

y_old = 9999999.0
y_new = -9999999.0

while abs(y_new - y_old) > eps:
    x_old = x_new
    y_old = y_new
    
    y_new = val3_1[len(val3_1) - 3] - (val3_1[len(val3_1) - 2] - val3_1[len(val3_1) - 3]) ** 2 / (val3_1[len(val3_1) - 1] - 2.0 * val3_1[len(val3_1) - 2] + val3_1[len(val3_1) - 3])
    
    x_new = g(x_old)
    
    val3_1.append(x_new)
    val3_2.append(y_new)    

####################
### Steffensen法 ###
####################
val4 = []

x_old = 99999999999.0
x_new = x_0

val4.append(x_new)

while abs(x_new - x_old) > eps:
    x_old = x_new
    x_new = (x_old * g(g(x_old)) - g(x_old) ** 2) / (g(g(x_old)) - 2 * g(x_old) + x_old)
    
    val4.append(x_new)

################
### Newton法 ###
################
val5 = []

x_old = 99999999999.0
x_new = x_0

val5.append(x_new)

while abs(x_new - x_old) > eps:
    x_old = x_new
    x_new = x_old - f(x_old)/f_prime(x_old)
    
    val5.append(x_new)
  • 結果
　　　 二分法
27回
　　　 反復法
33回
　　　 \mathrm{Aitken}加速
19回
　　　 \mathrm{Steffensen}の反復
8回
　　　 \mathrm{Newton}法
4回

今回の復習

今回はなし

参考文献

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