【Python】複素線積分の数値的解法 【数値計算】

はじめに

脳磁逆問題を解く際に複素線積分する必要があっただが、Pythonで複素線積分をしているものが見つからなかったので、記事を書こうと思いました。
まあ、複素積分って実質は線積分に過ぎないので、線積分が数値的に求められたら、それを複素数にして終わりです。

構成としては、

1. 数値積分の基礎理論を触れる
2. 閉曲線に沿った数値積分の理論について触れる
3. 実装

って感じでいきます。

理論

数値積分の基礎

数値積分には、代表的に台形則や、中点則、シンプトン則があるが、それぞれ、積分を「求積問題」として、求める面積を細かく区間で分けて、どのように近似していくかという概念からきてるものであった。 それぞれの近似の概念としては以下のような感じ

• 台形則...積分を狭い台形の足し合わせとして近似
• 中点則...積分を微小区間の中点を高さの代表値とする幅の狭い長方形の足し合わせとして近似
• シンプソン則...微小区間において関数$f(x)$を2次曲線で近似し、微小区間の端点とその中点とでの近似2次曲線上の点を用いて台形を作って積分を近似。(台形則の精度をあげている)

基本的に、下にいくほうが精度が高い(はず)

以下の図はシンプソン則の概念をわかりやすく図に起こしたもの。(重田 出 : 数値積分と数値微分（基礎）から引用)

図:シンプソ則について

閉曲線に沿った数値積分の理論

閉じた曲線$C$に沿う複素積分 $$ I := \oint_C f(z) \mathrm{d}z \tag{1} $$ について考えた時、まず、閉じた曲線を円に限る。(任意の閉曲線での線積分って円周での線積分で代用できなかったっけ?)
この時、文献1によると円に沿った線積分は上端も下端もなく、シンプソン則をはじめとした積分の上端と下端からの効果を補償したものだから、そんなことを考えない台形則で十分らしい。
ってな訳で、台形則だけで数値積分を行う。

という訳で、式(1)は

$$ I \simeq \sum_{k=0}^{N -1 } ( f( z_k ) + f(z_{ k -1 }) ) \frac{z_{k +1 } - z_k}{2} \tag{2} $$

で数値積分を行うが、この時、$z = z_0 + r\mathrm{e}^{i\theta}$は円周上と限られているので、極座標での数値積分を行うといい。
という訳で、結局考える積分は

$$ I = \int_0^{2\pi} f(z_0 + r\mathrm{e}^{i\theta}) ri\mathrm{e}^{i\theta} \mathrm{d}\theta \tag{3} $$

となる。これに台形則を適用すると、

$$ I \simeq \sum_{k=0}^{N -1 } f(z_0 + r\mathrm{e}^{i\theta_k}) ri\mathrm{e}^{i\theta_k} \Delta \theta \tag{4} $$

となる。ただし、

$$ \Delta \theta = \frac{2\pi}{N}, \hspace{0.5cm} \theta_k = k\Delta \theta \tag{5} $$

とした。

理論は以上。

数値複素積分を実際にやってみる

さて、理論も済んだし、複素積分をやっていきましょ。題名の通りPythonを使います。理由は僕が一番使い慣れてるから。
という訳で実際のコードは以下のよう。

import numpy as np

# 複素積分
_int_step_num = 10**3
_int_center_point = 0 + 0j
_int_r = 1.1

def numerical_int( func, z0=_int_center_point, r=_int_r, Ni=_int_step_num ):
    """"関数fを複素積分する。
    z０は積分半径の中心。rは積分経路の半径。Niは刻み数。
    数値積分の値を返す。"""
    
    I = 0 # 積分値
    delta_theta = 2*np.pi/Ni
    
    for k in range(Ni):
        theta_k = k*delta_theta
        if np.isnan(func(z0 + r* np.exp(complex(0, theta_k )))):
            continue
        I += func(z0 + r* np.exp(complex(0, theta_k )))*   r*1j * np.exp(complex(0, theta_k )) *delta_theta
        
    return I

刻み数は実験的にやったら、刻み数が$10^{3}$ぐらいでも精度良さげという適当な理由で決めた。

実験

$1/z$で実験してみた。厳密解は$2\pi i$。

図:実験結果

良さげ。

以上！！！

参考文献

• 1: 殿塚 勲 : 複素関数の積分に対する数値積分の試み
• 2: 重田 出 : 数値積分と数値微分（基礎）

最後に

これはアドベントカレンダーのProny法の実装に繋がる準備としてブログ化しました。アンド、将来僕がまた複素積分するときに困らないように。 今月は、これと、

• N次連立方程式を一気に求める方法
• Prony法の実装
• 多浪体験記①

の更新が目標。(できるとは言ってない。) 今月は5回の更新を目指すぞー！おー！

では、今日はここまで！！バイバイ！！！