Skip to content

Instantly share code, notes, and snippets.

@SoTaInverSpinel
Created September 19, 2018 11:54
Show Gist options
  • Star 0 You must be signed in to star a gist
  • Fork 0 You must be signed in to fork a gist
  • Save SoTaInverSpinel/b99093067b265c7065e8d536010965b3 to your computer and use it in GitHub Desktop.
Save SoTaInverSpinel/b99093067b265c7065e8d536010965b3 to your computer and use it in GitHub Desktop.
#2D-Turing pattern陰解法+水玉模様
#参考:三浦岳、発生の数理、京都大学学術出版会、2015、6.8
#julia v0.6
#using FFTW
using Plots
gr()
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.^2-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
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_inkaihou2.gif",fps=4)
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment