Schwarz-Christoffel formula¹ とは、複素平面上で単位開円板を多角形の上へ等角写像する関数の形を定めた公式です。今回これを Python で実装してみたのですが、実装の際に関数の多価性に起因してひとつ引っかかったポイントがあったので記事にしておこうと考えました。

さて、件の公式について具体的には次の定理の主張が成り立ちます。

定理

$w$-平面上の単位開円板 $|w|<1$ を、$z$-平面上の角 $a_k \pi\ (k=1,\dots,n)$ の多角形 $\Omega$ の上へ等角写像する関数 $F(w)$ は

F(w) = C_0 \int_0^w \prod_{k=1}^n (w' - w_k)^{-b_k} dw' + C_1

の形である。ただし、$b_k = 1 - a_k$（外角）で、$w_k$ は単位円周上の点、$C_0, C_1$ は定数である²。

証明の概略

ここで述べることは証明の概略（といえるかどうかも怪レい代物）に過ぎないので、詳しくは参考文献²を参照されてください。

Riemann の写像定理により、複素平面内の$\mathbb{C}$ と異なる任意の単連結領域は単位開円板と解析的同型であることが知られています。これより多角形 $\Omega$ を単位開円板に写す解析的同型写像 $f$ が存在し、とくに $f$ は等角写像となります。そこで $f$ の逆関数を $F$ とおきます。すると斯々然々の議論により $H(w) = F'(w) \prod_{k=1}^n (w - w_k)^{b_k}$ は単位閉円板の或る開近傍上で解析的かつ常に $\neq 0$ となります。するとこの近傍上で $\log H(w)$ の1価正則な枝がとれて、したがって $\Im \log H(w) = \arg H(w)$ は調和となります。そこで調和関数の最大値原理を適用するなどして、実は $H(w)$ は定数であることが示せます。したがって $F'(w) = C_0 \prod_{k=1}^n (w - w_k)^{-b_k}\ (\exists C_0 \in \mathbb{C})$ となります。これを積分して所期の公式を得ます。

実装

Google Colab 上で動くコードは次の通りです。コード後半はほとんど図示のための処理なので、定理の公式に関係するのは前半、とくに関数 dF や F の部分です。

import matplotlib.pyplot as plt
import numpy as np
import cmath
import scipy.integrate as integrate

# provide constants and parameters
PI = cmath.pi
C0 = 1
C1 = 0

wts = [PI*1/6, PI*5/6, PI*9/6]
ws = list(map(cmath.rect, np.ones_like(wts), wts))
bs = [2/3, 2/3, 2/3]
L = range(len(ws))

assert len(ws) == len(bs)
assert np.isclose(sum(bs), 2)

# define functions and data
def dF (w):
    rs = [abs(w - ws[k]) for k in L]
    # direct a branch cut towards outside the unit circle
    ts = [(cmath.phase(w - ws[k]) - wts[k]) % (2 * PI) + wts[k] for k in L]
    # ts = [cmath.phase(w - ws[k]) % (2 * PI) for k in L] # doesn't work :(
    prod_r = np.prod([rs[k] ** (-bs[k]) for k in L])
    sum_t = np.sum([ts[k] * (-bs[k]) for k in L])

    return prod_r * cmath.exp(1j * sum_t)

def F (w):
    re = integrate.quad(lambda t: (dF(w * t) * w).real, 0, 1)[0]
    im = integrate.quad(lambda t: (dF(w * t) * w).imag, 0, 1)[0]
    return C0 * (re + 1j * im) + C1

theta = np.sort(np.concatenate((np.linspace(0, 2 * PI, 200), wts)))
W0 = np.cos(theta) + 1j * np.sin(theta)

# define subplots
fig, axs = plt.subplots(1, 2, figsize=(16, 8))
axs = axs.flat

plt.suptitle('Schwarz-Christoffel formula', fontsize=16, y=0.94)
axs[0].set_xlabel('Re')
axs[0].set_ylabel('Im')
axs[0].set_aspect('equal')
axs[1].set_xlabel('Re')
axs[1].set_ylabel('Im')
axs[1].set_aspect('equal')

# draw scatter plots for each radius
for r in np.linspace(0, 1, 10):
    W = r * W0
    Z = np.array(list(map(F, W)))
    
    axs[0].scatter(W.real, W.imag, c=theta, s=5, alpha=r)
    axs[1].scatter(Z.real, Z.imag, c=theta, s=5, alpha=r)

# put labels on each vertex
for j in L:
    to = F(ws[j])
    axs[0].text(ws[j].real, ws[j].imag, f'$w_{j}$', fontsize=16)
    axs[1].text(to.real, to.imag, f'$F(w_{j})$', fontsize=16)

plt.show()

出力結果：

実装のポイント

上のコードは定理の公式を素朴に実装したものですが、実装に際して一箇所だけ注意すべきポイントがありました。それは関数 dF 内で # direct a branch cut towards outside the unit circle とコメントを付けている行、すなわち各 $k$ に対して $(w - w_k)^{-b_k}$ の適当な分枝をとる部分です。

Python の cmath ライブラリでは偏角関数の branch cut が負の実軸に伸びている³ため、べき演算を素朴に実装すると $w - w_k$ が負の実軸を通過するタイミングで $(w - w_k)^{-b_k}$ の（したがって $F'(w)$ の）値が不連続に変化してしまいます。そのようなまずい実装を参考までに件の行の直後にコメントアウトしてあり、その出力結果は次のようになります。

失敗例：

三角形がめくれてしまっているのがわかります。この現象は、左の $w$-平面において、点 $w$ が線分 $w_0$ to $w_1$ を上下に通過するタイミングで $(w - w_0)^{-b_0}$ の値が不連続に変化してしまうことが原因です。

そこで解決策として、上のコードのように各 $k$ に対して $(w - w_k)^{-b_k}$ の値を計算する際に branch cut を円の外側⁴に伸ばすことで対応してみました。これで不連続性の問題は解決され、無事に成功例のような図が出力されるようになります。

最後に

Schwarz-Christoffle formula に限らず、ふだん諸々の証明をするときには適当な分枝の存在さえわかればよい場合が多くて、あまり具体的な分枝を求める努力をしたことがなかったような気がします。しかし今回数値計算を行うにあたり、具体的に分枝を与えなければならない場面に遭遇し、せっかくだから記事に残してみようと思った次第です。おそらくきちんとした実装では（もっとよいやり方で）対策されていると思うので、ベストプラクティスは他の文献をあたっていただくとして、ここではとりあえず問題点とその簡易な解決策についての共有でした。

注釈

1. https://en.wikipedia.org/wiki/Schwarz%E2%80%93Christoffel_mapping ↩

2. L. V. Ahlfors. 複素解析. 笠原乾吉 訳 ↩ ↩²

3. https://docs.python.org/ja/3/library/cmath.html ↩

4. つまり、$w$ の関数 $(w - w_k)^{-b_k}$ は $w_k$ が分岐点となっていて、$w_k$ から $\arg w_k$ の向きに branch cut をとって分枝を定めるということです。 ↩