13 Carlsonの対称楕円積分

このところやっている数値計算ネタで、Mac とは直接的な関係のない話。

関数論のうちであまりよく知らないところの復習をしていて、 学生の頃から気になっていた楕円積分と対峙することとなった。

その数値計算について、知っていてしかるべきだったのに、全く知らなかった、 ということがあったので、自分の整理のためにも簡単にメモしておく。

まず微積分の復習から。 2変数の有理関数 $R(t,s)$と多項式 $p(x)$ に対して、不定積分 $\dsp\int R(x,\sqrt{p(x)})\Dx$ は、$p(x)$ の次数が$2$以下だったら初等関数の範囲で 求めることができる(微積分ネタとなる)が、次数が3以上の場合は、 特殊な例外をのぞいて求まらない。 特に $p(x)$ の次数が 3,4の場合は楕円積分と呼ぶ。 学生の頃の感想としては、 $\dsp\int R(x,\sqrt{ax^2+bx+c})\Dx$ の段階で、結構面倒な計算をやらされて、 正直ちょっと嫌だけれど、その先があるのかぁ、大変だなあ。

数学に少しはまりかけると、高木貞治「近世数学史談」とか読んで、 何か良く分からないけれど、すごそう、大変だったんだなあ、ということが知れる。 Legendre という数学者が楕円積分について記念碑的な著書を書いたけれど、 (楕円積分の逆関数である) 楕円関数の研究という大躍進はその後にやってきたとか。

授業でも、楕円関数や Riemann 面について教わる訳だけれど、 そういうのを何かに使うというのは経験がなかった。 時々、これは楕円積分、楕円関数を使って扱えるというのに遭遇した訳だけど、 勉強するコストを考えて後回しにしてきた面がある。

定積分の値が知りたいだけならば、数値積分という手がある。特異性があっても、 二重指数関数型数値積分公式とか使うとうまく計算できる、 というようなのは講義でも話したりする話題だ。 そういう立場からは、扱った例の中に楕円積分に分類されるものもありますね、 くらいのものだ。

今回は、与えられた領域の等角写像やその逆関数を計算する問題を考えていて、 「有名な」Schwarz-Christoffel 変換を初めてちゃんと学ぶこととなった。 カッコ付きなのは、取り上げていない関数論のテキストも結構あって、 取り上げていても証明略とか、説明がわかりにくいとかいうのが多かったので、 実は顔は知っているけれど、どんな人かはよく存じ上げません、 という感じ？

その勉強もどうやら一段落して (とても面白かった)、 最後にやはり上半平面 $H=\left\{z\in\C\relmiddle\vert \MyIm z>0 \right\}$ から長方形領域への双正則関数を求めるという定番問題の解を数値計算してみよう、 となった。 問題の設定のしかたによっては

$$F(z;k):= \int_0^z\frac{\D\zeta}{\sqrt{(1-\zeta^2)(1-k^2\zeta^2)}}$$

を計算することとなる。これはとても有名な、 Legendreの第1種楕円積分の標準形というものである。

それを計算するためのライブラリィがあったはずだ。えーと、、、 Mathematica だと EllipticF[$\phi$,m] というのがそれか。 これは

$$F(\phi,m):=\int_{0}^{\phi}\dfrac{\D t}{\sqrt{1-m\sin^2t}}$$

を計算する関数で、$F(z;k)$ とは

$$F(\sin\phi;k)=F(\phi,k^2)$$

という関係にある、と。記号が大変紛らわしいけれど、 これは歴史的経緯で仕方がないことのようだ。 そうすると、

$$F(z;k)=F(\arcsin z,k^2)$$

で計算すると良いのか。$\arcsin$ は多価関数だから、少し気持ちが悪いな。 まあやってみるか。

まずは論より run で、

F[z_,k_]:=EllipticF[ArcSin[z],k^2]

と定義して、曲線 $F[x,k]$ ($x\in\R$) を描いてみよう。 もし関数定義が正しければ、長方形の周を一周する曲線であるはずだ ( $F(-\infty;k)=F(+\infty;k)$ に注意)。

ParametricPlot[ReIm[F[t,0.5]],{t,-10,10}]

図: 長方形の周になっているかな？

あれ？$t=1$ で左折して(上に向かって)ほしいのに右折している。 方向音痴だ。 まあこういうことは覚悟していた。 手短に原因を説明すると、ArcSin[$z$]は、 $-i\Log[iz+\sqrt{1-z^2}]$ ($\sqrt{}$ はいわゆる主値) と定義してあって、 $(-\infty,-1]$ と $[1,+\infty)$ が分岐截線 (branch cut) だけれど、 右側の線上では下半平面と連続につながるようになっているから、である。

$\arcsin$ は、私の頭の中では $\dsp \int_0^z\dfrac{\D\zeta}{\sqrt{1-\zeta^2}}$ が第一近似で、 これは円弧二角形 $\left\{x+iy\relmiddle\vert-\pi/2\le x\le \pi/2,\ 0\le y<+\infty \right\}$ の Schwarz-Christoffel 写像と解釈できるから、 上半平面に制限して連続になってくれることを期待するけれど (Jordan領域の間の双正則関数は、閉包まで同相写像に拡張できるという、 Carathéodory の定理を思い出す)、 ライブラリィ関数としては、 $-i\Log[iz+\sqrt{1-z^2}]$ と定義してあるので、残念ながらそうはなっていないと。 まあ、奇関数にならないのも変な感じがするし、仕方がないのかな。 自作の myArcSin[] (上半平面に制限すると連続) で置き換えよう。

as[x_] := Which[x >= 1, Pi/2 + I  Log[x + Sqrt[x^2 - 1]], True, ArcSin[x]]
myArcSin[z_] := If[Im[z] != 0, ArcSin[z], as[z]]
F2[z_, k_] := EllipticF[myArcSin[z], k^2]
gcurve2 = ParametricPlot[ReIm[F2[t, 0.5]], {t, -5, 5}]

これでどうなる？

図: 長方形の周になっているかな？雪辱戦

うわあ。$t=1$ では無事左折するようになったけれど、 今度は $t=1/k$ で右折してしまった… これは $\arcsin$ のせいではなく、 EllipticF[] の値の選び方のせいである。

こちらも分岐截線のところで、崖の下と上とどちらを選ぶかの話で、 ArcSin[] の方は崖の上の方を選ぶべきだったけれど、 今度は崖の下の方を選ぶべきと。

もうひと頑張り… まあ、なんとか出来ました。どうやったかは内緒 (ちょっとその場しのぎの技を使いました)。

図: 長方形の周になっているかな？これでどうでしょう。

今は便利なツールがあって、なぜ考えていたようにならないのか、 すぐ画面で見えるようになっているから、1日の間に解決したけれど、 そうでなければ大変だったろうと思う。

意外と使いにくいな EllipticF[]。

数学のテキストを見ると $F(z;k)$ というのは良く出て来るけれど、 それをコンピューターで用意されている $F(\phi,m)$ (EllipticF[]) を使って計算するのは大変だ、 ということだ。 何か私の知らないうまい手があるのかしら？

それで見つけたのが、今回のお題の Carlson の対称楕円積分である。 たくさん種類があるけれど、代表的なのは

$$R_F(x,y,z):=\frac{1}{2} \int_0^{+\infty}\frac{\D t}{\sqrt{t+x}\sqrt{t+y}\sqrt{t+z}}.$$

Legendre の関数 $F(z;k)$ は、分母の零点が $\pm1$, $\pm 1/k$ で、 それを $k$ という1つのパラメーターで表したけれど、 $R_F(x,y,z)$ は分母の零点をそのまま ($-1$ はかけてあるか) $x$, $y$, $z$ と与えている。 変数の数が増えたけれど、対称性がある。

$R_F$ の他にどういう関数が用意されているかは、Wolfram MathWorld の

https://mathworld.wolfram.com/CarlsonEllipticIntegrals.html

が手取り早い。

ちなみに $R_F$があれば $F(\phi,m)$ は

$$F(\phi,m)=\sin\phi R_F(\cos^2\phi,1-m\sin^2\phi,1)$$

で計算できる。

