In [1]:
# import Pkg; Pkg.add("Flux")
using OrdinaryDiffEq, RecursiveArrayTools, LinearAlgebra, Test, SparseArrays, SparseDiffTools, Sundials
using Plots; pyplot()
using BenchmarkTools
using DifferentialEquations, DiffEqSensitivity
using Flux, DiffEqFlux, Optim
using CuArrays
using DataDrivenDiffEq
using ModelingToolkit
using DiffEqSensitivity

Reaction Diffusion Reformulation

In [2]:
## Diffusion only
In [3]:
# Generate the constants
# p = (0.14, 0.06, 0.035, 0.065) # a,α,ubar,β,D1,D2
p = (3.0, 3.0, 0.035, 0.060) # a,α,ubar,β,D1,D2
N = 100
Ax = Array(Tridiagonal([1.0 for i in 1:N-1],[-2.0 for i in 1:N],[1.0 for i in 1:N-1]))
Ay = copy(Ax)
Ax[2,1] = 2.0
Ax[end-1,end] = 2.0
Ay[1,2] = 2.0
Ay[end,end-1] = 2.0

function RD!(dr,r,p,t)
  D1, D2, F, K = p
  u = r[:,:,1]
  v = r[:,:,2]
  Du = D1*(Ay*u + u*Ax)
  Dv = D2*(Ay*v + v*Ax)
  dr[:,:,1] .= Du .- F*u
  dr[:,:,2] .= Dv .+ K*v
end


Du, Dv, F, K  = p
r0 = zeros(N,N,2)
r0[:,:,1]=ones(N,N)
r0[:,:,1] .+= 0.02 .*rand(N,N)
r0[:,:,2] .+= 0.02 .*rand(N,N)
N2=convert(Int,N/2)
r=6
s=convert(Int,N2-r)
e=convert(Int,N2+r)
r0[s:e, s:e,1] .= 0.5
r0[s:e, s:e,2] .= 0.5

tend=100;
numsaves=20

prob = ODEProblem(RD!,r0,(0.0,tend),p,saveat=1:tend/numsaves:tend)
sol= solve(prob,BS3());
display(plot(sol[:,:,1,end],st=:surface,camera=(-30,30),legend = false));
# display(plot(sol[:,:,2,end],st=:surface,camera=(-30,30),legend = false));
In [4]:
heatmap(sol[:,:,2,end] .> mean(sol[:,:,2,end]),legend = false)
Out[4]:
In [5]:
anim = @animate for i=1:length(sol.t)
    heatmap(sol[:,:,2,i] ,legend = false);
end

gif(anim, fps = 5)

┌ Info: Saved animation to 
│   fn = C:\Users\amira\Dropbox (MIT)\CBA\ReactionDiffusion\01_Code\200503_Search\tmp.gif
â”” @ Plots C:\Users\amira\.julia\packages\Plots\cc8wh\src\animation.jl:98
Out[5]:
In [6]:
anim = @animate for i=1:length(sol.t)
#     heatmap(sol[:,:,2,i],c = :redsblues,legend = false)
    l = @layout [a b]
    p1=plot(sol[:,:,2,i],st=:surface,camera=(-30,30),zlimit=(0,120),legend = false)
end

gif(anim, fps = 5)

┌ Info: Saved animation to 
│   fn = C:\Users\amira\Dropbox (MIT)\CBA\ReactionDiffusion\01_Code\200503_Search\tmp.gif
â”” @ Plots C:\Users\amira\.julia\packages\Plots\cc8wh\src\animation.jl:98
Out[6]:
In [7]:
anim = @animate for i=1:length(sol.t)
    l = @layout [a b]
    p1=plot(sol[:,:,1,i],st=:surface,camera=(-30,30),zlimit=(minimum(sol[:,:,1,:]),maximum(sol[:,:,1,:])),legend = false)
    p2=plot(sol[:,:,2,i],st=:surface,camera=(-30,30),zlimit=(minimum(sol[:,:,2,:]),maximum(sol[:,:,2,:])),legend = false)
    plot(p1, p2,title="Time $(sol.t[i])", layout = l)
end

gif(anim, "diff.gif",fps = 5)

┌ Info: Saved animation to 
│   fn = C:\Users\amira\Dropbox (MIT)\CBA\ReactionDiffusion\01_Code\200503_Search\diff.gif
â”” @ Plots C:\Users\amira\.julia\packages\Plots\cc8wh\src\animation.jl:98
Out[7]:

Variable Diffusion

In [76]:
# Generate the constants
# p = (0.14, 0.06, 0.035, 0.065) # a,α,ubar,β,D1,D2
p = (3.0, 3.0, 0.00, 0.000) # a,α,ubar,β,D1,D2
N = 100
Ax = Array(Tridiagonal([1.0 for i in 1:N-1],[-2.0 for i in 1:N],[1.0 for i in 1:N-1]))
Ay = copy(Ax)
Ax[2,1] = 2.0
Ax[end-1,end] = 2.0
Ay[1,2] = 2.0
Ay[end,end-1] = 2.0

function RD1!(dr,r,p,t)
  D1, D2, F, K = p
  u = r[:,:,1]
  v = r[:,:,2]
  Du = D1*(Ay*u + u*Ax)
  Dv = D2*(Ay*v + v*Ax)
  dr[:,:,1] .= 0
  dr[:,:,2] .= Dv .+ u*v
end


Du, Dv, F, K  = p
r0 = zeros(N,N,2)
# r0[:,:,2]=ones(N,N)

for i in 1:N
    for j in 1:N
#         r0[i,j,1]=(sin(i/N*2*pi)+cos(j/N*4*pi))/1000.0
        r0[i,j,1]=i/N/1000.0
    end
end
        
# r0[:,:,1] .+= 0.02 .*rand(N,N)
# r0[:,:,2] .+= 0.02 .*rand(N,N)
N2=convert(Int,N/2)
r=10
s=convert(Int,N2-r)
e=convert(Int,N2+r)
# r0[s:e, s:e,1] .= 0.5
r0[s:e, s:e,2] .= 0.5

tend=100;
numsaves=20

prob = ODEProblem(RD1!,r0,(0.0,tend),p,saveat=1:tend/numsaves:tend)
sol= solve(prob,BS3());
display(plot(sol[:,:,1,end],st=:surface,camera=(-30,30),legend = false););
display(plot(sol[:,:,2,end],st=:surface,camera=(-30,30),legend = false););
sol.t
Out[76]:
20-element Array{Float64,1}:
  1.0
  6.0
 11.0
 16.0
 21.0
 26.0
 31.0
 36.0
 41.0
 46.0
 51.0
 56.0
 61.0
 66.0
 71.0
 76.0
 81.0
 86.0
 91.0
 96.0
In [77]:
anim = @animate for i=1:length(sol.t)
    l = @layout [a b]
    p1=plot(sol[:,:,1,i],st=:surface,camera=(-30,30),zlimit=(minimum(sol[:,:,1,:]),maximum(sol[:,:,1,:])),legend = false)
    p2=plot(sol[:,:,2,i],st=:surface,camera=(-30,30),zlimit=(minimum(sol[:,:,2,:]),maximum(sol[:,:,2,:])),legend = false)
    plot(p1, p2,title="Time $(sol.t[i])", layout = l)
end

gif(anim, "diff2.gif",fps = 5)

┌ Info: Saved animation to 
│   fn = C:\Users\amira\Dropbox (MIT)\CBA\ReactionDiffusion\01_Code\200503_Search\diff2.gif
â”” @ Plots C:\Users\amira\.julia\packages\Plots\cc8wh\src\animation.jl:98
Out[77]:

Variable diffusion per time

In [139]:
# Generate the constants
# p = (0.14, 0.06, 0.035, 0.065) # a,α,ubar,β,D1,D2
p = (0.5, 0.1, 3.0, 0.06, 0.01,1.0) # a,α,ubar,β,D1,D2
N = 100
Ax = Array(Tridiagonal([1.0 for i in 1:N-1],[-2.0 for i in 1:N],[1.0 for i in 1:N-1]))
Ay = copy(Ax)
# Ax[2,1] = 2.0
# Ax[end-1,end] = 2.0
# Ay[1,2] = 2.0
# Ay[end,end-1] = 2.0

function RD3!(dr,r,p,t)
    D1, D2, D3, p1, p2,p3 = p
    u = r[:,:,1]
    v = r[:,:,2]
    c = r[:,:,3]
    Du = D1*(Ay*u + u*Ax)
    Dv = D2*(Ay*v + v*Ax)
    Dc = D3*(Ay*c + c*Ax)
    dr[:,:,1] .= Du .+ p1*u
    dr[:,:,2] .= Dv .+ p2*v
    dr[:,:,3] .= Dc .+ p3.*v.*u
end


Du, Dv,Dc, p1, p2,p3  = p


r0 = zeros(N,N,3)
# r0[:,:,2]=ones(N,N)

for i in 1:N
    for j in 1:N
        r0[i,j,2]=(sin(i/N*8*pi)+cos(j/N*8*pi))/100.0
    end
end
        
# r0[:,:,1] .+= 0.02 .*rand(N,N)
# r0[:,:,2] .+= 0.02 .*rand(N,N)
N2=convert(Int,N/2)
r=10
s=convert(Int,N2-r)
e=convert(Int,N2+r)
r0[s:e, s:e,1] .= 0.05
r0[s:e, s:e,3] .= 0.5

tend=100;
numsaves=20

prob = ODEProblem(RD3!,r0,(0.0,tend),p,saveat=1:tend/numsaves:tend)
sol= solve(prob,BS3());
display(plot(sol[:,:,1,end],st=:surface,camera=(-30,30),legend = false));
display(plot(sol[:,:,2,end],st=:surface,camera=(-30,30),legend = false));
display(plot(sol[:,:,3,end],st=:surface,camera=(-30,30),legend = false));
sol.t
Out[139]:
20-element Array{Float64,1}:
  1.0
  6.0
 11.0
 16.0
 21.0
 26.0
 31.0
 36.0
 41.0
 46.0
 51.0
 56.0
 61.0
 66.0
 71.0
 76.0
 81.0
 86.0
 91.0
 96.0
In [140]:
anim = @animate for i=1:length(sol.t)
    heatmap(sol[:,:,3,i] .> mean(sol[:,:,3,end]),legend = false)
end
gif(anim, "diff5_2d.gif",fps = 5)

┌ Info: Saved animation to 
│   fn = C:\Users\amira\Dropbox (MIT)\CBA\ReactionDiffusion\01_Code\200503_Search\diff5_2d.gif
â”” @ Plots C:\Users\amira\.julia\packages\Plots\cc8wh\src\animation.jl:98
Out[140]:
In [142]:
anim = @animate for i=1:length(sol.t)
    l = @layout [a b c]
    p1=plot(sol[:,:,1,i],st=:surface,camera=(-30,30),zlimit=(minimum(sol[:,:,1,:]),maximum(sol[:,:,1,:])),legend = false)
    p2=plot(sol[:,:,2,i],st=:surface,camera=(-30,30),zlimit=(minimum(sol[:,:,2,:]),maximum(sol[:,:,2,:])),legend = false)
    p3=plot(sol[:,:,3,i],st=:surface,camera=(-30,30),zlimit=(minimum(sol[:,:,3,:]),maximum(sol[:,:,3,:])),legend = false)

    plot(p1, p2,p3,title="Time $(sol.t[i])", layout = l)
end

gif(anim, "diff5.gif",fps = 5)

┌ Info: Saved animation to 
│   fn = C:\Users\amira\Dropbox (MIT)\CBA\ReactionDiffusion\01_Code\200503_Search\diff5.gif
â”” @ Plots C:\Users\amira\.julia\packages\Plots\cc8wh\src\animation.jl:98
Out[142]:

Non Linearity

