2010-12-09

ソリトンをRで描く（箱玉系の数理　第２章　セルオートマトン）

作者: 時弘哲治
出版社/メーカー: 朝倉書店
発売日: 2010/01/01
メディア: 単行本
購入: 1人クリック: 14回
この商品を含むブログ (6件) を見る

振幅に比例する速度で進む波が形を変えない。これをソリトン波という
ソリトン波は複数あって、ぶつかって、分かれていく
KdV(Korteweg-de Vries)方程式は、このソリトン波(複数のソリトン波の合成も含む)を解とする。ソリトン解という
KdV方程式は、複数のソリトン解をもつ
KdV方程式は、ソリトン解以外の解も持ち、周期波解
KdV方程式は $\frac{\partial u(x,t)}{\partial t} + 6u(x,t)\frac{\partial u(x,t)}{\partial x} + \frac{\partial^3 u(x,t)}{\partial x^3}=0$
KdV方程式には、係数として"6"などが出るが、これは、２次微分の２と３次微分の３の積から生じる６を解消することと関係している
複数(n個)のソリトン波の合成ソリトン波解は、次のような形をしている
- 個々のソリトン波は、速度変数pと、位相変数 $\theta$ を持つ
- 個々のソリトン波は $p_ix - p_i^3 t - \theta_i$ の形をしている。これが、ソリトン波の特長である「振幅に比例する速度で進む」ことを表した形(位置xと時刻tとの係数が独立ではない)となっている
- 複数のを重ね合わせて、KdV方程式の解となるための条件が以下の通りである
  - n個の個別の波のすべての組み合わせ(べき集合)に関する足し合わせとなる
  - べき集合は( $S=\{\phi,\{1\},\{2\},...,\{n\},\{1,2\},\{1,3\},....,\{1,2,...,n\}\}$ で、その要素数は $2^n$
  - べき集合のそれぞれに対応して、足し合わせてKdV方程式を満足するような相殺項が入った形をしている
- 式としては次の通り
  - $u(x,t)=2\frac{\partial^2 \log(\tau(x,t))}{\partial x^2}$
  - $\tau(x,t)=\sum_{s \in S} exp(\sum_{i \in s} \nu_i + \sum_{i,j \in s} A_{i,j})$
  - $\nu_i=p_i x - p_i^3 t -\theta_i$
  - $exp(A_{i,j})=\frac{(p_i-p_j)^2}{(p_i+p_j)^2}$
以下のRソース(末尾）は、作成した関数CalcS0123()を用いて、位置xと時刻tについて計算するもの
波は定速で進んでいるが、衝突すると、ちょっとずれる（位相が変化する）様子が見て取れる
保存量
- KdV方程式は微分に関して特徴的である
- 空間に関する偏導関数の多項式を２つとり、それぞれの時間微分と空間微分との和を取ると、それが０になるという特徴
- それのひとつが、「○○保存の法則」のようなもので「保存量」がいろいろあることを示す
- 「ポテンシャル」がいろいろと定義できることでもある
続きはこちら

n<-4 # ソリトン波の数
ps<-runif(n) # 速度と振幅の係数
thetas<-runif(n) # 位相の係数
# わかりやすいグラフになるように、値を調整
#ps<-rep(1,n)
#ps<-c(1,-2)
ps<-runif(n,-1,1)
ps<-c(-0.4,-0.3,0.2,0.5)
thetas<-runif(n,-1,1)*100

# 次空間を定める
xs<-seq(from=-500,to=500,by=5)
ts<-seq(from=0,to=2000,by=5)

CalcS012<-function(ps,thetas,xs,ts){ # 0,1,2次微分
	n<-length(ps) # ソリトンの数
# 組み合わせ相殺の項
	As<-(outer(ps,ps,FUN="-"))^2/(outer(ps,ps,FUN="+"))^2
	As[which(As!=0,arr.ind=TRUE)]<-log(As[which(As!=0,arr.ind=TRUE)])
# べき集合の要素を作る
	mus<-as.matrix(expand.grid(rep(list(c(0,1)),n)))
# tauの0次微分、1次微分、2次微分ごとに計算する
	ret0<-ret1<-ret2<-matrix(0,length(xs),length(ts))
	for(i in 1:length(mus[,1])){
		X<-mus[i,]*outer(ps,xs,FUN="*")
		T<-mus[i,]*outer(ps^3,ts,FUN="*")
		Th<-mus[i,]*thetas
		A<-as.matrix(outer(mus[i,],mus[i,],FUN="*"))
		A<-sum(A*As)/2
# 微分すると、exp(ax + bt +c)の a が出る。そのa=matX
		matX<-matrix(sum(ps*mus[i,]),length(xs),length(ts))
		tmp<-matrix(rep(apply(X,2,sum),length(ts)),length(xs),length(ts))-matrix(rep(apply(T,2,sum),length(xs)),length(xs),length(ts),byrow=TRUE)-sum(Th)+A
		tmpexp<-exp(tmp)
		ret0<-ret0+exp(tmp)
		ret1<-ret1+matX*exp(tmp)
		ret2<-ret2+matX*matX*exp(tmp)
	}
	ret4<-2/ret0*(ret2-ret1^2/ret0)
	return(ret4)
}

outout<-CalcS012(ps,thetas,xs,ts)

outout[which(outout==-Inf,arr.ind=TRUE)]<-min(outout[is.finite(outout)])
outout[which(outout==Inf,arr.ind=TRUE)]<-max(outout[is.finite(outout)])

#tmpoutout<-outout
#tmpoutout<-tmpoutout[which(tmpoutout!=-Inf,arr.ind=TRUE)]
#tmpoutout<-tmpoutout[which(tmpoutout!=Inf,arr.ind=TRUE)]

#maxtmpoutout<-max(tmpoutout)
#mintmpoutout<-min(tmpoutout)

#replaceoutout<-outout
#replaceoutout[which(replaceoutout==-Inf,arr.ind=TRUE)]<-mintmpoutout
#replaceoutout[which(replaceoutout==Inf,arr.ind=TRUE)]<-maxtmpoutout



filled.contour(outout,col=terrain.colors(15))

image(outout,col=terrain.colors(12))


persp(outout,phi=50,theta=30)
#image(outout)
library(rgl)
plot3d(col(outout),row(outout),outout,col="green")