Carlson の楕円積分は、登場したのはそんなに古くはないけれど、 新しくもない。1977年の論文が最初とされているが、 それ以前からの研究の延長上にあるので、文献リストは結構長い (以下は1963年からの論文と著書である -- [#!Carlson1963!#], [#!Carlson1964!#], [#!Carlson1977!#], [#!Carlson1977a!#], [#!Carlson1978!#], [#!Carlson1979!#], [#!Carlson-Notis1981!#], [#!Carlson1987!#], [#!Carlson1988!#], [#!Carlson-Gustafson1994!#], [#!Carlson1995!#], [#!Carlson2020!#])。

これらの論文は “読解中” である。 自分の問題にとにかく応用しよう、 という観点からは、 第三者による解説 Press-Teukolsky [#!Press-Teukolsky1990!#] を最初に読むのが分かりやすいと思われる。

実は Carlson の対称楕円積分は、知っている人には常識的なことらしい。 Mathematica にも Python にもそれを計算するための関数が用意されている ($R_F$ は、Mathematica の場合 CarlsonRF[], Python の場合 scipy.special.elliprf())。

以下は、 特殊関数についての有名な文献である Abramowitz-Stegun [#!Abramowitz-Stegun!#] の “現代化” といわれている NIST Handbook の§19.32 Conformal Map onto a Rectangle に載っていることである。

$x_1>x_2>x_3>0$ とするとき、

$$f(z):=R_F(z-x_1,z-x_2,z-x_3)$$

とおくと、 $f$ は上半平面 $H$ を $0=f(-\infty)=f(+\infty)$, $f(x_1)$, $f(x_2)$, $f(x_3)$ を頂点とする長方形の内部へ双正則に写す。

x1 = 3; x2 = 2; x3 = 1;
f[z_] := CarlsonRF[z - x1, z - x2, z - x3]
ParametricPlot[ReIm[f[t]],{t,-10,10},PlotPoints->500]

図: Carlson の $R_F$ による上半平面から長方形への写像 (面倒なので周のみ)

あっけないくらい簡単ですね。まあ変数が $z$ である分、 $F(z;k)$ に近いから、ということはできるか。

さて、私は古い人間なので、CやC++でガリガリとプログラムを書きたい。 それ以外にも、一体どうやって計算するのかに興味がある。 どこかにプログラムないかな… Boost にあるけれど、実変数専用だ (またかいな)。

仕方がないので、自分でなんとかしよう。 Legendre の楕円積分や Jacobiの楕円関数は、計算の仕方も非常に面白いけれど、 Carlson の楕円積分の計算も大変面白い。

上で紹介したPress-Teukolsky [#!Press-Teukolsky1990!#]に、 単精度実数版のFortran コードが載っているけれど、非常にコンパクトである。 そしてそれはほぼそのまま複素数版に書き換えることができる。

次は (内容は全く保証しないけれど)、C言語で $R_F$ を計算する関数である。

cdrf.c

/Users/mk/.tex-inputs/knowhow-2025-fig/cdrf.c

短くてシンプルなアリゴリズム。とても美しい (要点は[#!Press-Teukolsky1990!#]に書いてある)。

Carlson [#!Carlson1995!#] の、 第3節 Numerical checks の節にある数値例を試すプログラムをつけておく。

cdrf.c

/Users/mk/.tex-inputs/knowhow-2025-fig/testcdrf.c

コンパイル&実行

% gcc testcdrf.c cdrf.c
% ./a.out
I=0+1i
count=5
RF(1,2,0)=1.311028777146060+0.000000000000000 i, 誤差=4.440892e-16
count=6
RF(i,-i,0)=1.854074677301372+0.000000000000000 i, 誤差=4.791698e-16
count=5
RF(i-1,i,0)=0.796125865842339+-1.213856669836496 i, 誤差=2.482534e-16
count=4
RF(2,3,4)=0.584082841677152+0.000000000000000 i, 誤差=0.000000e+00
count=5
RF(i,-i,2)=1.044144565406436+0.000000000000000 i, 誤差=2.279405e-16
count=5
RF(i-1,i,1-i)=0.939120502186194+-0.532962520186353 i, 誤差=2.220446e-16

(いい加減なプログラムなのだけど、すごいなあ。)

(実は ChatGPT が「私が作りましょうか？」と言ってきたので、 任せてみたらやっぱり時間を浪費したとか (全然理解できていなくて、よくそれでそう言ってくるなあ、 とちょっと呆れた)、他にも色々な話はあるけれど、 すでに相当長くなったので、そういう枝葉は省略する。 カールソンって、子供の頃に「屋根の上のカールソン」というのがあったので、 珍しくない名前と思っていたけれど、 そういえば具体的な人の名前で出会ったのは初めてだなあ。)

(2025/6/14追記) Carlson は 1963 年に $R$ 関数という多変数超幾何級数を導入した (1893年の Lauricella による 3 変数超幾何級数 $F_D$ をもじって一般化したもの)。 それはある条件下で積分表示できて、 そのうちで楕円積分であるものが対称楕円積分ということだ。 Carlson の紹介は https://math.nist.gov/opsf/personal/billecarlson.pdf(これは亡くなったときに、 OP-SF NET - SIAM Activity Group on Orthogonal Polynomials and Special Functions の電子版ニュースレター -- に載ったものらしい) が分かりやすい。

青本・喜多 [#!_______________1994!#] には、 Appell-Lauricellaの超幾何関数のことは色々書いてあるけれど、 Carlson のことは文献表に [#!Carlson1977!#] が載っているだけだ。