読者です 読者をやめる 読者になる 読者になる

tsujimotterの下書きノート

このブログは「tsujimotterのノートブック」の下書きです。数学の勉強過程や日々思ったことなどをゆるーくメモしていきます。下書きなので適当です。

記事一覧はこちらです。このブログの趣旨はこちら

メインブログである「tsujimotterのノートブログ」はこちら

ワイエルシュトラスの楕円関数が計算できない

ワイエルシュトラスの楕円関数を計算しています。

目標はこんな図。


この計算は,以下で定義されるワイエルシュトラスの  \wp 関数を,変数  z の実部・虚部ともに  [-4, 4 ] の範囲で計算するというものです。

 \displaystyle \wp(z; \omega_1, \omega_2) = \frac{1}{z^2} + \sum_{m^2+n^2 \neq 0}\left\{ \frac{1}{(z+m\omega_1 + n\omega_2)^2} - \frac{1}{(m\omega_1 + n\omega_2)^2} \right\}

ここで, \omega_1, \omega_2 は half-period と呼ばれる定数で,  \wp 関数の周期を決めています。

この2つの定数は, g_2, g_3 という不変量と関連していて, g_2, g_3 が決まれば  \omega_1, \omega_2 も決まります(逆も然り)。

この図では, g_2, g_3 として,以下の複素数値を使っています。

 g_2 = 1+i
 g_2 = 2-3i


さぁこれを計算しようと思うわけですが,右辺は無限級数になっていますから,残念ながらそのまま計算することは出来ません。


そこで,以下に記す2通りの作戦を考えてみたのですが,試してみると見事失敗。。。


とりあえず,過程だけでも残しておきます。

作戦1:テータ関数を使う

Wikipedia のここにあるような「楕円テータ関数」をたくさん使って計算します。

テータ関数 - Wikipedia

使うのは,

 \vartheta(z, \tau), \vartheta_{01}(z, \tau), \vartheta_{10}(z, \tau), \vartheta_{11}(z, \tau)

です。ややこしいですが,Wikipediaの式をがんばって打ち込みます。テータ関数は,指数関数の級数で表されますから,収束が早いとも書いてあったので,これはいけると思いました。


楕円テータ関数からペー関数を作るための式は,以下のWikipediaに載っています。

ヴァイエルシュトラスの楕円函数 - Wikipedia

こんな式ですね。めっちゃ複雑。

 \displaystyle \wp(z; \tau) = \pi^2\vartheta^2(0, \tau)\vartheta^2_{10}(0, \tau) \frac{\vartheta^2_{01}(z, \tau)}{\vartheta^2_{11}(z, \tau)} - \frac{\pi^2}{3} \left[ \vartheta^4(0, \tau) + \vartheta^4_{10}(0, \tau)\right]


ただし,このままでは引数が 2つになってしまうので,次のように引数 3つバージョンに補正します。

 \displaystyle \wp(z; \omega_1, \omega_2) = \frac{ \wp(\frac{z}{\omega_1}; \frac{\omega_2}{\omega_1}) }{ \omega_1^2 }


一つ難点をあげるとすると,引数が half-period である  \omega_1, \omega_2 になっていること。
本当は  g_2, g_3 を使って計算したいのですが。

まぁとにかく,  \omega_1, \omega_2 を適当に入れてやってみようと思ったのですが,失敗でした。

どう失敗したかというと,NaN が出まくったのですね。たしかに,指数関数の計算が限りなく 0 に近づいてしまったところで,その値で割り算でもしたのでしょうね。まともに数値が出なかったので,この方法はやめです。

次に行きましょう。

作戦2:ローラン級数展開

MathWorldの以下のページにローラン級数展開の式がありました。式 (4) 以降ですね。

Weierstrass Elliptic Function -- from Wolfram MathWorld

 \displaystyle \wp(z) = \frac{1}{z^2} + \sum_{k=2}^{\infty} c_k z^{2k-2}

係数は以下のように漸化式で求められると。

 \displaystyle c_2 = \frac{g_2}{20}

 \displaystyle c_3 = \frac{g_3}{28}

 \displaystyle c_k = \frac{3}{(2k+1)(k-3)} \sum_{m=2}^{k-2} c_m c_{k-m}


今度こそ良い作戦だと,思いました。なにせ, g_2, g_3 から直接計算できますからね。
しかも,やっかいな指数関数が登場しない。

よしやってみよう,気合いを入れて書いたのが,以下のコード。

require 'Complex'

# 不変量
g2 = 1 + Complex::I
g3 = 2 - 3*Complex::I


# 係数の計算
c = Array.new
c.push 0           # c0 
c.push 0           # c1 
c.push (g2 / 20.0) # c2
c.push (g3 / 28.0) # c3

for k in 4..100 do 
	ck = 0
	for m in 2..(k-2) do 
		ck += c[m]*c[k-m]
	end
	ck = 3 * ck / ((2*k+1)*(k-3))
	c.push ck
end


# ワイエルシュトラスのペー関数
def weierstrassP(z, g2, g3, c)
	sum = 1.0/(z*z) 

	for k in 2..(100) do 
		sum += c[k] * (z**(2*k-2))
	end

	return sum
end


(-4.0).step(4, 0.1) do |y|
	(-4.0).step(4, 0.1) do |x|
		t  = x + Complex::I*y

		z = weierstrassP( t, g2, g3, c)  

		puts "#{t.real},#{t.imag},#{ (z.abs.nan?) ? 1.0e+99 : z.abs }"
	end
	puts ""
end


出た結果がこれです。

f:id:tsujimotter:20150220022159p:plain:w400

一応述べておくと,横軸  {\rm Re}(\tau) ,縦軸  {\rm Im}(\tau) ,カラーチャートは  |\wp(\tau)| を表しています。


何かがおかしい。。。ある円の中だけ収束していて,その周りがすべて発散してしまっている。
も,もしかして・・・,級数の収束半径の問題ですか。

あー,解析接続しないといけないのか。。。


というわけで,まだまだ先は遠いようです。


ここもデバッグの参考になりそう。

楕円関数 - 特殊関数 グラフィックスライブラリー Special Functions

いっそ,Mathematica 使うか。。。

http://reference.wolfram.com/language/tutorial/EllipticIntegralsAndEllipticFunctions.ja.html