結局、2ヶ月ほども掛けて出てきたのはありがちなグラフでした。ヤコビの楕円関数、sn(u)、cn(u)、dn(u)の実数部分を表計算ソフトで表示したので公開します。
やりたいことは複素平面での表示で、これはその1段階前の計算です。このグラフの形は数学愛好家なら既視感があるはずで、多分、ネットで見られると思います。
(参考文献1) 森口繁一他、岩波 数学公式 III、岩波書店 1960-1997
(参考文献2) 戸田盛和、楕円関数入門、日本評論社 2001
(参考文献3) 泉信一他、共立 数学公式 附・函数表 改訂増補、共立出版 1953-1969
楕円関数には母数(modulus) kという、元の楕円の形を決める係数があり、0≦k≦1で、第二種楕円積分では楕円の離心率(e)と同じです。0だと真円になり、数値を上げると縦方向が縮んで行き、1だとぺちゃんこになります。離心率は長軸上の焦点の位置でもあります。
ヤコビの楕円関数は第一種楕円積分の話なのですが、グラフを見ると何となく楕円上を周回する点の座標と関係しているような気がします。
この作図は参考文献1の50ページに載っている、楕円テータ関数を仲介する方法で、その4種の楕円テータ関数のグラフが下方に見えます。ほぼ正弦波に見えますが、わずかに歪んでいます。
セルA4のk2は母数の二乗の値で、ここを書き換えるとグラフの形が変わります。算出方法の制限から0.01≦k2≦0.99の値で試してみて下さい。その他の数値も受け付けますが、誤差がもろに出てきます。
4行目がグラフの計算に必要な定数値の計算です。
セルB4 =(1-(1-A4)^0.25)/2/(1+(1-A4)^0.25) は参考文献にあるl値で、その隣の、
セルC4 =B4+2*B4^5+15*B4^9+150*B4^13+1707*B4^17 のq値を計算するための中継地です。
このq値は本来は二重周期のKとK'から算出される値で、そのK,K'は第一種完全楕円積分で計算できます。つまり、k→K,K'→qの順で算出するのが正統ですが、第一種完全楕円積分(定積分)を数値計算する必要があります。これを避けて大幅に計算量を減らす工夫がl値の採用です。ただし、その代償としてk2=1に極端に近い値は採れません。1に近い部分はとても興味深い結果が得られて、参考文献2で詳しく紹介されていますが、その再現はできません。
また、k2=0だと本来は正弦波が出てくるはずですが、試してみると分かりますがそうはなっていません。0.01等の小さい値を入れる必要があります。
セルD4 =PI()/2*(1+2*(C4+C4^4+C4^9+C4^16+C4^25))^2 はそのK値です。グラフには使いませんが、本来は周期(4K)が延びるのにこのグラフでは周期を固定しているので目安として計算しました。k=0の時、Kはπ/2 (=1.570796...)となります。
セルE4 =2*PI()*(C4^(1/4)-3*C4^(9/4)+5*C4^(25/4)-7*C4^(49/4)+9*C4^(81/4)) は楕円テータ関数θ1(v)の微分関数のv=0での値です。これも今回のグラフでは使いませんが、他の楕円関数で必要になる場合があるので参考として表示しています。
セルF4 =2*(C4^(1/4)+C4^(9/4)+C4^(25/4)+C4^(49/4)+C4^(81/4))
セルG4 =1+2*(C4+C4^4+C4^9+C4^16+C4^25)
セルH4 =1+2*(-C4+C4^4-C4^9+C4^16-C4^25)
はそれぞれ、楕円テータ関数θ2(0)、θ3(0)、θ0(0) (=θ4(0))の値です。Kやθ'1(0)ともども、本来は無限級数ですが、グラフ表示に必要な精度の最初の5項で打ち切っています。
セルA10～セル110にはあらかじめ-50～50の1刻みの数値を入れています。これがグラフの横軸になります。
セルB10 =$B$8+$C$8*A10 はそれぞれの楕円関数の入力値vの計算で、これを下に続くセルB11～B110にコピーします。セルB8にインデックス0の時の値を入れ、セルC8に刻み値を指定します。
セルC10 =2*($C$4^(1/4)*SIN(PI()*B10)-$C$4^(9/4)*SIN(PI()*3*B10)+$C$4^(25/4)*SIN(PI()*5*B10))
セルD10 =2*($C$4^(1/4)*COS(PI()*B10)+$C$4^(9/4)*COS(PI()*3*B10)+$C$4^(25/4)*COS(PI()*5*B10))
セルE10 =1+2*($C$4*COS(2*PI()*B10)+$C$4^4*COS(4*PI()*B10)+$C$4^9*COS(6*PI()*B10))
セルF10 =1+2*(-$C$4*COS(2*PI()*B10)+$C$4^4*COS(4*PI()*B10)-$C$4^9*COS(6*PI()*B10))
には4種の楕円テータ関数、θ1(v)、θ2(v)、θ3(v)、θ0(v)の計算式が入っていて、
これをセルx11～x110にコピーします。本来はこれも無限級数です。最初の3項しか計算していませんが、ここでのグラフ描画には十分な精度のはずです。
セルG10 =$G$4/$F$4*C10/F10
セルH10 =$H$4/$F$4*D10/F10
セルI10 =$H$4/$G$4*E10/F10
が、今回の目的のヤコビの楕円関数の計算式です。これもセルx11～x110にコピーします。
あとは、必要なセルを選んで、メニューの挿入→折れ線グラフを指示すればグラフが描かれます。
Sheet1は示しませんが、Sheet2の行4の定数値の、k2が0.00～0.99の0.01刻みの表です(コピーして伸ばしたものです)。楕円テータ関数に興味のある方は作って眺めておくと良いと思います。
sn(u)は三角関数のsin(u)に相当し、cn(u)はcos(u)に相当します。dn(u)に相当物はありませんが、双曲線関数のsech(u)の繰り返しになっているそうです。sech(u)は懸垂線のcosh(u)の逆数の関数です。
