ミューラー法で方程式の数値解を求める-Pythonコード

ミューラー(Muller)法は、方程式の解を数値的に求める手法の一つです。ニュートン法などと違い、関数の導関数を与える必要がないこと、初期値を複素数にしなくても複素数解を得られることなどの長所があります。

1. ミューラー(Muller)法の導出
2. ミューラー(Muller)法のPythonプログラム

ミューラー(Muller)法の導出

Muller法は、

$$f(z) = 0$$

となる方程式の解$z$を求める手法の一つである。

まず最初に解の近似値として三つの値$z_1,z_2,z_3$を与える。次に、$ (z_1,f(z_1)),(z_2,f(z_2),(z_3,f(z_3)$の3点を通る二次関数、

$$ y = a z^2 + b z + c $$

を求める。すなわち以下の連立方程式を解いて、$a,b,c$を求める。

$$ \begin{eqnarray}
f(z_1) &=& a z_1^2 + b z_1 + c \\
f(z_2) &=& a z_2^2 + b z_2 + c \\
f(z_3) &=& a z_3^2 + b z_3 + c \\
\end{eqnarray} $$

これは行列を用いて下記のように求めることができる。

$$
\begin{pmatrix}
a \\
b \\
c \\
\end{pmatrix}
=
\begin{pmatrix}
z_1^2 & z_1 & 1 \\
z_2^2 & z_2 & 1 \\
z_3^2 & z_3 & 1 \\
\end{pmatrix}^{-1}
\begin{pmatrix}
f(z_1) \\
f(z_2) \\
f(z_3) \\
\end{pmatrix}
$$

これを解いて$a,b,c $を求め、二次方程式

$$ a z^2 + b z + c = 0$$

の解を求め、$z_3$に近いほうを解の近似値として$z_4$とする。再度$z_2,z_3,z_4$を用いて同様のことを行い$z_5$を求める。以降、これを繰り返し、$z_{k-2},z_{k-1},z_{k}$を用いて$z_{k+1}$を順次求めていく。

Muller法の場合には上記の二次方程式を解く時点で複素数が現れるので、初期値に複素数を設定しなくても複素数解が求められる。

ミューラー(Muller)法のPythonプログラム

以下のPythonプログラムでは$z^3+1=0$の解をMuller法で求めています。

# -*- coding: utf-8 -*-
"""
ミュラー(Muller)法で方程式の解を求める
"""

import sympy as sym
import numpy as np
import matplotlib.pyplot as plt

z = sym.Symbol('z')
f = z**3 + 1  # 根を求めたい関数

EPS = 1E-12  # 収束条件
LIMIT = 50  # 最大繰り返し回数


# ミュラー法
def muller_method(f,z0,z1,z2):  # f:関数, z0:初期値
    counter = 0
    while True:
        f0 = f.subs(z,z0)
        f1 = f.subs(z,z1)        
        f2 = f.subs(z,z2)
        
        A = np.array( [[z0**2,z0,1],
                       [z1**2,z1,1],
                       [z2**2,z2,1]])
        B = np.array([[f0],
                      [f1],
                      [f2]])
        inv_A = np.linalg.inv(A)  # 逆行列
        coef = np.dot(inv_A,B)  # [a,b,c] :(az^2 + bz + c)
        inroot = complex(coef[1]**2 - 4*coef[0]*coef[2])  # complex型に変換
        z3p = (-coef[1] + np.sqrt(inroot)) / (2*coef[0])  # プラスの解
        z3m = (-coef[1] - np.sqrt(inroot)) / (2*coef[0])  # マイナスの解
        z3p = complex(z3p)
        z3m = complex(z3m)
        z3 = z3p if abs(z2 - z3p)<abs(z2 - z3m) else z3m  # z2に近い方をz3にする
        
        if abs(complex(f.subs(z,z3)))<EPS:  # 収束条件
            break
        if counter == LIMIT:  # 計算回数の上限
            break
        
        z0, z1, z2 = z1, z2, z3  # (z1,z2,z3) -> (z0,z1,z2)に入れ替える 
        
        print("|root|=",abs(z3))        
        counter += 1        
    return z3


zroot = muller_method(f,-1,1,2)
print(zroot)

結果は以下のようになりました。

(0.4999999999999688-0.8660254037844678j)

解析解の一つ、

$$ z = \frac{1-\sqrt{3}i}{2}$$

に近い値が得られています。

ただしMuller法の場合、判定精度を小さくしすぎると解が発散(不安定)してしまうようです。上記のPythonコードで

EPS = 1E-16  # 収束条件

とした場合の、$f(z_k)$の挙動が以下のグラフのようになっており、計算が繰り返されるごとに発散気味になっています。

この対策法として個人的には、$z_{k-2},z_{k-1},z_{k}$と$ z_{k+1}$から新たな3つを選ぶときに、近い値同士で選んであげるようなアルゴリズムにすればいいのかなと思います。

前の記事

ニュートン法で関数のすべての複素根を求める-Pythonプログラム 2021.03.21
次の記事

ガウス求積で数値積分を行う-Pythonプログラム 2021.03.23

ミューラー法で方程式の数値解を求める-Pythonコード

ミューラー(Muller)法の導出

ミューラー(Muller)法のPythonプログラム

ニュートン法で関数のすべての複素根を求める-Pythonプログラム 2021.03.21

ガウス求積で数値積分を行う-Pythonプログラム 2021.03.23