|
| 1 | +using DrWatson |
| 2 | +@quickactivate "NonlinearDynamicsTextbook" |
| 3 | +include(srcdir("colorscheme.jl")) |
| 4 | +using Random, OrdinaryDiffEq |
| 5 | +import GLMakie |
| 6 | + |
| 7 | +function FitzHugh_Nagumo_ODE(du,u,p,t) |
| 8 | + a, b, d, epsilon, hsq6 = p |
| 9 | + uu = @view u[1,:,:] |
| 10 | + vv = @view u[2,:,:] |
| 11 | + |
| 12 | + # This code can be optimized a lot if time allows |
| 13 | + # Create padded u matrix to incorporate Newman boundary conditions - 9-point stencil |
| 14 | + uuu=[ [uu[2,2] uu[2,:]' uu[2,end-1] ] ; [ uu[:,2] uu uu[:,end-1] ] ; [ uu[end-1,2] uu[end-1,:]' uu[end-1,end-1] ] ] |
| 15 | + diff_term = d .* ( |
| 16 | + 4 .* (uuu[2:end-1,1:end-2] .+ uuu[2:end-1,3:end] .+ |
| 17 | + uuu[3:end,2:end-1] .+ uuu[1:end-2,2:end-1] ) .+ |
| 18 | + uuu[3:end,3:end] .+ uuu[3:end,1:end-2] .+ uuu[1:end-2,3:end] .+ |
| 19 | + uuu[1:end-2,1:end-2] .- 20 .*uu |
| 20 | + ) ./ hsq6 |
| 21 | + |
| 22 | + du[1,:,:] = a .* uu .*(1 .-uu).*(uu.-b) .- vv .+ diff_term |
| 23 | + du[2,:,:] = epsilon .* (uu .- vv ) |
| 24 | +end |
| 25 | + |
| 26 | +# -------------------------- main program ----------------------------- |
| 27 | + |
| 28 | +@time begin |
| 29 | + |
| 30 | +# parameters |
| 31 | +L = 300 # size of domain |
| 32 | +N = 150 # no. of grid points |
| 33 | +h = L/N # spatial steps size |
| 34 | +hsq6 = 6*h*h |
| 35 | + |
| 36 | +# FHN parameters |
| 37 | +a = 3. |
| 38 | +b = 0.2 |
| 39 | +d = 1. # diffusion |
| 40 | +epsilon = 0.01 |
| 41 | + |
| 42 | +# open figure |
| 43 | +kplt = 0 # no. of subplot |
| 44 | + |
| 45 | +allsols = [] |
| 46 | +allts = [] |
| 47 | + |
| 48 | +for ic = 1:2 # initial conditions ( = rows of the plot ) |
| 49 | + |
| 50 | +println("ic=",ic) |
| 51 | + |
| 52 | +uu = zeros(N,N) |
| 53 | +vv = zeros(N,N) |
| 54 | +Random.seed!(13371) |
| 55 | +uu = uu .+ 0.01*randn(N,N) # random perturbation |
| 56 | + |
| 57 | +# local stimulius = concentric waves |
| 58 | +if ic == 1 |
| 59 | + uu[40:41,40:41] .= 0.999 |
| 60 | + uu[75:76,75:76] .= 0.999 |
| 61 | +end |
| 62 | + |
| 63 | +# local stimulius = spiral waves |
| 64 | +if ic == 2 |
| 65 | + uu[1:end,10:11] .= 0.999 |
| 66 | +end |
| 67 | + |
| 68 | +# initial values |
| 69 | +u0 = zeros(2,N,N) |
| 70 | +u0[1,:,:] = uu |
| 71 | +u0[2,:,:] = vv |
| 72 | + |
| 73 | +p = a, b, d, epsilon, hsq6 |
| 74 | + |
| 75 | +if ic == 1 |
| 76 | + saveat = 0.0:2:600 |
| 77 | + tspan = (saveat[1], saveat[end]) |
| 78 | + |
| 79 | + prob = ODEProblem(FitzHugh_Nagumo_ODE, u0, tspan, p) |
| 80 | + sol = solve(prob, Tsit5(); reltol=1e-6, abstol=1e-6, maxiters = 1e7, saveat) |
| 81 | + tvec = sol.t |
| 82 | + uout = [u[1, :, :] for u in sol.u] |
| 83 | + push!(allsols, uout) |
| 84 | + push!(allts, tvec) |
| 85 | + |
| 86 | +elseif ic == 2 |
| 87 | + saveat = 0.0:2.0:300.0 |
| 88 | + tspan = (saveat[1], saveat[end]) |
| 89 | + prob = ODEProblem(FitzHugh_Nagumo_ODE, u0, tspan, p) |
| 90 | + sol = solve(prob, Tsit5(); reltol=1e-6, abstol=1e-6, maxiters = 1e7, saveat) |
| 91 | + |
| 92 | + uout = [u[1, :, :] for u in sol.u] |
| 93 | + # second perturbation line |
| 94 | + u = sol.u[end] |
| 95 | + u[1, 55:73,1:70] .= 0.999 |
| 96 | + saveat = 300.0:2:900 |
| 97 | + tspan = (saveat[1], saveat[end]) |
| 98 | + prob = ODEProblem(FitzHugh_Nagumo_ODE, u, tspan, p) |
| 99 | + sol = solve(prob, Tsit5(); reltol=1e-6, abstol=1e-6, maxiters = 1e7, saveat) |
| 100 | + |
| 101 | + tvec = sol.t |
| 102 | + uout2 = [u[1, :, :] for u in sol.u] |
| 103 | + append!(uout, uout2) |
| 104 | + push!(allsols, uout) |
| 105 | + push!(allts, tvec) |
| 106 | + |
| 107 | +end |
| 108 | +end # ic loop |
| 109 | +end # cputime |
| 110 | + |
| 111 | + |
| 112 | +# %% Animate |
| 113 | +names = ("concentric", "plane_to_spiral") |
| 114 | + |
| 115 | +for i in (1,2) |
| 116 | + name = names[i] |
| 117 | + uout = allsols[i] |
| 118 | + tvec = allts[i] |
| 119 | + |
| 120 | + heatobs = GLMakie.Observable(uout[1]) |
| 121 | + tobs = GLMakie.Observable(0.0) |
| 122 | + titobs = GLMakie.lift(t -> "t = $(t)", tobs) |
| 123 | + |
| 124 | + fig = GLMakie.Figure(resolution = (600, 550)) |
| 125 | + ax = fig[1,1] = GLMakie.Axis(fig; title = titobs) |
| 126 | + hmap = GLMakie.heatmap!(ax, heatobs; colormap = :tokyo, colorrange = (-0.2, 1)) |
| 127 | + cb = GLMakie.Colorbar(fig[1, 2], hmap; width = 20) |
| 128 | + display(fig) |
| 129 | + |
| 130 | + GLMakie.record( |
| 131 | + fig, projectdir("animations", "11", "fitzhugh_$(name).mp4"), |
| 132 | + 1:length(tvec); framerate = 15 |
| 133 | + ) do i |
| 134 | + |
| 135 | + tobs[] = tvec[i] |
| 136 | + heatobs[] = uout[i] |
| 137 | + end |
| 138 | + |
| 139 | +end |
0 commit comments