In [98]:
# Generate the constants
# p = (0.14, 0.06, 0.035, 0.065) # a,α,ubar,β,D1,D2
p = (0.1, 0.1, 3.0, 0.06, 0.06,0.1) # a,α,ubar,β,D1,D2
N = 100
Ax = Array(Tridiagonal([1.0 for i in 1:N-1],[-2.0 for i in 1:N],[1.0 for i in 1:N-1]))
Ay = copy(Ax)
# Ax[2,1] = 2.0
# Ax[end-1,end] = 2.0
# Ay[1,2] = 2.0
# Ay[end,end-1] = 2.0

function RD4!(dr,r,p,t)
    D1, D2, D3, p1, p2,p3 = p
    u = r[:,:,1]
    v = r[:,:,2]
    c = r[:,:,3]
    Du = D1*(Ay*u + u*Ax)
    Dv = D2*(Ay*v + v*Ax)
    Dc = D3*(Ay*c + c*Ax)
#     dr[:,:,1] .= Du .+ p1*u
#     dr[:,:,2] .= Dv .+ p2*v
#     dr[:,:,3] .= Dc .+ p3.*v.*u +sin.(u).+sin.(v)
    dr[:,:,1] .= Du .+ p1*u
    dr[:,:,2] .= Dv .+ p2*v
    dr[:,:,3] .= Dc .+ p3.*v.*u
end


Du, Dv,Dc, p1, p2,p3  = p


r0 = zeros(N,N,3)
# r0[:,:,2]=ones(N,N)

for i in 1:N
    for j in 1:N
        r0[i,j,2]=(sin(i/N*9*pi)+cos(j/N*5*pi)+i/N+j/N)/100.0
    end
end
        
# r0[:,:,1] .+= 0.02 .*rand(N,N)
# r0[:,:,2] .+= 0.02 .*rand(N,N)
N2=convert(Int,N/2)
r=1
s=convert(Int,N2-r)
e=convert(Int,N2+r)
r0[s:e, s:e,1] .= 0.05
r0[s:e, s:e,3] .= 0.5

tend=100;
numsaves=5

prob = ODEProblem(RD4!,r0,(0.0,tend),p,saveat=1:tend/numsaves:tend)
sol= solve(prob,BS3());
display(plot(sol[:,:,1,end],st=:surface,camera=(-30,30),legend = false));
display(plot(sol[:,:,2,end],st=:surface,camera=(-30,30),legend = false));
display(plot(sol[:,:,3,end],st=:surface,camera=(-30,30),legend = false));
sol.t;
In [99]:
anim = @animate for i=1:length(sol.t)
    heatmap(sol[:,:,1,i] .> mean(sol[:,:,1,i]),legend = false)
end
gif(anim, "diff5_2d2.gif",fps = 5)

┌ Info: Saved animation to 
│   fn = C:\Users\amira\Dropbox (MIT)\CBA\ReactionDiffusion\01_Code\200503_Search\diff5_2d2.gif
â”” @ Plots C:\Users\amira\.julia\packages\Plots\cc8wh\src\animation.jl:98
Out[99]:
In [97]:
anim = @animate for i=1:length(sol.t)
    l = @layout [a b c]
    p1=plot(sol[:,:,1,i],st=:surface,camera=(-30,30),zlimit=(minimum(sol[:,:,1,:]),maximum(sol[:,:,1,:])),legend = false)
    p2=plot(sol[:,:,2,i],st=:surface,camera=(-30,30),zlimit=(minimum(sol[:,:,2,:]),maximum(sol[:,:,2,:])),legend = false)
    p3=plot(sol[:,:,3,i],st=:surface,camera=(-30,30),zlimit=(minimum(sol[:,:,3,:]),maximum(sol[:,:,3,:])),legend = false)

    plot(p1, p2,p3,title="Time $(sol.t[i])", layout = l)
end

gif(anim, "diff8.gif",fps = 5)

┌ Info: Saved animation to 
│   fn = C:\Users\amira\Dropbox (MIT)\CBA\ReactionDiffusion\01_Code\200503_Search\diff8.gif
â”” @ Plots C:\Users\amira\.julia\packages\Plots\cc8wh\src\animation.jl:98
Out[97]:
In [ ]:
In [ ]: