ヤコビの楕円関数のグラフを書こう！！ [数値解析]
ヤコビの楕円関数のグラフを書こう！！
ヤコビの楕円関数
あるいは、
のグラフは、表計算ソフトがあれば、簡単にかけてしまう。
ヤコビの楕円関数は、（１）のようにF(k,φ)と書くのが正式なんだけれど、母数kは0<k²<1を満たす定数だから（２）と書いても混乱は生じないだろう。
ということで、ヤコビの楕円関数をF(φ)で表すことにする。
（２）をφで微分すると、
という微分方程式が得られる。
この微分方程式の初期条件は
というわけで、１階の常微分方程式の初期値問題に書き換えることができて、オイラー法やルンゲ・クッタ法でこれを数値積分することによって、ヤコビの楕円積分の（近似）値を求めることができる。
この数値積分は最も簡単なオイラー法を用いれば十分で、
と計算できる。
特に、
と分割の幅が等間隔の場合、
と計算できる。
ここで、
ということで、表計算ソフトを用いて、h=π/10、母数k=0.5として計算した結果は次の通り。
k=0.5のときは、ヤコビの楕円関数はほとんどF(φ)=φと直線で近似できてしまうことがわかる。
だって、このヤコビの楕円関数は振り子の等時性に深く関係している。これが直線で近似できないと、振り子時計なんてものは生まれなかった！！
kを大きくし、k=0.9にしたとき、ヤコビの楕円関数のグラフは次のようになる。
これくらいkを大きくすると、ヤコビの楕円関数は直線からずれてくるけれど、それでも直線的な振る舞いをすることがわかるにゃ。
グラフにすると、この関数は意外につまらない。
常微分方程式の初期値問題の数値解法に関する記事はこれまでにずいぶんと書いたので、そろそろ、これくらいの計算ならば自力で計算できるようになってほしいと思うのだけれど、この計算に使用したスプレッドシートを公開したので、これを使って、あれこれと遊んで欲しい。
ネムネコと同じように、くたばれマイクロソフト、アップルのペンギンOS仲間は、ods形式のコチラがいいと思う。
このスプレッドシートは、kの値を変化させると、それに合せて、再計算、グラフを再描画してくれるにゃ。
k=0.7に変更すると、
ということで、恒例の自画自賛ソングを♪ これは重要な儀式だから省くわけにはいかない！！
上の計算ではh=π/10としているけれど、 この値を変えても、それに合わせて、再計算、グラフを再描画してくれる。お前らにできるだけ負担がかからないように、配慮してあるにゃ。
そして、ネムネコが考えるに、
を計算し、このグラフをかいてくれるスプレッドシートは、ddt³さんが作ってくれるんじゃ〜ないか。
たぶん、ddt³さんは台形公式を使って所定の精度で計算できる素敵なスプレッドシートを作ってくれると思う。
k=0.9のとき、単純なオイラー法の場合、h=π/10という粗い分割だと、誤差が大きすぎて、
とおいたとき、
という関係を満たさない(^^ゞ
オイラー法だと 
だから、（４）という数学的な関係を満たさない(^^ゞ
近似計算だとは言え、こういった基本的な数学的関係を満たさないものはまずい。
お昼の前にバタバタと１分ほどで作ったスプレッドシートで、計算結果の検討が足りなかった。これくらいの計算ならばオイラー法で十分だろうという見込みが甘かった。
「んっ」と、このことに気づいた時には、スプレッドシートと記事をアップしていたので、すべては後の祭りであった。
この動画のフレンズのように、ここで倒れるわけにはいかない。
そこで、積分の中点公式を使って、
と近似し、漸化式⑨を
とした改良した方がいいみたいですね。
⑨³は２次ルンゲ・クッタ法に相当するもので、単純なオイラー法よりも精度がよく、h=π/10という粗い分割でも、
となって、
という関係を満足する。
こちらの方がいいね。スプレッドシートの修正も必要最小限で済むし・・・。
より高精度の４次のルンゲ・クッタ法で計算した結果は、K=2.2805247と、２次のルンゲ・クッタ法の計算結果とほとんど同じ値なので、この近似計算は十分に信用に足るものと考える。
k=0.9のとき、K≒2.280549くらいらしいから。
だ・か・ら、
あらためて
というわけで、改良版のスプレッドシートを改めて公開！！
もっと、本格的な数値計算をしたいという物好きのために、ルンゲ・クッタ法を用いたFortranプラグラムを紹介するにゃ。
parameter(n=20)
real x(0:n), y(0:n)
pi=acos(-1.0)
h=2*pi/n
do i=0, n
x(i)=i*h
end do
y(0)=0.
do i = 1, n
dk1=h*f(x(i-1), y(i-1))
dk2 = h*f(x(i-1)+h/2, y(i-1)+dk1/2)
dk3 = h*f(x(i-1)+h/2, y(i-1)+dk2/2)
dk4 = h*f(x(i-1)+h, y(i-1)+dk3)
y(i)=y(i-1)+(dk1+2*dk2+2*dk3+dk4)/6.
end do
do i=n,1,-1
write(*,100) -x(i), -y(i)
end do
do i = 0 , n
write(*,100) x(i), y(i)
end do
100 format(f10.7,1x, f10.7)
end
function f(x,y)
real k
k=0.9
f=1./sqrt(1-k**2*(sin(x))**2)
end
function文のkの値を変えるといいにゃ。
2018-09-06 12:59
nice!(1)
コメント(0)
コメント 0