Skip to content

Instantly share code, notes, and snippets.

Embed
What would you like to do?
#2D-Turing pattern陰解法
#参考:三浦岳、発生の数理、京都大学学術出版会、2015、6.6
#julia v0.6
#using FFTW
using Plots
gr(size=(320,300),titlefont=font("sans-serif",9))
dx=0.04
N=convert(Int,1/dx) #gridSize=25
dt=0.02
dp=0.0002
dq=0.01
f(u,v)=u+dt*(0.6*u-v-u.^3)
g(u,v)=v+dt*(1.5*u-2*v)
kernelU=zeros(N,N)
kernelU[1,1]=1+4*dt*dp/dx/dx
kernelU[2,1] =-dt*dp/dx/dx
kernelU[end,1]=-dt*dp/dx/dx
kernelU[1,end]=-dt*dp/dx/dx
kernelU[1,2] =-dt*dp/dx/dx
kernelV=zeros(N,N)
kernelV[1,1]=1+4*dt*dq/dx/dx
kernelV[2,1] =-dt*dq/dx/dx
kernelV[end,1]=-dt*dq/dx/dx
kernelV[1,end]=-dt*dq/dx/dx
kernelV[1,2] =-dt*dq/dx/dx
ku=1./(fft(kernelU))
kv=1./(fft(kernelV))
function oneStep(u,v)
u=ifft(ku.*fft(f(u,v)))
v=ifft(kv.*fft(g(u,v)))
return u,v
end
uI=rand(N,N)*0.01
vI=rand(N,N)*0.01
u=uI
v=vI
#=
for i in 1:100
u,v=oneStep(u,v)
end
u
=#
anim = @animate for m in 0:3000
t=round(m*dt,3)
heatmap(real(u),
title="t = $t",
aspect_ratio=:equal
)
u,v=oneStep(u,v)
end every 25
@time gif(anim,"tmp/6_6_2D_pattern_inkaihou_320.gif",fps=4)
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment