In [58]:
using PyPlot;
PyPlot.svg(true);
In [59]:
function csd_jac(f,g,u,v,h=im*1e-12)
Dkᵤ = zeros(Float64,(2,2))
Dkᵤ[1,1] = real(f(u+h,v)/h)
Dkᵤ[1,2] = real(f(u,v+h)/h)
Dkᵤ[2,1] = real(g(u+h,v)/h)
Dkᵤ[1,2] = real(g(u,v+h)/h)
return Dkᵤ
end
Out[59]:
In [60]:
# This function defines the cellular dynamics and the forward derivatives
function fitzhugh_nagumo(u,ϵ=0.01,α=1.5,β=0.1,γ=0.5,t=1)
f = (u,v)->(u - v - u.^3)
g = (u,v)->(ϵ*(sign(u).*abs(u).^(1.0+γ) - α*v + β))
u,u′= reim(u[:])
kᵤ = zeros(u) + zeros(u′)
kᵤ[1] = f(u[1],u[2])
kᵤ[2] = g(u[1],u[2])
Dkᵤ = csd_jac(f,g,u[1],u[2])
if t<0
Dkᵤ = transpose(Dkᵤ)
end
kᵤ += im*Dkᵤ*u′
return kᵤ
end
Out[60]:
In [61]:
# This function defines the cellular dynamics and the forward derivatives
function karma(u,ϵ=0.01,u⁺=1.5415,M=4.0,β=1.389,s=1.0,t=1)
Θ = (x)->(1.0 + tanh(s*x))/2.0
f = (u,v)->((u⁺ - v^M)*(1.0 - tanh(u-3.0))*u*u/2.0 - u)
g = (u,v)->(ϵ*(β*Θ(u-1.0) + Θ(v-1.0)*(v-1.0) - v))
u,u′= reim(u[:])
kᵤ = zeros(u) + zeros(u′)
kᵤ[1] = f(u[1],u[2])
kᵤ[2] = g(u[1],u[2])
Dkᵤ = csd_jac(f,g,u[1],u[2])
if t<0
Dkᵤ = transpose(Dkᵤ)
end
kᵤ += im*Dkᵤ*u′
return kᵤ
end
Out[61]:
In [62]:
# This function defines the cellular dynamics and the forward derivatives
function barkley(u,ϵ=0.08,a=0.75,b=0.0006,c=1.0,t=1)
f = (u,v)->((u.*(1.0-u).*(u - (v+b)./a)))
g = (u,v)->(ϵ.*(u.^c-v))
u,u′ = reim(u[:])
kᵤ = zeros(u) + zeros(u′)
kᵤ[1] = f(u[1],u[2])
kᵤ[2] = g(u[1],u[2])
Dkᵤ = csd_jac(f,g,u[1],u[2])
if t<0
Dkᵤ = transpose(Dkᵤ)
end
kᵤ += im*Dkᵤ*u′
return kᵤ
end
Out[62]:
In [63]:
# This function defines the cellular dynamics and the forward derivatives
function aliev_panfilov(u,ϵ=0.1,a=0.01,b=1.0,c=0.05,k=4.0,t=1)
f = (u,v)->(k.*u.*(1.0 - u).*(u - a) - u.*v)
g = (u,v)->((ϵ + (1.0-ϵ).*(tanh(b.*(c-u))-1.0)).*(k*u - v))
u,u′ = reim(u[:])
kᵤ = zeros(u) + zeros(u′)
kᵤ[1] = f(u[1],u[2])
kᵤ[2] = g(u[1],u[2])
Dkᵤ = csd_jac(f,g,u[1],u[2])
if t<0
Dkᵤ = transpose(Dkᵤ)
end
kᵤ += im*Dkᵤ*u′
return kᵤ
end
Out[63]:
In [64]:
function fourier_operators(Ω,L)
ndim = length(Ω)
Q = Array(Array{Int,1},ndim)
q = [0,0,0]
for n=1:ndim
if iseven(Ω[n])
Q[n] = [0:Ω[n]/2-1,0,1-(Ω[n]/2):-1]
elseif isodd(Ω[n])
Q[n] = [0:(Ω[n]-1)/2-1,0,1-(1+Ω[n])/2:-1]
end
end
∇ = Array(Array{Complex{Float64},ndim},ndim)
Δ = zeros(Ω)
q = [0,0,0]
for n=1:ndim
∇[n] = Complex{Float64}[-im*(2π/L[n])*q[n] for q[1]=Q[1], q[2]=Q[2], q[3]=Q[3]]
Δ += ∇[n].*∇[n]
end
return ∇,Δ
end
Out[64]:
In [391]:
function initialize_field(nvar,Ω)
u = Array(Array{Complex{Float64},length(Ω)},nvar)
nrm = Array(Complex{Float64},2)
u[1] = Complex{Float64}[ exp(sin(x)).*sin(x).*cos(z)
for x=(1:Ω[1])*(2π/Ω[1]), y=(1:Ω[2])*(2π/Ω[2]), z=(1:Ω[3])*(2π/Ω[3])]
u[2] = Complex{Float64}[ exp(sin(x-π)).*exp(sin(2*y-π)).*exp(sin(z))
for x=(1:Ω[1])*(2π/Ω[1]), y=(1:Ω[2])*(2π/Ω[2]), z=(1:Ω[3])*(2π/Ω[3])]
for n=1:nvar
nrm[n] = norm(real(u[n][:]),Inf) + im*norm(imag(u[n][:]),Inf)
u[n] = real(u[n])./max(1.,real(nrm[n])) + im.*imag(u[n])./max(1.,imag(nrm[n]))
end
return u
end
Out[391]:
In [392]:
function translate!(u,φ,Ω,∇)
a = zeros(Int32,4)
û = Array(Array{Complex{Float64},3},nvar)
û′= Array(Array{Complex{Float64},3},nvar)
for a[4]=1:2
û[a[4]] = fft(real(u[a[4]]))
û′[a[4]]= fft(imag(u[a[4]]))
for a[1]=1:Ω[1]
for a[2]=1:Ω[2]
for a[3]=1:Ω[3]
Tφ = exp(φ[1].*∇[1][a[1],a[2],a[3]])
Tφ .*= exp(φ[2].*∇[2][a[1],a[2],a[3]])
Tφ .*= exp(φ[3].*∇[3][a[1],a[2],a[3]])
û[a[4]][a[1],a[2],a[3]] = Tφ .* û[a[4]][a[1],a[2],a[3]]
û′[a[4]][a[1],a[2],a[3]] = Tφ .* û′[a[4]][a[1],a[2],a[3]]
end
end
end
u[a[4]] = real(ifft((û[a[4]]))) + im*real(ifft((û′[a[4]])))
end
return u
end
Out[392]:
In [393]:
function shift_angles(u,Ω,L,∇)
φ = zeros(Float64,length(Ω))
a = zeros(Int32,4)
nvar = length(u)
û = Array(Array{Complex{Float64},3},nvar)
û′= Array(Array{Complex{Float64},3},nvar)
for a[4]=1:2
û[a[4]] = fft(real(u[a[4]]))
û′[a[4]]= fft(imag(u[a[4]]))
end
if Ω[1]>1
if abs(û[1][2,1,1]) > Ω[1]*eps(abs(û[1][2,1,1]))
φ[1] = angle(û[1][2,1,1])
end
end
if Ω[2]>1
if abs(û[1][1,2,1]) > Ω[1]*eps(abs(û[1][1,2,1]))
φ[2] = angle(û[1][1,2,1])
end
end
if Ω[3]>1
if abs(û[1][1,1,2]) > Ω[1]*eps(abs(û[1][1,1,2]))
φ[3] = angle(û[1][1,1,2])
end
end
return φ
end
Out[393]:
In [394]:
function spectral_slice!(u,Ω,L,∇)
φ = shift_angles(u,Ω,L,∇)
# rescale angles to L-space
x = φ;
for n=1:3
x[n] = φ[n].*(L[n]/(2π));
end
# Now rigidly shift the solution by x = -c*hₜ,
# so the phase of the first Fourier modes vanish
translate!(u,x,Ω,∇)
# c should be normalized by the time-step h, however,
# so that I don't need to pass in h to this function,
# we'll handle this detail in the integration code.
return u,x
end
Out[394]:
In [395]:
function rkn(u0,f,h,A=[[0//1 1//2 1//2 1//1],[1//6 1//3 1//3 1//6]])
k = zeros(size(u0))
d = zeros(size(u0))
N = size(A,2)
M = 1
n = 1
k = f(u0+A[1,n]*h*k)
d += h*A[2,n]*k
δ = norm(real(k[:]),Inf)
if δ > 0.1
M = min(10, δ/0.1)
end
h = h/M
for n=2:N
k = f(u0+A[1,n]*h*k)
d += h*A[2,n]*k
end
for m=2:M
for n=1:N
k = f(u0+A[1,n]*h*k)
d += h*A[2,n]*k;
end
end
return d
end
Out[395]:
In [396]:
function reaction_diffustion_etd2rk2(u, f, δ, Ω, L, T=10.0, h=0.1)
# Solves:
# u_t = u_xx + f(u,v), v_t = δ v_xx + g(u,v)
∇,Δ = fourier_operators(Ω,L)
A = [[0//1 2//3],[1//4 3//4]]
fast_fft = plan_fft!(u[1])
fast_ifft = plan_ifft!(u[1])
N = int(ceil(T/h))
h = T/N
c = zeros(Float64,(N,length(Ω)));
#u,c[1,:] = spectral_slice!(u,Ω,L,∇)
a = zeros(Int32,4);
for n=1:N
# advance by h using L
for a[4]=1:2
u[a[4]] = fast_fft(u[a[4]])
u[a[4]] = exp(h * δ[a[4]] * Δ).*u[a[4]]
u[a[4]] = fast_ifft(u[a[4]])
end
for a[1]=1:Ω[1]
for a[2]=1:Ω[2]
for a[3]=1:Ω[3]
u0 = [u[1][a[1],a[2],a[3]],u[2][a[1],a[2],a[3]]];
u0 += rkn(u0,f,h,A);
for a[4]=1:2
u[a[4]][a[1],a[2],a[3]] = u0[a[4]]
end
end
end
end
u,c[n,:] = spectral_slice!(u,Ω,L,∇)
c[n,:] = c[n,:]/h
if isnan(norm(u[1][:]))
println("Divergence!")
break
end
end
return u,c
end
Out[396]:
In [397]:
type stateVector
end
type field
Ω
L
u
end
function field2vector(u::field)
end
function vector2field(x::stateVector)
end
function show(x::stateVector)
try
PyPlot.plot3D(x.u[1][1],x.u[1][2],x.u[1][3],".k");
end
end
function show(x::field)
if x.Ω[1]==1
try
PyPlot.plot(x.u[1],x.u[2],".k");
end
elseif x.Ω[2]==1
try
PyPlot.plot((1:x.Ω[1])*(x.L[1]/x.Ω[1]), x.u[1][:,1,1], "-k");
hold(true);
PyPLot.plot((1:x.Ω[1])*(x.L[1]/x.Ω[1]), x.u[2][:,1,1], "-g");
end
elseif x.Ω[3]==1
try
PyPlot.imshow(real(u[1][:,:,1]),cmap=PyPlot.ColorMap("Spectral"));
end
else
try
for n=1:x.Ω[3]
PyPlot.plot3D((n-1)+real(u[1][:,:,n]));
hold(true);
end
end
end
end
Out[397]:
In [398]:
Ω = (96,48,12)
L₀ = 192.0
L = (L₀*Ω[1]/maximum(Ω),L₀*Ω[2]/maximum(Ω),L₀*Ω[3]/maximum(Ω))
x = (1:Ω[1])*L[1]/Ω[1]
nvar = 2
u = initialize_field(nvar, Ω);
PyPlot.plot(x,real(u[1][:,1,1]),"-k",x,real(u[2][:,1,1]),"-g");
∇,Δ = fourier_operators(Ω,L);
u,c = spectral_slice!(u,Ω,L,∇);
PyPlot.plot(x,real(u[1][:,1,1]),"--k",x,real(u[2][:,1,1]),"--g");
In [399]:
δ=[9//4,3//2 * 9//4]
qq = int(16)
T = 2.0
h = 0.25
C = zeros(int(ceil(T/h))*qq,3);
PyPlot.plot(x,real(u[1][:,1,1]),"-k",x,real(u[2][:,1,1]),"-g");
for n=1:qq
println(n," of ",qq)
u,c = reaction_diffustion_etd2rk2(u,fitzhugh_nagumo,δ,Ω,L,T,h)
C[1+ceil(T/h)*(n-1):ceil(T/h)*n,:] = c
PyPlot.plot(x,2*(n)+real(u[1][:,1,1]),"-k",x,2*(n)+real(u[2][:,1,1]),"-g")
end
In [400]:
PyPlot.plot(h:h:qq*T, cumsum(C[:,1]*h/L[1]), linestyle="-");
PyPlot.plot(h:h:qq*T, cumsum(C[:,2]*h/L[2]), linestyle="-");
PyPlot.plot(h:h:qq*T, cumsum(C[:,3]*h/L[3]), linestyle="-");
In [401]:
PyPlot.plot(h:h:qq*T, (C[:,1]*h/L[1]), linestyle="-");
PyPlot.plot(h:h:qq*T, (C[:,2]*h/L[2]), linestyle="-");
PyPlot.plot(h:h:qq*T, (C[:,3]*h/L[3]), linestyle="-");
In [413]:
for n=[Ω[3],round(Ω[3]/2),1]
PyPlot.surf((n-1)/Ω[3] + real(u[1][:,:,n]),rstride=2,cstride=2,cmap=PyPlot.ColorMap("Spectral"))
end
In [403]:
PyPlot.plot3D(cumsum(C[:,1]*h),cumsum(C[:,2]*h),cumsum(C[:,3]*h),"-k");
In [404]:
PyPlot.plot3D(C[:,1]*h,C[:,2]*h,C[:,3]*h,"-k");
In [405]:
for n=[1,round(Ω[2]/2),Ω[2]]
PyPlot.plot(x,(n-1)/Ω[2] + real(u[1][:,n,1]),"-k",x,real(u[2][:,n,1]),"-g");
end
In [406]:
PyPlot.imshow(real(u[1][:,:,1]));
In [407]:
function vorticity(u,v,Ω,∇)
û = fft(real(u))
v̂ = fft(real(v))
∇u = Array(Array{Complex{Float64},length(Ω)},length(Ω))
∇v = Array(Array{Complex{Float64},length(Ω)},length(Ω))
for n=1:length(Ω)
∇u[n] = real(ifft(∇[n].*û))
∇v[n] = real(ifft(∇[n].*v̂))
end
ω = Array(Array{Complex{Float64},length(Ω)},length(Ω))
ω[1] = ∇u[2].*∇v[3] - ∇u[3].*∇v[2]
ω[2] = ∇u[3].*∇v[1] - ∇u[1].*∇v[3]
ω[3] = ∇u[1].*∇v[2] - ∇u[2].*∇v[1]
ω = real(ω)
return ω;
end
Out[407]:
In [408]:
ω = vorticity(u[1],u[2],Ω,∇)
size(ω),size(ω[size(ω)[1]])
PyPlot.imshow(ω[3][:,:,1]);
Out[408]: