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

試行錯誤

R 常微分方程式
  • こちらで常微分方程式を０からやってみている
  • 常微分方程式があって、変数t(時間のようなもの)があるときにy(t)という関数があって、そのr階の導関数をy^{r}(t)としたときに、y^n(t)=F(y(t),y^1(t),...,y^{n-1}(t))と表されるとしよう
  • この常微分方程式の関数Fにはtが入っていない。これは、t(時間)については、「平等な世界」を仮定していることに相当する
  • y^n(t)=F(y(t),y^1(t),...,y^{n-1}(t))はたった一つの方程式であるが、y^{r+1}(t)=\frac{d y^{r}(t)}{dt}という関係式がn-1個あって、定義から前提としていることを忘れない
  • さて、今、y^n(t)はy^r(t);r=0,1,2,...,n-1の関数であるから、y^{r}(t)にそれぞれ軸を与えた\mathbf{R}^{n}空間を考えて、その空間上の点(y_0,y_1,...,y_{n-1})にy^n(t)の値を与えるスカラー関数f(\mathbf{y})を考えよう
  • このn次元空間は(y_1,y_2,...,y_{n-},f(\mathbf{y}=(y_0,y_1,...,y_{n-1}))というベクトルが与えられたベクトル場と考えることもできる
  • n次元のベクトル場の空間なので、任意の点から軌道が生じる
  • この軌道が発散しないようにしたい
  • 一つの方法として、原点から点の位置への位置ベクトル(y_0,y_1,...,y_{n-1})とこの点の運動ベクトルとが必ず直交するようなベクトル場であれば、原点に落ち込むこともない代わりに、発散もしないだろう
  • そのような条件は(y_0,y_1,...,y_{n-1}) * (y_1,y_2,...,y_{n-1},,f(\mathbf{y}))=0
  • これはf(\mathbf{y})=\frac{1}{y_{n-1}} \sum_{i=0}^{n-2} y_i \times y_{i+1};y_{n-1}\ne 0だろう
  • Rで書けば以下のようになるが、これって、なんだろう？？？
  • この話は、こちらとこちらに続く
n<-3

Niter<-10000
dx<-1/1000

Ys<-matrix(0,Niter,n+1)

loop<-TRUE
while(loop){
	Ys[1,]<-c(rnorm(n),0)*0.001
	if(sum(Ys[1,1:(n-2)]*Ys[1,2:(n-1)])<0)loop<-FALSE
}
#Ys[1,]<-c(-1,1,0.5,0)
shogen<-sign(Ys[1,])
TmpFx<-function(ys,Equal=FALSE,a=0.1){
	tmp<-0
	ret<-0
	for(i in 1:(length(ys)-1)){
		tmp<-tmp+ys[i]*ys[i+1]
	}
	ret<--tmp/ys[length(ys)]
	if(!Equal){
		ret<-ret+sign(ys[length(ys)])*a
	}
	ret
}

Ys[1,n+1]<-TmpFx(Ys[1,1:n])
alpha<-0.01
for(i in 2:Niter){
	V<-c(Ys[i-1,2:(n+1)],0)
	Vnorm<-sqrt(sum(V^2))
	Ys[i,]<-Ys[i-1,]+dx/max(c(Vnorm))*V
	Ys[i,n+1]<-TmpFx(Ys[i,1:n],Equal=TRUE)
	zeroq<-sum(Ys[i-1,2:(n+1)]*Ys[i-1,1:n])
	if(abs(zeroq)>0.0001)print(zeroq)
	tmpshogen<-sign(Ys[i,])
	#if(sum((shogen-tmpshogen)^2)!=0){
	#	print("chage shogen")
	#	tmpshogen2<-abs(shogen-tmpshogen)*(-0.5)+1
	#	tmpshogen3<-(tmpshogen2-1)*(-1)
	#	Ys[i,]<-Ys[i,]*tmpshogen2+c(rep(alpha,n)*tmpshogen3[1:n],0)
	#}
	
}
matplot(Ys[,1:n],type="l")

plot(as.data.frame(Ys))

library(rgl)
plot3d(Ys[,1:3],col=rainbow(1000),xlim=range(Ys[,1:3]),ylim=range(Ys[,1:3]),zlim=range(Ys[,1:3]))