ワイエルシュトラスの楕円関数を計算しています。
目標はこんな図。
"Weierstrass elliptic function P" by Fibonacci - Own work. Licensed under CC BY-SA 3.0 via Wikimedia Commons.
この計算は,以下で定義されるワイエルシュトラスの 関数を,変数 の実部・虚部ともに ] の範囲で計算するというものです。
ここで, は half-period と呼ばれる定数で, 関数の周期を決めています。
この2つの定数は, という不変量と関連していて, が決まれば も決まります(逆も然り)。
この図では, として,以下の複素数値を使っています。
さぁこれを計算しようと思うわけですが,右辺は無限級数になっていますから,残念ながらそのまま計算することは出来ません。
そこで,以下に記す2通りの作戦を考えてみたのですが,試してみると見事失敗。。。
とりあえず,過程だけでも残しておきます。
作戦1:テータ関数を使う
Wikipedia のここにあるような「楕円テータ関数」をたくさん使って計算します。
使うのは,
です。ややこしいですが,Wikipediaの式をがんばって打ち込みます。テータ関数は,指数関数の級数で表されますから,収束が早いとも書いてあったので,これはいけると思いました。
楕円テータ関数からペー関数を作るための式は,以下のWikipediaに載っています。
こんな式ですね。めっちゃ複雑。
ただし,このままでは引数が 2つになってしまうので,次のように引数 3つバージョンに補正します。
一つ難点をあげるとすると,引数が half-period である になっていること。
本当は を使って計算したいのですが。
まぁとにかく, を適当に入れてやってみようと思ったのですが,失敗でした。
どう失敗したかというと,NaN が出まくったのですね。たしかに,指数関数の計算が限りなく 0 に近づいてしまったところで,その値で割り算でもしたのでしょうね。まともに数値が出なかったので,この方法はやめです。
次に行きましょう。
作戦2:ローラン級数展開
MathWorldの以下のページにローラン級数展開の式がありました。式 (4) 以降ですね。
係数は以下のように漸化式で求められると。
今度こそ良い作戦だと,思いました。なにせ, から直接計算できますからね。
しかも,やっかいな指数関数が登場しない。
よしやってみよう,気合いを入れて書いたのが,以下のコード。
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
出た結果がこれです。
一応述べておくと,横軸 ,縦軸 ,カラーチャートは を表しています。
何かがおかしい。。。ある円の中だけ収束していて,その周りがすべて発散してしまっている。
も,もしかして・・・,級数の収束半径の問題ですか。
あー,解析接続しないといけないのか。。。
というわけで,まだまだ先は遠いようです。
ここもデバッグの参考になりそう。
いっそ,Mathematica 使うか。。。
http://reference.wolfram.com/language/tutorial/EllipticIntegralsAndEllipticFunctions.ja.html