試行錯誤 - ryamadaのコンピュータ・数学メモ

こちらで常微分方程式を０からやってみている
常微分方程式があって、変数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]))