遺伝学・遺伝統計学関連の姉妹ブログ『ryamadaの遺伝学・遺伝統計学メモ』
京都大学大学院医学研究科ゲノム医学センター統計遺伝学分野のWiki
講義・スライド
医学生物学と数学とプログラミングの三重学習を狙う学習ツール
駆け足で読む○○シリーズ
ぱらぱらめくるシリーズ
カシオの計算機
オンライン整数列大辞典

結び目をパラメタ化

R 偏微分 複素数
  • 結び目は２次元のわっかを３次元に置いた状態
  • だから、「閉じた曲線」として表せる
  • 閉じた曲線については、ここ数日、やってきた通りで、パラメタを置くと、周期的に同じ位置に戻ってきてしまうような座標定義のこと
  • また、曲線の位置を陽に表さずに、偏微分で表す方法もあった
  • 三葉結び目を考えてみる(Mathematicaのウェブレファレンスの図は次の通り)

http://reference.wolfram.com/mathematica/ref/Files/KnotData.ja/O_1.gif
http://www.genome.med.kyoto-u.ac.jp/StatGenet/lectures/2010/Sanyou.png

  • これの特徴はなんだろう
    • ３方向について平等(対称)
    • ３方向をx,y,zとして、その一つを考えた上で、位相を\frac{1}{3}*2\piずつずらしてやればよいので、xのみを考える
    • x方向は、図の縦軸とする
    • xについては、移動の方向は、増える方向、減る方向、増える方向、減る方向、と方向を変えた上でもとに戻る
    • 方向転換点は、左右中央に２か所、中央からはずれた左右対称な位置に２か所ある
    • 移動速度は、結び目を一方方向に進んでいるとすれば、常に正
    • その速さの絶対値は大きい時と小さいときがあるが、方向転換する４か所のうち、左右中央の２か所のその上下中央よりの点を境に、増減が逆転する。左右中央の２か所の上側の点でも、増減が逆転している
    • つまり、(方向増・速さ絶対値増)→(方向減・速さ絶対値増)→(方向増・速さ絶対値減)→(方向減・速さ絶対値減)の４相がある
    • ２つの異なる周期があって、それが２対１の周波数になっていることがわかる
  • \frac{dx}{dt}=(cos(t)+k)*cos(2t)は、このような位置変化に対応する
    • ここで、kが1より大なら、速さの絶対値が正なので、条件を満たす
  • このようにパラメタ化してみると、\frac{dx}{dt}を速さの絶対値項と方向項とを三角関数で表すことで、それらに周期性をもたせるときには、２つの項の周期(の比)と、速さの絶対値項の大小を決めるk（速さの絶対値項と言ったが、この項が０になったり負になったりしてもよいようにkの定義域を拡張してやると、より多様な曲線に対応付けすることができるので、そうしてしまおう)
  • y,zの相をずらせば
  • \frac{dy}{dt}=(cos(t+\frac{1}{3}2\pi)+k)*cos(2(t+\frac{1}{3}2\pi))
  • \frac{dz}{dt}=(cos(t+\frac{2}{3}2\pi)+k)*cos(2(t+\frac{2}{3}2\pi))
  • 掲載図は、結び目ができるようにk=1.5とし、２つの周期をいろいろに変えて描いたもの
  • これをRで書けば
library(rgl)
# 周期は0-(2*pi)をX回繰り返しとする
X<-5
# １周あたりの点の数をnpt
npt<-1000
t<-seq(from=0,to=1,length.out=npt)*2*pi*X

# x,y,zの座標は、直径の周期C1と角の周期C2とで決める
# 2piの間にC1回、C2回繰り返しさせる
#C1とC2の最小公倍数*2piだけ時間が経過すると、「元に戻る」
# C1=C2=0のときは直線
# C1=0,C2!=0のときは円
# C1=C2>0のときはらせん
# C1*2=C2 のときは三葉結び目
# C1*k1=C2*k2 ;k1,k2は異なる自然数 のときは複雑な結び目
# C1とC2に最小公倍数がないと、周回して面をなす
C1<-1
C2<-2.5
# x,y,zは相互に対等なので、0/3*2pi, 1/3*2pi,2/3*2pi のずれを入れる
# 直径は周期的で (k-1)と(k+1)の間を上下するものとする
# k=1では三葉結び目は結ばれずに１個の集中点をもち
# k<1では結ばれない
# k>1で、結ばれる
k<-1.5

x<-y<-z<-rep(0,length(t))
# 初期値は適当
x[1]<-runif(1)
y[1]<-runif(1)
z[1]<-runif(1)

for(i in 2:length(t)){
	x[i]<-x[i-1]+(cos(C1*(t[i-1]+0))+k)*cos(C2*(t[i-1]+0))
	y[i]<-y[i-1]+(cos(C1*(t[i-1]+2/3*pi))+k)*cos(C2*(t[i-1]+2/3*pi))
	z[i]<-z[i-1]+(cos(C1*(t[i-1]+4/3*pi))+k)*cos(C2*(t[i-1]+4/3*pi))
}

plot3d(x,y,z,col=rainbow(1000),cex=0.1)
ylim<-c(min(x,y,z),max(x,y,z))
plot(x,type="l",col=1,ylim=ylim)
par(new=TRUE)
plot(y,type="l",col=2,ylim=ylim)
par(new=TRUE)
plot(z,type="l",col=3,ylim=ylim)

# いろいろなC1,C2で描いてみる
#C1s<-0:4
#C2s<-C1s
#C12s<-expand.grid(C1s,C2s)
C12s<-matrix(0,9,2)
C12s[1,]<-c(0,0) # 直線
C12s[2,]<-c(0,1) # 円
C12s[4,]<-c(1,1) # らせん
C12s[5,]<-c(1,2) # 三葉結び目
C12s[6,]<-c(1,4) # 複雑な結び目
C12s[7,]<-c(1,2.1) # 複雑な結び目
C12s[8,]<-c(1,sqrt(2)) # 元に戻らず面をなす。結び目はできない?
C12s[9,]<-c(1,sqrt(5)) # 元に戻らず面をなす。結び目ができる
par(mfcol=c())
for(c in 1:length(C12s[,1])){
	#C1とC2の最小公倍数*2piだけ時間が経過すると、「元に戻る」
C1<-C12s[c,1]
C2<-C12s[c,2]


x<-y<-z<-rep(0,length(t))
# 初期値は適当
x[1]<-runif(1)
y[1]<-runif(1)
z[1]<-runif(1)

for(i in 2:length(t)){
	x[i]<-x[i-1]+(cos(C1*(t[i-1]+0))+k)*cos(C2*(t[i-1]+0))
	y[i]<-y[i-1]+(cos(C1*(t[i-1]+2/3*pi))+k)*cos(C2*(t[i-1]+2/3*pi))
	z[i]<-z[i-1]+(cos(C1*(t[i-1]+4/3*pi))+k)*cos(C2*(t[i-1]+4/3*pi))
}
open3d()
plot3d(x,y,z,col=rainbow(1000),cex=0.1)
#ylim<-c(min(x,y,z),max(x,y,z))
#plot(x,type="l",col=1,ylim=ylim)
#par(new=TRUE)
#plot(y,type="l",col=2,ylim=ylim)
#par(new=TRUE)
#plot(z,type="l",col=3,ylim=ylim)
}
  • 微分から出発したが
  • \frac{dx}{dt}=(\cos(C1 \times t)+k)*\cos(C2 \times t)は積分できて
  • \frac{dx}{dt}=\cos(C1 \times t)\cos(C2 \times t) +k \cos(C2\times t)
  • [tex:\frac{dx}{dt}=\frac{1}{2}(\cos*1+k\cos(C2 \times t)]から
  • [tex:x=\frac{1}{2} (\frac{1}{C1+C2} \sin*2 + \frac{k}{C2} \sin(C2 t)+Const]
  • Rで書けば
# 周期は0-(2*pi)をX回繰り返しとする
X<-4
# １周あたりの点の数をnpt
npt<-1000
t<-seq(from=0,to=1,length.out=npt)*2*pi*X
# x,y,zの座標は、直径の周期C1と角の周期C2とで決める
# 2piの間にC1回、C2回繰り返しさせる
#C1とC2の最小公倍数*2piだけ時間が経過すると、「元に戻る」
# C1=C2=0のときは直線
# C1=0,C2!=0のときは円
# C1=C2>0のときはらせん
# C1*2=C2 のときは三葉結び目
# C1*k1=C2*k2 ;k1,k2は異なる自然数 のときは複雑な結び目
# C1とC2に最小公倍数がないと、周回して面をなす

C1<-1
C2<-2.4
# x,y,zは相互に対等なので、0/3*2pi, 1/3*2pi,2/3*2pi のずれを入れる
# 直径は周期的で (k-1)と(k+1)の間を上下するものとする
# k=1では三葉結び目は結ばれずに１個の集中点をもち
# k<1では結ばれない
# k>1で、結ばれる
k<-2

x<-1/(2*(C1+C2))*sin((C1+C2)*t)+1/(2*(C1-C2))*sin((C1-C2)*t)+k/C2*sin(C2*t)
y<-1/(2*(C1+C2))*sin((C1+C2)*(t+2/3*pi))+1/(2*(C1-C2))*sin((C1-C2)*(t+2/3*pi))+k/C2*sin(C2*(t+2/3*pi))
z<-1/(2*(C1+C2))*sin((C1+C2)*(t+4/3*pi))+1/(2*(C1-C2))*sin((C1-C2)*(t+4/3*pi))+k/C2*sin(C2*(t+4/3*pi))

plot3d(x,y,z,col=rainbow(1000))
  • 三角関数は複素数でe^{i\theta}=\cos(\theta) + i \sin(\theta)と表せることを使えば、x,y,zの陽なパラメタ表示も複素係数の指数関数の実部(もしくは虚部)として表せる
myi<-complex(real=0,imaginary=1)
xz<-Im(1/2*1/(C1+C2)*exp(myi*(C1+C2)*t))+Im(1/2*1/(C1-C2)*exp(myi*(C1-C2)*t))+Im(k/C2*exp(myi*(C2)*t))
yz<-Im(1/2*1/(C1+C2)*exp(myi*(C1+C2)*(t+2/3*pi)))+Im(1/2*1/(C1-C2)*exp(myi*(C1-C2)*(t+2/3*pi)))+Im(k/C2*exp(myi*(C2)*(t+2/3*pi)))
zz<-Im(1/2*1/(C1+C2)*exp(myi*(C1+C2)*(t+4/3*pi)))+Im(1/2*1/(C1-C2)*exp(myi*(C1-C2)*(t+4/3*pi)))+Im(k/C2*exp(myi*(C2)*(t+4/3*pi)))
par(mfcol=c(1,3))
plot(1/(2*(C1+C2))*sin((C1+C2)*t),Im(1/2*1/(C1+C2)*exp(myi*(C1+C2)*t)))
plot(1/(2*(C1-C2))*sin((C1-C2)*t),Im(1/2*1/(C1-C2)*exp(myi*(C1-C2)*t)))
plot(k/C2*sin(C2*t),Im(k/C2*exp(myi*(C2)*t)))
par(mfcol=c(1,1))

plot3d(xz,yz,zz,col=rainbow(1000))

*1:C1+C2)t)+\cos(C1-C2)t

*2:C1+C2)t)+\frac{1}{C1-C2}\sin((C1-C2)t