存在から発展へ

シミュレーション 熱力
  • こちら のつづき
  • 反応拡散方程式
    • \begin{matrix} \frac{dX}{dt} &=& A + X^2Y -BX -X + D_X\frac{\par^2X}{\par r^2}\\ \frac{dY}{dt} &=& BX-X^2Y + D_Y\frac{\par^2Y}{\par r^2 } \end{matrix}
  • 下を解として代入する（固定端）
    • \begin{matrix} X &=& A +X_0\sin\frac{n \pi r}{L} \\ Y &=& \frac{B}{A} +Y_0 \sin\frac{n\pi r}{L} \end{matrix}
    • ここで X,Y どちらも同じ n が使われていることに注意しておく
  • 次式が得られる
    • D_x D_x (\frac{n \pi}{L})^4 +\{ A^2D_x+(1-B)D_y \}(\frac{n \pi}{L})^2 +A^2 =0
  • \frac{n}{L} について地道に解いてみる（L = 1 と考える）
    • 図5.5 をみるかぎりでは 8 くらいの値になるはずであるが大きくずれてしまった..
    • rについて任意スケールとなってはいるが、 \frac{r}{L} だから任意といっているはずだと思う..
Dx<-0.0016
Dy<-0.008
A<-2
B<-4.17

# 係数の計算
a<-Dx*Dy*pi^4
b<-(A^2*Dx+(1-B)*Dy)*pi^2
c<-A^2
d<-b^2-4*a*c  #  これは判別式
d>0
n1<-sqrt((-b+sqrt(d))/(2*a))
n2<-sqrt((-b-sqrt(d))/(2*a))
n1;n2

# 実行

> d>0
[1] TRUE
> n1<-sqrt((-b+sqrt(d))/(2*a))
> n2<-sqrt((-b-sqrt(d))/(2*a))
> n1;n2
[1] 11.14744
[1] 5.081013

バッチモードでRを動かす

R 
/* シェルスクリプト */
#!/bin/sh
R --vanilla --slave --args < arg.R > arg.txt  0 1 2 a b c

# Rのスクリプト
x<-commandArgs()
print(x)
  • 実行してできた "arg.txt"の中身
 [1] "/Library/Frameworks/R.framework/Resources/bin/exec/x86_64/R"
 [2] "--vanilla"                                                  
 [3] "--slave"                                                    
 [4] "--args"                                                     
 [5] "0"                                                          
 [6] "1"                                                          
 [7] "2"                                                          
 [8] "a"                                                          
 [9] "b"                                                          
[10] "c"
  • 引数が数字の場合は arg の中から対応する部分を取り出して数字に変換する
# Rのスクリプト
x <- as.numeric(commandArgs()[5:7])
  • --arg 以下（？） 与えた引数のみを取り出すときは trailingOnly = TRUE とすればよい
x<-commandArgs(trailingOnly = TRUE)

# 実行
[1] "0" "1" "2" "a" "b" "c"

分岐

シミュレーション

KABIRA2011-07-30

  • 存在から発展へ
    • p120
    • パラメータBの値をかえたときの挙動
  • 定常状態になるまで時間がかかる
    • T = 10 ほど計算させる必要がある
  • 領域の大きさの影響もないとは言い切れない
    • 以下のソースは図 5.10
      • 1.18 \times 1.18 の領域内で計算するもの
      • 本では R = 0.5861 での計算
      • （ 1.18 \sim 2R）
    • 図 5.11 のような回転解がみれていない
  • マスを増やして計算していると次のようなエラーがでた
    • 解読中...

R(1270,0xa020a4e0) malloc:

mmap(size=2441216) failed (error code=12)
error: can't allocate region
set a breakpoint in malloc_error_break to debug
library(rgl)

A<-2
B<-5.4
X0<-A
Y0<-B/A

Lx<-1.18
dx<-0.02
Nx<-floor(Lx/dx)

