たくさんの固定点を持つロトカ=ヴォルテラ様な連立偏微分方程式

• 昨日の続き
• 実数解を持つ実係数関数の作り方についてはこちら
• と表すこととし、を実解を持つ多項式とすることとする
• ,と表せたとする
• 保存量は次を満足する
  • 
  • 
• 積分すると
  • 
  • 
• Rで保存量の等値線を描いてみる
• 状態推移させれば、等値線上を巡回する(ただし、離散的に位置をとると、誤差が生じて発散する。ただし、その誤差は、自然現象(生命現象)でも偶然項として起きえて、そのときは、等値線を乗り換えるはずだから、そちらの方が「自然」とも言えるが今はそのことは気にしないことにする)

# 虚数解のベクトルxsから変数の次数別の複素係数を求める関数
mypoly2<-function(xs){
	library(gtools)
	n<-length(xs)
	ret<-rep(0,n+1)
	for(i in 1:n){
		cmb<-combinations(n,length(xs)-i+1)
		print(cmb)
		for(j in 1:length(cmb[,1])){
			tmp<-1
			for(k in cmb[j,]){
				tmp<-tmp*(-xs[k])
			}
			ret[i]<-ret[i]+tmp
		}
	}
	ret[length(ret)]<-1
	ret
}
k<-5 # 次数
# 次数の数だけ複素数解
xs<-complex(real=rnorm(k),imaginary=rnorm(k))
# 一部は実数解
xs[1:3]<-Re(xs[1:3])
# 複素係数
ks<-mypoly2(xs)
# 実係数
ksRe<-Re(ks)
# 任意の実数について計算
x<-rnorm(1)
# 複素係数関数
print(sum(ks*x^(0:length(xs))))
# 実係数関数 実部は一緒
print(sum(ksRe*x^(0:length(xs))))
# 検算
tmp<-1
for(i in 1:length(xs)){
	tmp<-tmp*(x-xs[i])
}
print(tmp)
# 任意の複素数について計算
x<-complex(real=rnorm(1),imaginary=rnorm(1))
# 複素係数関数
print(sum(ks*x^(0:length(xs))))
# 実係数関数 異なる
print(sum(ksRe*x^(0:length(xs))))
# 変数の実部だけを用いても(もちろん)異なる
print(sum(ksRe*Re(x)^(0:length(xs))))
# 検算
tmp<-1
for(i in 1:length(xs)){
	tmp<-tmp*(x-xs[i])
}
print(tmp)
# 複素解ごとに関数の値を計算
for(i in 1:length(xs)){
	x<-xs[i]
# 複素係数関数
print(sum(ks*x^(0:length(xs))))
# 実係数関数
print(sum(ksRe*x^(0:length(xs))))
# 検算
tmp<-1
for(i in 1:length(xs)){
	tmp<-tmp*(x-xs[i])
}
print(tmp)
}

###############
# fx=vx/x^(kx) 
# gy=uy/y^(ky)

#nx,nyはvx,uyの次数
nx<-7
ny<-7

kx<-2
ky<-2

# 複素解
xs<-complex(real=rnorm(nx),imaginary=rnorm(nx))
ys<-complex(real=rnorm(ny),imaginary=rnorm(ny))

# 実解数 (この例では、すべて実数解にしてある)
nrx<-nx # nrx<= nx
nry<-ny # nry<= ny

# 固定点を集めて、後で、保存量計算する領域を限定しやすくする
xs[1:nrx]<-runif(nrx)*0.5+0.25
ys[1:nry]<-runif(nry)*0.5+0.25

# 最大次数の係数を1として係数ベクトルを作ったので、すべての係数を定数倍する
Cx<-3
Cy<-3

# 係数の実数化
fkx<-Re(mypoly2(xs))*Cx
fky<-Re(mypoly2(ys))*Cy

# 多項式を次数とその係数のベクトルで表すこととする
fkx2<-cbind((0:nx)-kx,fkx)
fky2<-cbind((0:ny)-ky,fky)
# 負の次数を含む多項式の値を返す
mypoly3<-function(x,k){
	ret<-0
	if(is.matrix(k)){
		ret<-sum(k[,2]*x^k[,1])
	}else{
		ret<-sum(k*x^(0:(length(k)-1)))
	}
	ret
}
# 負の次数を含む多項式の積分の値を返す
mypoly3integral<-function(x,k){
	ret<-0
	if(is.matrix(k)){
		nonzero<-which(k[,1]!=-1)
		zero<-which(k[,1]==-1)
		ret<-sum(k[nonzero,2]*x^(k[nonzero,1]+1)/(k[nonzero,1]+1))
		ret<-ret+k[zero,2]*log(x)
	}else{
		ret<-sum(k*x^(1:length(k))/(1:length(k)))
	}
	ret
}

# 保存量Vの計算格子
minxr<-min(Re(xs[1:nrx]))
maxxr<-max(Re(xs[1:nrx]))
minyr<-min(Re(ys[1:nry]))
maxyr<-max(Re(ys[1:nry]))

# 実数解の周辺のみを調べる
kizami<-1/200
kizamiNum<-200
X<-seq(from=minxr*0.9,to=maxxr/0.9,length.out=kizamiNum)
Y<-seq(from=minyr*0.9,to=maxyr/0.9,length.out=kizamiNum)
#X<-Y<-seq(from=kizami,to=1,by=kizami)
#X<-Y<-seq(from=0.25,to=0.75,by=kizami)
XY<-expand.grid(X,Y)
V<-rep(0,length(XY[,1]))

for(i in 1:length(V)){
	V[i]<-mypoly3integral(XY[i,1],fkx2)-mypoly3integral(XY[i,2],fky2)
}

V<-V-min(V)+1
image(X,Y,matrix(V,ncol=length(X)),col = topo.colors(50))
# 等高線を入れる
contour(X,Y,matrix(log(V),ncol=length(X)),nlevels=100,add=TRUE)
# 実数解の位置を示すために白線を入れる(白線の交点が実数解)
for(i in 1:nx){
	abline(v=xs[i],col="white")
}
for(i in 1:ny){
	abline(h=ys[i],col="white")
}