- Dirichlet問題と言うのがある
- 周辺の値が固定されているときに内部で調和する関数はどんなものか、という問題
- これをシミュレーションで解いてみることにする(正方形という簡単な領域での単純な例)
- 以下は、こちらにあるJuliaのソース(Juliaのインストールなどはこちら)と、それを翻訳したRのソース
- どういうシミュレーションか、というと:
- 正方形領域の周囲に固定値を与える
- 正方形の内部の点については、
- そこから酔歩をさせて、初めて行き着いた周辺の値が何かを何度も記録し、その平均を持って、そこの値とする
- どういう意味かというと、固定値を持つ周辺から、正方形内部の点に対して平等に影響が及ぶとする。その平等な影響の及び方が酔歩であるとする
- このシミュレーションがRのベクトル処理とかでやったりするのに向かないのは、「酔歩がいつ辺縁に到着するかがわからないから、それを目途にできないため」
- さて、十分に平均がとれるまで回数を回すとするとN=5000とかになるが、それだとjuliaは終わるがRで終わらないので、N=100でやってみるとすると
- juliaが0.4秒くらい
- Rが18秒強
function rytest(N)
nx, ny = 25, 25
f = zeros(nx,ny)
f[1,:], f[nx,:], f[:,1], f[:,ny] = 2, 2, -2, -2
C = zeros(nx,ny)
for x = 1:nx
for y = 1:ny
for k = 1:N
xx = x
yy = y
while ((xx > 1) && (xx < nx) && (yy > 1) && (yy < ny))
b = rand(2)
if (b[1] < .5)
if (b[2] < .5)
xx -= 1
else
xx += 1
end
elseif (b[2] < .5)
yy -= 1
else
yy += 1
end
end
C[x,y] += f[xx,yy]
end
end
end
draw(PNG("dirichlet.png", 10cm, 10cm), spy(C/N))
return 0
end
@time rytest(100)
rytest <- function(N){
nx <- ny <- 25
f <- matrix(0,nx,ny)
f[1,] <- f[nx,] <- 2
f[,1] <- f[,ny] <- -2
z <- matrix(0,nx,ny)
for(i in 1:nx){
for(j in 1:ny){
for(k in 1:N){
xx <- i
yy <- j
while((xx>1) & (xx<nx) & (yy>1) & (yy<ny)){
r <- runif(2)
if(r[1]>0.5){
if(r[2]>0.5){
xx <- xx + 1
}else{
xx <- xx -1
}
}else{
if(r[2]>0.5){
yy <- yy + 1
}else{
yy <- yy -1
}
}
}
z[i,j] <- f[xx,yy]
}
}
}
image(z/N)
return(0)
}
system.time(rytest(100))
- 1次元GFFのjuliaシミュレーションとそのR翻訳版
Pkg.add("Gadfly")
Pkg.add("Cairo")
using Compose, Gadfly, Cairo
L = 500
Z = randn(L-1)
lambda = 2*sin(pi*[1:L-1]/(2*L)).^2
evec = sqrt(2)*sin(pi*kron([1:L-1],[1:L-1]')/L)
GFF1 = zeros(L-1)
for x = 1:L-1
for y = 1:L-1
for n = 1:L-1
GFF1[x] += (lambda[n])^(-1/2)*evec[n,x]*evec[n,y]*Z[y]
end
end
end
p = plot(x = [1:L-1], y = GFF1,
Geom.line,
Guide.xlabel(""),
Guide.ylabel(""),
Theme(default_color=color("black")))
draw(PNG("gff1.png", 10cm, 10cm), p)
L <- 500
Z <- rnorm(L-1)
xy <- expand.grid(1:(L-1),1:(L-1))
xy2 <- matrix(xy[,1]*xy[,2],L-1,L-1)
lambda <- 2*sin(pi*(1:(L-1))/(2*L)) ^2
evec <- sqrt(2)*sin(pi*xy[,1]*xy[,2]/L)
evec2 <- matrix(evec,L-1,L-1)
GFF1 <- rep(0,L-1)
for(i in 1:(L-1)){
for(j in 1:(L-1)){
for(k in 1:(L-1)){
GFF1[i] <- GFF1[i] + lambda[k] ^ (-1/2)*evec2[k,i]*evec2[k,j]*Z[j]
}
}
}
plot(1:(L-1),GFF1,type="l")
library(rgl)
plot3d(xy[,1],xy[,2],z)