更新日: 2026年2月8日

X. 複素積分の定理公式を数値積分で確認（２）

5. Cauchyの積分定理の確認

Cauchyの積分定理というのは以下の通りです。

単連結領域\(D\)において \(f(z)\) が正則であれば、\(D\)内のすべての単純閉曲線\(C\)に対して
\[ \oint_C f(z) \, dz = 0 \tag{7} \label{g} \]
となる。
積分記号に○が付いているのは閉曲線に沿った周回積分だからです。単連結領域\(D\)というのは、\(D\)内のすべての閉曲線が\(D\)内の点だけを囲むような領域、つまり穴が開いていない領域のことです。単純閉曲線というのは自分自身と交わらないループ状の曲線です。4節で確認した正則関数の積分値が経路によらないという結果は数学的にはCauchyの積分定理から証明されますが、Cauchyの積分定理についても複素数値積分で確認してみます。

下図はZIPでダウンロードできるcurve2num.xlsmのSheet2になりますが、今回用いる積分経路はこの閉曲線\(C\)です。「4.1 フリーハンド曲線の頂点の座標取得」で紹介したVBAマクロ「Sub 第1曲線頂点()」をそのまま使って頂点の座標を取得しています。ここでは、被積分関数

\[ f(z) = \frac{1}{z} \tag{8} \label{e} \]
に対してCauchyの積分定理を確認しますが、閉曲線\(C\)の内側に原点 \(0 + 0 \, i\) を含まないというのがポイントになります。この閉曲線\(C\)は「6. Cauchyの積分公式の確認」でも、「7. 留数定理の確認」でも使いますが、前者で使う理由は \(-1 + 1 \, i\) を含むから、後者は \(-2 + 0 \, i\) と \(0 + 2 \, i\) を含むからで、三者三様に異なります。

画像を新しいタブで開くと高解像度で見ることができます。

式\eqref{e}の \(f(z)\) は原点においてのみ微分できません。従って、内側に原点を含まない閉曲線\(C\)に沿って積分すれば、結果は0になるはずです。この複素数値積分をするPythonプログラムは以下の通りです。基本的には「4.2 取得した座標に沿った数値積分」にあるcplx_intgrl_2.pyと同じ構成ですが、積分経路が一つだけなので変数mに対するループはありません。それ以外の変更点は、4行目でsheet2を参照しており、21行目の被積分関数が式\eqref{e}になっていることです。26行目にある計算結果を見ると有効数字8～9桁で0 + 0jとなっており、Cauchyの積分定理通りになっています。

コードの表示にはPrism.jsを使用しています。
prism.cssの設定は団塊爺ちゃんの備忘録を参考にしました。 

import openpyxl

wb = openpyxl.load_workbook('(path to the xlsm file)/curve2num.xlsm')
sheet = wb['Sheet2']
i0 = 12
div = 1000

vertNum = sheet.cell(row=1, column=i0+2).value
vertArX = []
vertArY = []
for n in range(vertNum):
    vertArX.append(sheet.cell(row=n+3, column=i0+1).value)
    vertArY.append(sheet.cell(row=n+3, column=i0+2).value)

Sum = 0. + 0.j
for k in range(vertNum - 1):
    z1 = complex(vertArX[k], vertArY[k])
    z2 = complex(vertArX[k+1], vertArY[k+1])
    for p in range(div):
        z3 = z1 + (z2 - z1) * (0.5 + p) / div
        Sum += (1 / z3) * (z2 - z1) / div
print(f"{Sum.real:.10f} + {Sum.imag:.10f}j")


# ----- Output -----
# -0.0000000096 + -0.0000000359j

6. Cauchyの積分公式の確認

Cauchyの積分公式というのは以下の通りです。

単連結領域\(D\)において正則な関数 \(f(z)\) について
\[ f(z_{0}) = \frac{1}{2\pi i} \, \oint_C \frac{f(z)}{z - z_{0}} \, dz \tag{9} \label{f} \]
となる。ただし、\(z_{0}\)は\(D\)内の任意の点、\(C\)は\(z_{0}\)を取囲む\(D\)内の任意の単純閉曲線である。
式\eqref{g}と式\eqref{f}の積分は何か紛らわしいですよね。両式とも \(f(z)\) は単連結領域\(D\)において正則ですが、式\eqref{f}では被積分関数が \(f(z)\) を \((z - z_{0})\) で割ったものになっているのがミソです。\((z - z_{0})\) で割った関数全体は、\(z = z_{0}\) において微分不可能になります。そういう微分不可能な点が内側にある閉路で周回積分すると、その積分値は微分不可能な点\(z_{0}\)における \(f(z)\) の値 \(f(z_{0})\) の \(2\pi i\) 倍になるというのです。この摩訶不思議な公式の数学的な証明は別途教科書で学んでいただくとして、ここでは例によって複素数値積分による確認を行います。

積分経路には5節の閉曲線\(C\)を用いるので、その内側にある点を選んで \(z_{0} = -1 + 1 \, i\) とします。\(f(z)\) には3節と4節で用いた式\eqref{h}を用いますが、ここに再掲します。

\[ f(z) = z^{2} + 3z + 1 \tag{3} \label{h} \]
このお膳立てで式\eqref{f}右辺にある複素数値積分を行いますが、左辺に \(2\pi i\) を移項した \(2\pi i\,f(z_{0})\) についても別途数値計算して積分値と比較します。以下がそのプログラムです。24行目までの積分プログラムは今までに出てきたものと同じですが、23行目の被積分関数は式\eqref{h}を \((z - z_{0})\) で割ったものになっています。一方、27行目は \(2\pi i\,f(z_{0})\) を数式通りに計算しているだけです。両者の計算結果を比較すると、有効数字8～9桁で一致しているのが分ります。 


import openpyxl
import math

wb = openpyxl.load_workbook('(path to the xlsm file)/curve2num.xlsm')
sheet = wb['Sheet2']
i0 = 12
div = 1000

vertNum = sheet.cell(row=1, column=i0+2).value
vertArX = []
vertArY = []
for n in range(vertNum):
    vertArX.append(sheet.cell(row=n+3, column=i0+1).value)
    vertArY.append(sheet.cell(row=n+3, column=i0+2).value)

Sum = 0. + 0.j
z0 = complex(-1., 1.)
for k in range(vertNum - 1):
    z1 = complex(vertArX[k], vertArY[k])
    z2 = complex(vertArX[k+1], vertArY[k+1])
    for p in range(div):
        z3 = z1 + (z2 - z1) * (0.5 + p) / div
        Sum += (pow(z3, 2) + 3 * z3 + 1) / (z3 - z0) * (z2 - z1) / div
print(f"Integral:  {Sum.real:.10f} + {Sum.imag:.10f}j")

# Calculation for 2πj∙f(z0)
dpj_fz0 = complex(0., 2. * math.pi) * (pow(z0, 2) + 3 * z0 + 1)
print(f"2πj∙f(z0): {dpj_fz0.real:.10f} + {dpj_fz0.imag:.10f}j")


# ----- Output -----
# Integral:  -6.2831852832 + -12.5663706052j
# 2πj∙f(z0): -6.2831853072 + -12.5663706144j

Cauchyの積分公式から導かれるGoursatの定理についても数値積分で確認します。Goursatの定理は以下の通りです。\(n = 0\) の時がCauchyの積分公式です（拡張された階乗の定義で \(0 \, ! = 1\)）。

単連結領域\(D\)において正則な関数 \(f(z)\) について、\(D\)内の点 \(z = z_{0}\) における\(n\)階微分の値 \(f^{(n)}(z_{0})\) は
\[ f^{(n)}(z_{0}) = \frac{n \, !}{2\pi i} \, \oint_C \frac{f(z)}{(z - z_{0})^{n+1}} \, dz \qquad (n = 0, 1, 2, \dots) \tag{10} \label{o} \]
となる。

確認は1階微分などとケチ臭いことを言わず、2階微分にしてみます。

\[ f^{(2)}(z_{0}) = \frac{1}{\pi i} \, \oint_C \frac{f(z)}{(z - z_{0})^{3}} \, dz \tag{11} \label{i} \]
となるわけですが、積分経路には5節の閉曲線\(C\)を用い、 \(z_{0} = -1 + 1 \, i\) とします。\(f(z)\) には式\eqref{h}を用いることにすると、2階導関数は簡単に求められて
\[ f^{(2)}(z) = 2 \tag{12} \]
です。従って、式\eqref{i}右辺の積分値は \(2\pi i\) となるはずです。以下のプログラムで数値積分すると有効数字7桁で一致していますが、精度は若干悪くなるようです。 


import openpyxl
import math

wb = openpyxl.load_workbook('(path to the xlsm file)/curve2num.xlsm')
sheet = wb['Sheet2']
i0 = 12
div = 1000

vertNum = sheet.cell(row=1, column=i0+2).value
vertArX = []
vertArY = []
for n in range(vertNum):
    vertArX.append(sheet.cell(row=n+3, column=i0+1).value)
    vertArY.append(sheet.cell(row=n+3, column=i0+2).value)

Sum = 0. + 0.j
z0 = complex(-1., 1.)
for k in range(vertNum - 1):
    z1 = complex(vertArX[k], vertArY[k])
    z2 = complex(vertArX[k+1], vertArY[k+1])
    for p in range(div):
        z3 = z1 + (z2 - z1) * (0.5 + p) / div
        Sum += (pow(z3, 2) + 3 * z3 + 1) / pow((z3 - z0), 3) \
               * (z2 - z1) / div
print(f"Integral: {Sum.real:.10f} + {Sum.imag:.10f}j")

# Calculation for 2πj
dpj = complex(0., 2. * math.pi)
print(f"2πj     : {dpj.real:.10f} + {dpj.imag:.10f}j")


# ----- Output -----
# Integral: 0.0000001012 + 6.2831854978j
# 2πj     : 0.0000000000 + 6.2831853072j

7. 留数定理の確認

留数定理というのは以下の通りです。

単連結領域\(D\)に単純閉曲線\(C\)があり、\(f(z)\) は\(C\)で囲まれた領域で有限個の点 \(c_{1}, c_{2}, \dots c_{n}\) を除いて正則であるとき
\[ \frac{1}{2\pi i} \, \oint_C f(z) \, dz = \sum_{k=1}^{n} \, \mathrm{Res} \, (f(z), \, c_{k}) \tag{13} \label{j} \]
である。ただし、\(\mathrm{Res}(f(z), \, c_{k})\) は \(f(z)\) の \(z = c_{k}\) での留数とする。
では、留数とは何なのでしょう。上記の有限個の点 \(c_{1}, c_{2}, \dots c_{n}\) のことを特異点と呼びますが、特異点が1個しかない場合の式\eqref{j}がそのまま留数の定義になります。\(f(z)\) が\(C\)で囲まれた領域で1点 \(c\) を除いて正則であるとき、\(f(z)\) の \(z = c\) での留数 \(\mathrm{Res} \, (f(z), \, c)\) は
\[ \mathrm{Res} \, (f(z), \, c) \equiv \frac{1}{2\pi i} \, \oint_C f(z) \, dz\ \tag{14} \label{k} \]
と定義されます。式\eqref{j}と式\eqref{k}の言っていることを合体すると、単純閉曲線\(C\)の内側に複数の特異点がある場合、\(C\)に沿った周回積分は個々の特異点を囲む単純閉曲線 \(C_{1}, C_{2}, \dots\) に沿った周回積分の和になるということです。

7.1 留数自体も周回積分で計算

まずは留数を式\eqref{k}の周回積分で計算して、全体の周回積分が個々の周回積分の和になるということを確かめてみます。被積分関数は複数の特異点を持つ必要があるので

\begin{align} f(z) &= \frac{z \, e^{\pi z}}{z^{4} - 16} \\ &= \frac{z \, e^{\pi z}}{(z - 2) \, (z + 2) \, (z - 2i) \, (z + 2i)} \tag{15} \label{l} \end{align}
とします。分母を展開した式から明らかなように、\(f(z)\) の特異点は \(\pm \, 2\), \(\pm \, 2 \, i\) の4つです。これをcurve2num.xlsmのSheet2にある閉曲線\(C\)に沿って積分すると、\(C\)の内側にある特異点は\(-2\)と\(2 \, i\)の2つになります。この2つの特異点の留数を求めるためには各特異点のみを囲む閉曲線\(C_{1}\), \(C_{2}\)が必要になりますが、\(C_{1}\), \(C_{2}\)はcurve2num.xlsmのSheet3（下図）で座標を取得します。

画像を新しいタブで開くと高解像度で見ることができます。

実際に周回積分を行うプログラムは以下のようになっています。10行目でSheet2を参照し、11～25行目で閉曲線\(C\)に沿った式\eqref{l}の積分を計算しています。28行目ではSheet3を参照し、30～47行目で閉曲線\(C_{1}\), \(C_{2}\)に沿った積分を計算します。後者では閉曲線が2つあるので「4.2 取得した座標に沿った数値積分」のcplx_intgrl_2.pyの時のように変数mによるforループを使います。後者の2つの積分値は48行目でSumArというリストに格納し、49～50行目でその和を計算して表示します。53行目のC:と56行目のC1+C2:を比較すると、小数点以下7桁まで一致しています。 


import openpyxl
import math
import cmath

wb = openpyxl.load_workbook('(path to the xlsm file)/curve2num.xlsm')
i0 = 12
div = 1000

# Cに沿った周回積分を計算
sheet = wb['Sheet2']
vertNum = sheet.cell(row=1, column=i0+2).value
vertArX = []
vertArY = []
for n in range(vertNum):
    vertArX.append(sheet.cell(row=n+3, column=i0+1).value)
    vertArY.append(sheet.cell(row=n+3, column=i0+2).value)
Sum = 0. + 0.j
for k in range(vertNum - 1):
    z1 = complex(vertArX[k], vertArY[k])
    z2 = complex(vertArX[k+1], vertArY[k+1])
    for p in range(div):
        z3 = z1 + (z2 - z1) * (0.5 + p) / div
        Sum += z3 * cmath.exp(math.pi * z3) / (pow(z3, 4) - 16.) \
               * (z2 - z1) / div
print(f"C:    {Sum.real:.10f} + {Sum.imag:.10f}j")

# C1, C2に沿った周回積分を計算
sheet = wb['Sheet3']
SumAr = []
for m in range(2):
    col_m1 = i0 + m * 2 + 1
    col_m2 = i0 + m * 2 + 2
    vertNum = sheet.cell(row=1, column=col_m2).value
    vertArX = []
    vertArY = []
    for n in range(vertNum):
        vertArX.append(sheet.cell(row=n+3, column=col_m1).value)
        vertArY.append(sheet.cell(row=n+3, column=col_m2).value)
    Sum = 0. + 0.j
    for k in range(vertNum - 1):
        z1 = complex(vertArX[k], vertArY[k])
        z2 = complex(vertArX[k+1], vertArY[k+1])
        for p in range(div):
            z3 = z1 + (z2 - z1) * (0.5 + p) / div
            Sum += z3 * cmath.exp(math.pi * z3) / (pow(z3, 4) - 16.) \
                   * (z2 - z1) / div
    print(f"C{m + 1}:    {Sum.real:.10f} + {Sum.imag:.10f}j")
    SumAr.append(Sum)
Sum = sum(SumAr)
print(f"C1+C2: {Sum.real:.10f} + {Sum.imag:.10f}j")

# ----- Output -----
# C:    -0.0000000245 + -0.3919657209j
# C1:    0.0000000000 + 0.0007333430j
# C2:    0.0000000009 + -0.3926990832j
# C1+C2: 0.0000000009 + -0.3919657401j

7.2 Laurent展開に基づく検算

ここまでの議論では留数をわざわざ定義する意味が分らないのですが、留数は周回積分以外にも様々な解析的手法で計算することができます。その基本になるのが、留数は \(f(z)\) をLaurent展開した時の \((z - c)^{-1}\) の項の係数 \(a_{-1}\) であるということです。正則関数は実関数と同様にTaylor展開できますが、1点\(c\)を除いて正則な関数はその特異点\(c\)のまわりでLaurent展開できます。

\(0 < | \, z - c \, | < R\) で正則な複素関数 \(f(z)\) は以下のようにべき級数展開できる。
\[ f(z) = \sum_{n=-\infty}^{\infty} \, a_{n} \, (z - c)^{n} \tag{16} \label{n} \]
ただし，各係数 \(a_{n}\) は
\[ a_{n} = \frac{1}{2\pi i} \, \oint_{|z-c|=r} \frac{f(z)}{(z - c)^{n+1}} \, dz \tag{17} \label{m} \]
で計算できる（\(r\)は \(0 < r < R\) を満たす実数なら何でもよい）。
式\eqref{m}で \(n=-1\) とした時の右辺は式\eqref{k}の右辺と同じになるので、確かに
\[ \mathrm{Res} \, (f(z), \, c) = a_{-1} \tag{18} \label{p} \]
になっています。逆に式\eqref{p}を留数の定義としている教科書もあります。

Laurent展開とTaylor展開の違いは、式\eqref{n}における\(\Sigma\)の下限が0ではなく\(-\infty\)になっており、負べきの項があるという点です。また、Taylor展開の\(n\)乗の項の係数は

\[ a_{n} = \frac{f^{(n)}(c)}{n \, !} \tag{19} \label{q} \]
なので、式\eqref{m}とは異なります。しかし、式\eqref{o}（Goursatの定理）を思い返すと、もし \(f(z)\) が点\(c\)で微分可能ならLaurent展開の \(a_{0}, a_{1}, a_{2}, \dots\) は式\eqref{q}になると得心がいきます。つまりLaurent展開は点\(c\)が特異点である場合にも適用できるようにTaylor展開を拡張したものであり、点\(c\)が特異点でない場合にはそのままTaylor展開になるということです。

特異点は、Laurent展開した時に負べきの項がどこまで続くかによって分類できます。負べきの項が \((z - c)^{-m}\) までしかない（\(a_{-m}\) より先の係数が全て0となる）場合、その特異点は\(m\)位の極と言います。負べきの項が無限に続く場合は真性特異点と呼びます。例えば、式\eqref{l}の4つの特異点は全て1位の極です。特異点 \(-2\) について言うと、式\eqref{l}の2番目の等式から \((z + 2)\) を取除くと残った関数は \(z = -2\) において微分可能です。従って、残った関数はTaylor展開可能であり、元の関数 \(f(z)\)はこれを \((z + 2)\) で割ったものになるので、負べきの項は1つだけになります。

さて、上記の \((z + 2)\) を取除いた式において \(z = -2\) とするとどうなるでしょう。Taylor展開で残る項は \((z + 2)^{0}\) の係数である定数項のみとなり、これがLaurent展開の \(a_{-1}\) となります。従って、 \((z + 2)\) を取除いた式に \(z = -2\) を代入した結果が \(z = -2\) での留数になります。

\begin{align} \mathrm{Res} \, (f(z), \, -2) &= \frac{z \, e^{\pi z}}{(z - 2) \, (z^{2} + 4)} \; \Biggr|_{z=-2} \\[5pt] &= \frac{e^{-2\pi}}{16} \tag{20} \end{align}
周回積分と比較しやすくするために、これに \(2\pi i\) をかけると
\begin{align} 2\pi i \, \mathrm{Res} \, (f(z), \, -2) &= \frac{\pi \, e^{-2\pi}}{8} \, i \\[5pt] &= 0.0007333430 \, i \tag{21} \label{r} \end{align}
となります。式\eqref{r}の右辺は純虚数で、虚数単位を除いた式を計算をするには関数電卓（アプリ）で十分ですが、一応Pythonで計算しています。

同様に \(z = 2i\) での留数も計算すると

\begin{align} \mathrm{Res} \, (f(z), \, 2i) &= \frac{z \, e^{\pi z}}{(z^{2} - 4) \, (z + 2i)} \; \Biggr|_{z=2i} \\[5pt] &= - \, \frac{1}{16} \tag{22} \end{align}
\begin{align} 2\pi i \, \mathrm{Res} \, (f(z), \, 2i) &= - \, \frac{\pi}{8} \, i \\[5pt] &= -0.3926990817 \, i \tag{23} \label{s} \end{align}
となります。式\eqref{r}の数値を前項のPythonプログラムresidue_theorem.pyで計算した周回積分の値（54行目のC1:）と比較すると、小数点以下10桁まで完全に一致しています。式\eqref{s}の数値を55行目のC2:と比較すると、小数点以下8桁まで一致しています。定義通りの周回積分で数値計算した留数と、Laurent展開の \(a_{-1}\) として解析的に計算した留数が一致するというのは数学的には当り前の話ですが、プログラマーとしては嬉しい結果です。

「理屈はそうだけど、その通りにはいかないよね」というのが現実世界だという話をもう少し掘り下げます。これは何も技術者が「理屈」を軽視しているということではありません。現実世界が理屈通りにいかないのは、その理屈では考慮されていない要因が色々とあるからです。数学の世界は土俵となる世界そのものを数学者が構築しています。ですから「そこまでは考えていなかったなあ」という不意打ちを食らう心配がないのです。いくら用意周到に考えても現実世界には考慮しきれない要素が多数存在します。

とはいえ技術者が取扱う自然科学の分野では「設計」という概念が成立します。最初から設計通りにはいかないにしても、設計と現実の差分を生じる原因を分析して、さらに設計に取込んでいくことで最終目標に近付くことができます。翻って人生はどうでしょう。「人生設計」という言葉はありますが、人生を設計通りに生きることはできません。人の世には世の中を構成する人それぞれの考えや思いがあって、物理法則が支配する自然界以上に複雑です。設計と現実の差分を生じる原因が分ったとしても、年齢を元に戻すことはできません。

「いくつになっても人生はやり直せる」というのは単なる気休めの呪文だと思います。駄待ち狐が「会社の選択を間違えた」と悟ったのは30代後半ですが、その時はもう家庭の事情で身動きが取れなくなっていました。会社人生のロスタイムで最後の挑戦をしようとした時には完全に戦力外扱いでした。ただ、やり直すことはできなくても、今までの経験を基に導き出した指針を残りの人生にフィードフォワードしていくことはできます。「駄待ちの人生も捨てたもんじゃなかったなあ」と笑いながら人生の最期を迎えたいものです。

複素積分の定理公式を数値積分で確認（１）に戻る