経時変化

• マスを仕切って計算していく手法を考えたが（こちらやこちら）
• こちらのような誤差が気になるので、マスを区切らずにたくさんの点の集合を扱うことを考える
• 微分形式に戻って
• それぞれの点が以下の式に従って点の変位を計算していく
  • 
• ソースでは次の式に従っていて一点に収束する
  • 
• この書き方では誤差が数値計算の桁をあげることでおさえられる
  • マスで計算する場合はマスの細かさにあたるので、あまり細かくするとマスの数が増えることになる
  • 個々の振る舞いはこちらのほうが扱いやすそう
  • 一方で確率や分布の計算には工夫が必要となりそう

• 色の付け方に関しては こちら も参照

library(rgl)

#x1<-1
#x2<-0.01
#X<-c(x1,x2)

dt<-0.1
T<-20
Nt<-T/dt
A<-matrix(0,3,Nt+1)
#A[,1]<-c(X,0)
B<-c()

f1<-function(X){
	y<- 1-X[2]
	return(y)
	}
f2<-function(X){
	y<- X[1]^2-X[2]
	return(y)
	}

table<-data.matrix(expand.grid(0:20/10,0:20/10))
dimnames(table)<-NULL
k<-nrow(table)

for(i in 1:k){

X<-table[i,]
A[,1]<-c(X,0)
for(nt in 1:Nt){
dX<-c(f1(X),f2(X))*dt
X<-X+dX
#print(X)
A[,nt+1]<-c(X,nt)
}
B<-cbind(B,A)
}


range<-c(min(B[1:2,]),max(B[1:2,]))
plot3d(B[1,],B[2,],B[3,],xlab="x1",ylab="x2",zlab="time",xlim=range,ylim=range,zlim=c(1,Nt),col=rainbow(Nt+1)[B[3,]+1],type="p")


for(j in 1:(Nt+1)){
	plot(B[1,which(B[3,]==j-1)],B[2,which(B[3,]==j-1)],pch=16,xlab="x1",ylab="x2",xlim=range,ylim=range,col=rainbow(Nt+1)[j])
	par(new=TRUE)
	}