コーシーの積分公式を使って複素周回積分で関数の解を求める-Pythonコードも紹介

コーシーの積分公式を用いることで、関数の解を囲む領域での複素周回積分により、その解の値を求めることができます。今回は、その導出方法とPythonでのプログラムを紹介します。

複素周回積分による解の導出方法

解\(z_1,z_2,\cdots,z_n \)を持つ\(n\)次多項式\(f(z)\)を考える。

$$ f(z) =\prod_{i=1}^n (z-z_i) $$

\(f(z)\)の1次導関数は以下のようになる。

$$ f'(z) = \sum_{i=1}^n \prod_{j=1,j\neq i }^n (z -z_j) $$

これらを用いることで、次式のようになることが分かる。

$$ \frac{f'(z)}{f(z)} = \sum_{i=1}^n \frac{1}{z-z_i} $$

複素平面上の閉曲線\(C\)に沿って周回積分を行う。閉曲線\(C\)内の領域に\(m\)個の\(f(z)\)の解があるとすると、コーシーの積分公式より以下の式が成り立つ。

$$ \begin{eqnarray}
\frac{1}{2\pi i} \oint_C \frac{f'(z)}{f(z)} dz&=& \sum_{i=1}^n \frac{1}{2\pi i} \oint_C \frac{1}{z-z_i} dz\\
&=& \sum_{i=1}^m 1 \\
&=& m
\end{eqnarray} $$

よって上記の積分を実行することで閉曲線\(C\)内にある解の個数が求められる。

続いて次の積分を実行することで、閉曲線\(C\)内にある解のべき乗和を求めることができる。

$$ \begin{eqnarray}
\frac{1}{2\pi i} \oint_C \frac{f'(z)}{f(z)} z^k dz&=& \sum_{i=1}^n \frac{1}{2\pi i} \oint_C \frac{z^k}{z-z_i}dz \\
&=& \sum_{i=1}^m z_i^k \\
&=& z_1^k + z_2^k +\cdots +z_m^k \\
&=&p_k(z_1,z_2,\cdots,z_m) \ \
\end{eqnarray} $$

よって、\(k=1,2,\cdots,m\)におけるべき乗和\( \ p_k(z_1,\cdots,z_m) \ \)を計算することで、ニュートンの恒等式を使って\(m\)次多項式の係数を求めることができる。

関連記事:
多項式の係数と根のべき乗和との関係【ニュートンの恒等式】のPythonプログラム

そしてその多項式をニュートン法やラゲール法などの求根アルゴリズムを使って解くことにより、閉曲線\(C\)内にある解を求めることができる。

関連記事:
ニュートン法による複素求根アルゴリズム~Pythonコード付き~
ラゲール法による求根アルゴリズム～Pythonコード付き～

Pythonプログラム

以下のPythonプログラムでは、まず関数

$$ f(z) = z (z-0.5-0.5i)(z+0.5+0.5i)$$

について、閉曲線\(C : z_1 \rightarrow z_2 \rightarrow z_3 \rightarrow z_ 4 \rightarrow\ z_1 \)で\(\frac{1}{2\pi i}\frac{f'(z)}{f(z)}\)の複素平面上での周回積分積分を行います。

$$ \begin{eqnarray}
z_1 &=& -1-i \\
z_2 &=& 1-i \\
z_3 &=& 1+i \\
z_4 &=& -1+i
\end{eqnarray} $$

そしてその結果より、閉曲線\(C\)内における解の個数\(m\)を求め、その後\(\frac{1}{2\pi i}\frac{f'(z)}{f(z)} z^k \ (k=1,2,\cdots, m)\) についての周回積分を実施しています。これによりべき乗和\(p_k\)を求めることができるので、ニュートンの恒等式を用いて\(m\)次多項式の係数に変換し、それを解くことで閉曲線\(C\)内にある関数\(f(z)\)の解を求めます。

積分に関してはガウス求積により求めています。

関連記事:
ガウス求積で数値積分を行う-Pythonプログラム

# -*- coding: utf-8 -*-

import sympy as sym
import numpy as np


# n次ルジャンドル多項式のゼロ点x_iと重みw_iを求める
def legendre_zeros_weights(n:int):
    x = sym.Symbol('x')
    roots = sym.Poly(sym.legendre(n, x)).all_roots()  # n次ルジャンドル多項式の根を求める
    x_i = np.array([rt.evalf(20) for rt in roots]) 
    w_i = np.array([(2*(1-rt**2)/(n*sym.legendre(n-1, rt))**2).evalf(20) for rt in roots])  # 重み
    #print("x_i=",x_i,"\n","w_i=",w_i)
    return x_i, w_i  # numpy.ndarrayを返す


# ガウス求積により関数funcの[-1,1]における定積分を求める
def gauss_quadrature(func,x,zeros_n=15):  
    # 関数(f(x))、変数(x)、積分点数(zeros_n)

    zeros, weights = legendre_zeros_weights(zeros_n)  
    integ = 0
    for i in range(zeros_n):
        integ += weights[i]*func.subs(x,zeros[i]).evalf()
    return integ  # 積分値


# 複素平面内における周回積分を実行する (経路z1→z2→z3→z4→z1)
def contour_integration(func, x, z1, z2, z3, z4, zeros_n=15):
    # 引数 : 関数f(x), 変数x, 積分経路(矩形の頂点z1~z4), ガウス求積における積分点数 (整数)
    
    integ = 0  # 積分結果
    
    
    # z1 → z2
    mid = (z1 + z2)/2
    delta = (z2 -z1)/2  # 積分範囲を mid + delta*t (t:-1→1)　に変換
    integ += gauss_quadrature(delta/(2*np.pi*1j)*func.subs(x,mid+delta*x), x, zeros_n)
    
    # z2 → z3
    mid = (z2 + z3)/2
    delta = (z3 -z2)/2    
    integ += gauss_quadrature(delta/(2*np.pi*1j)*func.subs(x,mid+delta*x), x, zeros_n)
    
    # z3 → z4
    mid = (z3 + z4)/2
    delta = (z4 -z3)/2    
    integ += gauss_quadrature(delta/(2*np.pi*1j)*func.subs(x,mid+delta*x), x, zeros_n)

    # z4 → z1
    mid = (z4 + z1)/2
    delta = (z1 -z4)/2    
    integ += gauss_quadrature(delta/(2*np.pi*1j)*func.subs(x,mid+delta*x), x, zeros_n)
    
    return complex(integ)


# ニュートンの恒等式 (べき乗和から多項式の係数を求める)
def newton_identities(n,p): 
    # 引数: 次数(n), べき乗和のリスト(p)
    e = [0]*(n+1)  # 対称式
    e[0] = 1
    poly = [1]  # 最高次の係数は1 
    for k in range(1,n+1):
        for i in range(1,k+1):
            e[k] += (-1)**(i-1) * e[k-i] * p[i]  # ニュートンの恒等式
        e[k] /= k
        poly.append((-1)**k * e[k])  
        
    return poly


def root_find(func, z, z1, z2, z3, z4, zeros_n=15):
    # 複素関数func(z)の z1→z2→z3→z4(反時計回り)の中にある零点を求める　

    def g(z, k:int):  # zはSymbolでも数値でも可, k:べき乗の指数
        dfunc = func(z).diff(z)  # 1階微分
        return (dfunc/func(z)*z**k)  # 周回積分するとべき乗和になる (コーシーの積分公式より)
    
    p = []  # べき乗和のリスト (ニュートン多項式の係数)

    # 0次のべき乗和
    integ0 = contour_integration(g(z,0),z, z1, z2, z3, z4, zeros_n)  
    print("解の個数:", integ0)
    n = round(integ0.real)  # 積分経路内の解の個数
    p.append(n)  # 0番目は0乗(=解の個数)

    # 1次以降のべき乗和 
    for i in range(1,n+1):  # 1 ~ n
        p.append(contour_integration(g(z,i),z, z1, z2, z3, z4, zeros_n))  # i次のべき乗和をリストpに加える

    poly = newton_identities(n,p)  # 多項式の係数のリスト
    f_poly = np.poly1d(poly)  # 多項式
    return f_poly.r  # 根を得る


# メイン実行部
def main():
    
    z = sym.Symbol('z')
    
    # 解を求めたい関数
    def f(z): 
        return z*(z-(0.5+0.5j))*(z-(-0.5-0.5j))  # 積分したい関数
    
    # 周回積分の範囲 z1→z2→z3→z4　(複素平面で反時計回り)
    z1 = -1-1j
    z2 = 1-1j
    z3 = 1+1j
    z4 = -1+1j
    
    roots = root_find(f,z, z1, z2, z3, z4, zeros_n=15)    
    print("解: ",['{:.10f}'.format(rt) for rt in roots])
    

if __name__ == '__main__':
    main()

関数の解と積分経路は以下のグラフのようになります。

実行結果

解の個数: (3.0000000314797544-2.117582368135751e-22j)
解:  ['0.5000000039+0.5000000039j', '-0.5000000039-0.5000000039j', '0.0000000000+0.0000000000j']

最初に設定した解、\(0, 0.5+0.5i, 0.5-0.5i\)に近い値が求められていることが分かります。