Ly<-1.18
dy<-0.02
Ny<-floor(Ly/dy)

T<-1
dt<-0.001
Nt<-T/dt

# 拡散係数
Dx<-0.008
Dy<-0.004

# 補正
Dx<-Dx/dx^2*dt
Dy<-Dy/dy^2*dt


# 円を作る
K<-K2<-matrix(0,Nx,Ny)
for(nx in 1:Nx){
for(ny in 1:Ny){
if((nx-(Nx+1)/2)^2/(Nx/2)^2+(ny-(Ny+1)/2)^2/(Ny/2)^2 < 1){K[nx,ny]<-1}
}}

# 内部は1、境界は2
K2[-1,][which(K[-Nx,]-K[-1,] == -1)]<-2
K2[-Nx,][which(K[-Nx,]-K[-1,] == 1)]<-2
K2[,-1][which(K[,-Ny]-K[,-1]  == -1)]<-2
K2[,-Ny][which(K[,-Ny]-K[,-1]  == 1)]<-2
K2[1,]<-(K[1,]>0)*2;K2[Nx,]<-(K[Nx,]>0)*2;K2[,1]<-(K[,1]>0)*2;K2[,Ny]<-(K[,Ny]>0)*2
K[which(K2 == 2)]<-2


# 各点のX,Yの濃度
C<-array(c(X0),c(Nx,Ny,2))

# C[floor((Nx+1)/2):floor((Nx)/2+1),floor((Ny+1)/2):floor((Ny)/2+1),2]<-Y0+seq(from=0,to=0.1,length=4)
C[floor((Nx+1)/2):floor((Nx)/2+1),floor((Ny+1)/2):floor((Ny)/2+1),1]<-X0
C[,,2]<-rep(Y0,Nx*Ny)+runif(Nx*Ny)*0.1

C[,,1]<-C[,,1]*(K>0)
C[,,2]<-C[,,2]*(K>0)

C[,,1][which(K==2)]<-X0
C[,,2][which(K==2)]<-Y0

# 拡散係数の計算
DX<-array(0,c(Nx-1,Ny,2))
DY<-array(0,c(Nx,Ny-1,2))
DX[,,1]<-(K[-1,]*K[-Nx,]>0)*Dx
DX[,,2]<-(K[-1,]*K[-Nx,]>0)*Dy
DY[,,1]<-(K[,-1]*K[,-Ny]>0)*Dx
DY[,,2]<-(K[,-1]*K[,-Ny]>0)*Dy


# ここからがシミュレーション	
for(t in 1:Nt)
{		
	if(t %% 100 ==1){
	M<-max(C[,,1])
	m<-min(C[,,1][which(C[,,1]!=0)])
		if(M == m){
		Col<-1
		}else{
		Col<-(C[,,1]-m)/(M-m)*100*(K>0)+1
		}
	plot3d(slice.index(C[,,1],1)*Lx/Nx,slice.index(C[,,1],2)*Ly/Ny,C[,,1],alpha=ifelse(K==0,0,1),xlab="x",ylab="y",zlab="X conc.",col=rainbow(101)[Col],main=paste("time = ",(t-1)*dt,sep=""))
#	rgl.snapshot(paste(1000+(t-1)))
	}	
	
	# ブリュッセレータの計算
	tempC<-(K>0)*(A+C[,,1]^2*C[,,2]-B*C[,,1]-C[,,1])*dt
	C[,,2]<-C[,,2]+(K>0)*(B*C[,,1]-C[,,1]^2*C[,,2])*dt
	C[,,1]<-C[,,1]+tempC
	
	
	
	# 拡散の計算
	dC_x<-C[-1,,]-C[-Nx,,]
	dC_y<-C[,-1,]-C[,-Ny,]

	Ix<- -DX*dC_x
	Iy<- -DY*dC_y

	C[-1,,]<-C[-1,,]+Ix
	C[-Nx,,]<-C[-Nx,,]-Ix
	C[,-1,]<- C[,-1,]+Iy
	C[,-Ny,]<-C[,-Ny,]-Iy
	
	# 固定端
#	C[,,1][which(K==2)]<-X0
#	C[,,2][which(K==2)]<-Y0
}

