ryamadaのコンピュータ・数学メモ

2011-09-19

円だけど等速ではない

R 楕円関数

数学セミナー2011年10月号には、ヤコビの楕円関数の連載記事がある
三角関数 $x=\cos(t),y=\sin(t)$ は円を表す
$\frac{dx}{dt}=y,\frac{dy}{dt}=-x$ という関係にもある
xとyの追いつ追われつを表す
ヤコビの楕円関数では
- $x=sn(t,k),y=cn(t,k),z=dn(t,k)$ のようにkを母数としたtの関数が３つできる
- $\frac{dx}{dt}=yz,\frac{dy}{dt}=-xz,\frac{dz}{dt}=-k^2xy$ という関係がある
- また $x^2+y^2=1$ のままなので、 $x=\sin(f(t)),y=\cos(f(t))$ で $f(t)$ は単調な関数らしい
Rでは、ellipticパッケージがあって、ヤコビの楕円関数があるので、それで描いてみる(初めは微分方程式を離散的に追跡しようかとも思ったが)
母数を時間で変化させると、かなり変な軌道も作れることも示した

library(elliptic)
# Jacobi form of the elliptic functions
N<-10000
t<-seq(from=0,to=1,length=N)*2*pi*4
k<-0.7
x<-sn(t,k)
y<-cn(t,k)
z<-dn(t,k)
XYZ<-cbind(x,y,z)
matplot(XYZ,type="l")
plot3d(XYZ,xlim=c(-1,1),ylim=c(-1,1),zlim=c(-1,1))

k2<-1/k

x<-sn(t,k2)
y<-cn(t,k2)
z<-dn(t,k2)
XYZ<-rbind(XYZ,cbind(x,y,z))
matplot(XYZ,type="l")
plot3d(XYZ,xlim=c(-1,1),ylim=c(-1,1),zlim=c(-1,1),col=c(rep(1,N),rep(2,N)))



N<-1000
t<-seq(from=0,to=1,length=N)*2*pi*4

for(i in 1:N){
	tmpk<-i/(N+1)
	x[i]<-sn(t[i],tmpk)
	y[i]<-cn(t[i],tmpk)
	z[i]<-dn(t[i],tmpk)

}
XYZ<-cbind(x,y,z)
matplot(XYZ,type="l")
plot3d(XYZ,xlim=c