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

方程式の解を求める方法の一つ、ニュートン法(Newton’s method)またはニュートン・ラフソン法(Newton-Raphson method)を使って、複数の解(関数のゼロ点)を求める方法とPythonプログラムを紹介します。

1. ニュートン法で多数の解を求める方法
2. 複数の複素根を求めるPythonコード
3. まとめ

ニュートン法で多数の解を求める方法

ニュートン法とは、

$$f(z) = 0$$

となる方程式の解$z$を求める手法で、具体的には$ z_0 $を解の十分近くの値として適当に設定し、

$$ z_{n+1} = z_n – \frac{f(z_n)}{f'(z_n)} \tag{1}$$

のように、$ z_0, z_1, z_2\cdots $を逐次的に求めていき、解の近似値を得るアルゴリズムです。

関数$f(z)$が複数の根(零点)を持っている場合には、以下のような方法でそれらを求めることができます。

まず、方程式$f(z)=0$にニュートン法を適用し、一つの解$\alpha_1$が求められたとします。そのとき、$f(z)$から次の関数$f_1(z)$を作ります。

$$ f_1(z) = \frac{f(z)}{z-\alpha_1} $$

この$f_1(z)$を使い、次の方程式を考えます。

$$ f_1(z) = 0 $$

この方程式について再度ニュートン法を適用することで、先ほどとは別の解$\alpha_2$を得ることができます。さらに他の解が存在する場合には、

$$ f_2(z) = \frac{f_1(z)}{z-\alpha_2} $$

として、方程式$f_2(z)=0$にニュートン法を適用します。以下、これを繰り返していくことで関数$f(z)$の複数の根を求めることができます。

代数方程式の場合には、複素数の範囲での解の個数は多項式の次数$n$に一致するので、上記の処理を$n$回繰り返せばすべての解を求めることができます。一方、超越方程式( $ \sin z, \mathrm{e}^{z} $などを含むもの)については解の個数は一般には分からないので、あらかじめグラフをプロットするなどして解の個数を調べておく必要があります。

複数の複素根を求めるPythonコード

以下のプログラムでは、

$$f(z) = z^3 – 1 $$

に対して複素数の範囲でニュートン法を適用し、すべての根を求めています。

# -*- coding: utf-8 -*-
"""
Newton法で方程式の複数の根を求める
"""

import sympy as sym
import cmath

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


# ニュートン法
def newton_method(f,z0):  # f:関数, z0:初期値
    df = sym.diff(f)  # 関数の1階導関数

    counter = 0
    while True:
        if df.subs(z,z0)==0:  # 導関数が0のときはニュートン法を適用できない
            print("傾き0のため計算停止")
            break

        z_new = complex(z0 - f.subs(z,z0)/df.subs(z,z0)) # z1 = z0 - f(z0)/f'(z0)
        if abs(complex(f.subs(z,z_new))) < 1E-15:  # 収束条件
            break
        if counter == 100:  # 計算回数の上限
            break

        z0 = z_new
        counter += 1
    return z_new


#  以下で、ニュートン法を複数回繰り返す。
f = f0
for num in range(100):  # 適当な回数繰り返す
    zroot = newton_method(f,1+1j)  # 実数解のみ求めたい場合は、初期値を実数にする
    if cmath.isinf(zroot) or cmath.isnan(zroot) :  # 解が出尽くしたら無限かnanが戻り値になる
        break
    print(zroot, ", |f0(z)|=",abs(complex(f0.subs(z,zroot))))
    f = f / (z-zroot)

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

(1+5.058454998933075e-31j) , |f0(z)|= 1.5175364996799225e-30
(-0.5+0.8660254037844386j) , |f0(z)|= 1.5052626332710354e-16
(-0.49999999999999983-0.8660254037844385j) , |f0(z)|= 6.953149470236021e-16

一方、$f(z) = z^3 – 1 $の解析解は、

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

です。

Pythonでもsympy.solve()を使って下記のように求められます。

# -*- coding: utf-8 -*-
import sympy as sym
import cmath

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

roots = sym.solve(f0) 
roots_deci = map(complex,roots)  # 小数点表示にする

print(list(roots_deci))

計算結果

[(1+0j), (-0.5-0.8660254037844386j), (-0.5+0.8660254037844386j)]

ニュートン法で求めた数値解と解析解はほぼ一致していることが分かります。

まとめ

ニュートン法を使って、方程式の複数の複素数解を求める方法を紹介しました。

求根アルゴリズムについては以下の関連記事も書いているので、是非ご覧ください。

前の記事

Jenkins-Traubアルゴリズムで多項式の根を求める～Pythonコード付き～　　 2021.03.20
次の記事

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

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

ニュートン法で多数の解を求める方法

複数の複素根を求めるPythonコード

まとめ

Jenkins-Traubアルゴリズムで多項式の根を求める～Pythonコード付き～ 2021.03.20

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

Jenkins-Traubアルゴリズムで多項式の根を求める～Pythonコード付き～　　 2021.03.20