#色で表示するだけにする
library(rgl)

A<-2
B<-5.4
X0<-A
Y0<-B/A

Lx<-1.18
dx<-0.02
Nx<-floor(Lx/dx)

Ly<-1.18
dy<-0.02
Ny<-floor(Ly/dy)

T<-1
dt<-0.001
Nt<-T/dt

# 拡散係数
Dx<-0.008
Dy<-0.004

# 補正
Dx<-Dx/dx^2*dt
Dy<-Dy/dy^2*dt


# 円を作る
K<-K2<-matrix(0,Nx,Ny)
for(nx in 1:Nx){
for(ny in 1:Ny){
if((nx-(Nx+1)/2)^2/(Nx/2)^2+(ny-(Ny+1)/2)^2/(Ny/2)^2 < 1){K[nx,ny]<-1}
}}

# 内部は1、境界は2
K2[-1,][which(K[-Nx,]-K[-1,] == -1)]<-2
K2[-Nx,][which(K[-Nx,]-K[-1,] == 1)]<-2
K2[,-1][which(K[,-Ny]-K[,-1]  == -1)]<-2
K2[,-Ny][which(K[,-Ny]-K[,-1]  == 1)]<-2
K2[1,]<-(K[1,]>0)*2;K2[Nx,]<-(K[Nx,]>0)*2;K2[,1]<-(K[,1]>0)*2;K2[,Ny]<-(K[,Ny]>0)*2
K[which(K2 == 2)]<-2


# 各点のX,Yの濃度
C<-array(c(X0),c(Nx,Ny,2))

# C[floor((Nx+1)/2):floor((Nx)/2+1),floor((Ny+1)/2):floor((Ny)/2+1),2]<-Y0+seq(from=0,to=0.1,length=4)
C[floor((Nx+1)/2):floor((Nx)/2+1),floor((Ny+1)/2):floor((Ny)/2+1),1]<-X0
C[,,2]<-rep(Y0,Nx*Ny)+runif(Nx*Ny)*0.1

C[,,1]<-C[,,1]*(K>0)
C[,,2]<-C[,,2]*(K>0)

C[,,1][which(K==2)]<-X0
C[,,2][which(K==2)]<-Y0

# 拡散係数の計算
DX<-array(0,c(Nx-1,Ny,2))
DY<-array(0,c(Nx,Ny-1,2))
DX[,,1]<-(K[-1,]*K[-Nx,]>0)*Dx
DX[,,2]<-(K[-1,]*K[-Nx,]>0)*Dy
DY[,,1]<-(K[,-1]*K[,-Ny]>0)*Dx
DY[,,2]<-(K[,-1]*K[,-Ny]>0)*Dy


# ここからがシミュレーション	
for(t in 1:Nt)
{		
	if(t %% 100 ==1){
	M<-max(C[,,1])
	m<-min(C[,,1][which(C[,,1]!=0)])
		if(M == m){
		Col<-1
		}else{
		Col<-(C[,,1]-m)/(M-m)*100*(K>0)+1
		}
	plot3d(slice.index(K,1)*Lx/Nx,slice.index(K,2)*Ly/Ny,0,alpha=ifelse(K==0,0,1),xlab="x",ylab="y",zlab="X conc.",col=rainbow(101)[Col],main=paste("time = ",(t-1)*dt,sep=""))
	rgl.viewpoint(theta=0,phi=0)
#	rgl.snapshot(paste(1000+(t-1)))
	}	
	
	# ブリュッセレータの計算
	tempC<-(K>0)*(A+C[,,1]^2*C[,,2]-B*C[,,1]-C[,,1])*dt
	C[,,2]<-C[,,2]+(K>0)*(B*C[,,1]-C[,,1]^2*C[,,2])*dt
	C[,,1]<-C[,,1]+tempC
	
	
	
	# 拡散の計算
	dC_x<-C[-1,,]-C[-Nx,,]
	dC_y<-C[,-1,]-C[,-Ny,]

	Ix<- -DX*dC_x
	Iy<- -DY*dC_y

	C[-1,,]<-C[-1,,]+Ix
	C[-Nx,,]<-C[-Nx,,]-Ix
	C[,-1,]<- C[,-1,]+Iy
	C[,-Ny,]<-C[,-Ny,]-Iy
	
	# 固定端
#	C[,,1][which(K==2)]<-X0
#	C[,,2][which(K==2)]<-Y0
}

時間を閉じる

基礎 熱力
  • こちら のコメントから
  • 時間的に閉じるとはどういうかとか
  • 空間の閉じた・開いたについては
    • 空間内外の間で量のやり取りのないとき 閉じた系
    • 空間内外の間で量のやり取りがあれば 開いた系
  • 時間的に閉じるとしたら？
    • ある時刻と別の時刻との間で量の移動がないということ？
    • 時間を切り出したり平行移動できない系が閉じていない系？
    • だろうか
  • 集合の閉じているかいないか
    • ある元を考え、この元の存在している空間を用意しておく
    • さらにこの空間（のベキ集合）からこの空間へのある写像が定義されている
    • この空間内のある部分集合について、この集合（のベキ）の写像による像がふたたびこの集合に包まれるならば、この写像に関してこの集合は閉じている
  • 時間発展においては軌跡が閉じているかどうか
    • ある変数が存在し、この変数の時間に関する微分方程式がある
    • ある初期値から出発し、時間発展を計算しながら空間内の部分集合を生成するとする
    • この集合がとじているかどうか
    • カオスならばこのようにしてつくられる集合は閉じない？
    • 例としてあがっているトーラス上に限定されている場合、それは閉じているのだろうか
    • おそらくエルゴード的
      • 集合自体は閉じていない
      • 軌跡がトーラス上にあることから"トーラスにのる"という性質がすべての時刻で成り立っている
      • これは何かしらの保存量にあたるはず
  • この意味で時間的に閉じるかどうかというのなら、何かしらの時間を生成する方法が与えられれば可能かも

分岐：ブリュッセレータ

シミュレーション R

KABIRA2011-07-27

  • こちら から
    • p120の図5.7
  • 中心部にちいさな揺らぎを与えてその後の時間発展をみる
    • 空間の大きさに影響を受けるようで、半径 R = 0.2 の円の中だけでシミュレーションしてみたが、円柱対称な波が非常に小さく、平面の様になってしまう
    • R = 0.4 ほどの大きさで計算しておいて R \leq 0.2 の範囲を見てみるとたしかに 図5.7 の様な定常状態になるようである
library(rgl)

A<-2
B<-4.6
X0<-A
Y0<-B/A

Lx<-0.82
dx<-0.02
Nx<-floor(Lx/dx)

Ly<-0.82
dy<-0.02
Ny<-floor(Ly/dy)

T<-30
dt<-0.001
Nt<-T/dt

# 拡散係数
Dx<-0.0016
Dy<-0.005

# 補正
Dx<-Dx/dx^2*dt
Dy<-Dy/dy^2*dt


# 円を作る
K<-K2<-matrix(0,Nx,Ny)
for(nx in 1:Nx){
for(ny in 1:Ny){
if((nx-(Nx+1)/2)^2/(Nx/2)^2+(ny-(Ny+1)/2)^2/(Ny/2)^2 < 0.9){K[nx,ny]<-1}
}}

# 内部は1、境界は2
K2[-1,][which(K[-Nx,]-K[-1,] == -1)]<-2
K2[-Nx,][which(K[-Nx,]-K[-1,] == 1)]<-2
K2[,-1][which(K[,-Ny]-K[,-1]  == -1)]<-2
K2[,-Ny][which(K[,-Ny]-K[,-1]  == 1)]<-2
K2[1,]<-(K[1,]>0)*2;K2[Nx,]<-(K[Nx,]>0)*2;K2[,1]<-(K[,1]>0)*2;K2[,Ny]<-(K[,Ny]>0)*2
K[which(K2 == 2)]<-2


# 各点のX,Yの濃度
C<-array(c(X0),c(Nx,Ny,2))

C[floor((Nx+1)/2):floor((Nx)/2+1),floor((Ny+1)/2):floor((Ny)/2+1),2]<-Y0+seq(from=0,to=0.1,length=4)
C[floor((Nx+1)/2):floor((Nx)/2+1),floor((Ny+1)/2):floor((Ny)/2+1),1]<-X0

C[,,1]<-C[,,1]*(K>0)
C[,,2]<-C[,,2]*(K>0)

C[,,1][which(K==2)]<-X0
C[,,2][which(K==2)]<-Y0

# 拡散係数の計算
DX<-array(0,c(Nx-1,Ny,2))
DY<-array(0,c(Nx,Ny-1,2))
DX[,,1]<-(K[-1,]*K[-Nx,]>0)*Dx
DX[,,2]<-(K[-1,]*K[-Nx,]>0)*Dy
DY[,,1]<-(K[,-1]*K[,-Ny]>0)*Dx
DY[,,2]<-(K[,-1]*K[,-Ny]>0)*Dy

# ブリュッセレータを計算する関数
f1<-function(X){
	y<- A+X[1]^2*X[2]-B*X[1]-X[1]
	return(y)
	}
f2<-function(X){
	y<- B*X[1]-X[1]^2*X[2]
	return(y)
	}

# ここからがシミュレーション	
for(t in 1:Nt)
{		
	if(t %% 100 ==1){
	plot3d(slice.index(C[,,1],1)*Lx/Nx,slice.index(C[,,1],2)*Ly/Ny,C[,,1],alpha=ifelse(K==0,0,1),xlab="x",ylab="y",zlab="X conc.",col=rainbow(101)[C[,,1]/max(C[,,1])*100+1],main=paste("time = ",(t-1)*dt,sep=""))
	}	
	
	# ブリュッセレータの計算
	for(nx in 1:Nx){
	for(ny in 1:Ny)
	if(K[nx,ny] >= 1){
	X<-C[nx,ny,]
	dX<-c(f1(X),f2(X))*dt
	C[nx,ny,]<-X+dX
	}
	}
	
	# 拡散の計算
	dC_x<-C[-1,,]-C[-Nx,,]
	dC_y<-C[,-1,]-C[,-Ny,]

	Ix<- -DX*dC_x
	Iy<- -DY*dC_y

	C[-1,,]<-C[-1,,]+Ix
	C[-Nx,,]<-C[-Nx,,]-Ix
	C[,-1,]<- C[,-1,]+Iy
	C[,-Ny,]<-C[,-Ny,]-Iy
	
	# 固定端
#	C[,,1][which(K==2)]<-X0
#	C[,,2][which(K==2)]<-Y0
}

分岐：ブリュッセレータ

R シミュレーション

KABIRA2011-07-26

  • 存在から発展へ
  • 拡散係数を補正する処理を入れる
    • p121 図5.6
    • 拡散係数 D_x, \, D_y とシミュレーションする空間 K の境界の計算も修正
  • 中心部に小さな揺らぎをあたえる
  • はじめ全体が振動するがやがて対称性の破れた定常状態となる
    • Nx, \, Ny がともに奇数だと完璧に対称になってしまうので、対称性の崩れはみえない
library(rgl)

A<-2
B<-4.6
X0<-A
Y0<-B/A

Lx<-0.2
dx<-0.01
Nx<-floor(Lx/dx)

Ly<-0.2
dy<-0.01
Ny<-floor(Ly/dy)

T<-30
dt<-0.001
Nt<-T/dt

# 拡散係数
Dx<-0.00325
Dy<-0.0162

# 補正
Dx<-Dx/dx^2*dt
Dy<-Dy/dy^2*dt


# 円を作る
K<-K2<-matrix(0,Nx,Ny)
for(nx in 1:Nx){
for(ny in 1:Ny){
if((nx-(Nx+1)/2)^2/(Nx/2)^2+(ny-(Ny+1)/2)^2/(Ny/2)^2 < 1){K[nx,ny]<-1}
}}

# 内部は1、境界は2
K2[-1,][which(K[-Nx,]-K[-1,] == -1)]<-2
K2[-Nx,][which(K[-Nx,]-K[-1,] == 1)]<-2
K2[,-1][which(K[,-Ny]-K[,-1]  == -1)]<-2
K2[,-Ny][which(K[,-Ny]-K[,-1]  == 1)]<-2
K2[1,]<-(K[1,]>0)*2;K2[Nx,]<-(K[Nx,]>0)*2;K2[,1]<-(K[,1]>0)*2;K2[,Ny]<-(K[,Ny]>0)*2
K[which(K2 == 2)]<-2


# 各点のX,Yの濃度
C<-array(c(X0),c(Nx,Ny,2))

C[floor((Nx+1)/2):floor((Nx)/2+1),floor((Ny+1)/2):floor((Ny)/2+1),2]<-Y0+seq(from=0,to=0.1,length=4)
C[floor((Nx+1)/2):floor((Nx)/2+1),floor((Ny+1)/2):floor((Ny)/2+1),1]<-X0

C[,,1]<-C[,,1]*(K>0)
C[,,2]<-C[,,2]*(K>0)

C[,,1][which(K==2)]<-X0
C[,,2][which(K==2)]<-Y0

# 拡散係数の計算
DX<-array(0,c(Nx-1,Ny,2))
DY<-array(0,c(Nx,Ny-1,2))
DX[,,1]<-(K[-1,]*K[-Nx,]>0)*Dx
DX[,,2]<-(K[-1,]*K[-Nx,]>0)*Dy
DY[,,1]<-(K[,-1]*K[,-Ny]>0)*Dx
DY[,,2]<-(K[,-1]*K[,-Ny]>0)*Dy

# ブリュッセレータを計算する関数
f1<-function(X){
	y<- A+X[1]^2*X[2]-B*X[1]-X[1]
	return(y)
	}
f2<-function(X){
	y<- B*X[1]-X[1]^2*X[2]
	return(y)
	}

# ここからがシミュレーション	
for(t in 1:Nt)
{		
	if(t %% 1000 ==1){
	plot3d(slice.index(C[,,1],1)*Lx/Nx,slice.index(C[,,1],2)*Ly/Ny,C[,,1],alpha=ifelse(K==0,0,1),xlab="x",ylab="y",zlab="X conc.",col=rainbow(101)[C[,,1]/max(C[,,1])*100+1],main=paste("time = ",(t-1)*dt,sep=""))
	}	
	
	# ブリュッセレータの計算
	for(nx in 1:Nx){
	for(ny in 1:Ny)
	if(K[nx,ny] >= 1){
	X<-C[nx,ny,]
	dX<-c(f1(X),f2(X))*dt
	C[nx,ny,]<-X+dX
	}
	}
	
	# 拡散の計算
	dC_x<-C[-1,,]-C[-Nx,,]
	dC_y<-C[,-1,]-C[,-Ny,]

	Ix<- -DX*dC_x
	Iy<- -DY*dC_y

	C[-1,,]<-C[-1,,]+Ix
	C[-Nx,,]<-C[-Nx,,]-Ix
	C[,-1,]<- C[,-1,]+Iy
	C[,-Ny,]<-C[,-Ny,]-Iy
	
	# 固定端
#	C[,,1][which(K==2)]<-X0
#	C[,,2][which(K==2)]<-Y0
}