Kuramoto model.jl

### A Pluto.jl notebook ###
# v0.12.21

using Markdown
using InteractiveUtils

# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
    quote
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : missing
        el
    end
end

# ╔═╡ 8b45e880-86b1-11eb-16ca-893bbbe09a0c
using Distributions, DifferentialEquations, Plots, PlutoUI

# ╔═╡ 5a829f70-86b3-11eb-0c16-4b417e1145d9
begin
	A = @bind astr Select(string.(0:0.1:2); default="1.3") 
	md"""
	# Kuramoto model
	
	```math
	\frac{d\theta_i}{dt} = \omega_i + \frac{K}{N}\sum_{j\ne i}\sin(\theta_j - \theta_i)
	\quad (i=1,2,\ldots,N)
	```
	
	Parameter: ``a`` = $A　``(K = a K_{\mathrm{critical}})``
	"""
end

# ╔═╡ 992e8bf0-86b1-11eb-00b1-0f1436ad9495
function kuramoto!(dθ, θ, param, t)
	N = length(θ)
	K, ω = param 
    for i in 1:N
		dθ[i] = ω[i] + K*mean(sin(θ[j] - θ[i]) for j in 1:N)
	end
end;

# ╔═╡ b0368af0-86b1-11eb-309d-2323de539323
begin
	m, n = 16, 16
	d = Normal(0, 2)
	v = 1.0
	K_c = 2/π/pdf(d, 0)
	tmax = 25.0
end;

# ╔═╡ b701a730-86ba-11eb-1293-0bae6bf3a318
	a = parse(Float64, astr);

# ╔═╡ c9808aae-86b1-11eb-05e9-5fb82716a5a2
begin
	θ₀ = 2π*rand(m, n)
	tspan = (0.0, tmax)
	K = a*K_c
	ω = rand(d, m, n) .+ v
	param = (K, ω)
	prob = ODEProblem(kuramoto!, θ₀, tspan, param)
	sol = solve(prob)
end;

# ╔═╡ 1903eb40-86b2-11eb-11a2-91fdac60377c
anim = @animate for t in [fill(0, 10); 0:0.1:tmax; fill(tmax, 10)]
    plot(; size=(250, 270))
	title!("K = $(a)K_c: t = $t"; titlefontsize=10)
    heatmap!(sin.(sol(t))'; c=:bwr, colorbar=:none, frame=:none)
end;

# ╔═╡ 74870b00-86b2-11eb-0344-9d96d65df59f
gif(anim, "kuramoto$(a).gif", fps=20)

# ╔═╡ de8be7a0-86b7-11eb-1ca8-157bf909d2b9
animh = @animate for t in [fill(0, 10); 0:0.1:tmax; fill(tmax, 10)]
    plot(; size=(400, 270))
	title!("Kuramoto model at K = $(a)K_c: t = $t"; titlefontsize=12)
	bin = range(0, 2π; length=41)
    histogram!(mod2pi.(vec(sol(t))); norm=true, alpha=0.5, label="", bin)
	plot!(; xlabel="θ", ylabel="density")
	plot!(; xlim=(-0.1, 2π+0.1), xtick=(0:π/2:2π, ["0", "π/2", "π", "3π/2", "2π"]))
	plot!(; ylim=(0.0, 1.2), ytick=0:0.1:1.2)
end;

# ╔═╡ de65c200-86b7-11eb-013b-bb8d1ac4a451
gif(animh, "kuramoto$(a)h.gif", fps=20)

# ╔═╡ Cell order:
# ╟─5a829f70-86b3-11eb-0c16-4b417e1145d9
# ╟─74870b00-86b2-11eb-0344-9d96d65df59f
# ╟─de65c200-86b7-11eb-013b-bb8d1ac4a451
# ╠═8b45e880-86b1-11eb-16ca-893bbbe09a0c
# ╠═992e8bf0-86b1-11eb-00b1-0f1436ad9495
# ╠═b0368af0-86b1-11eb-309d-2323de539323
# ╠═b701a730-86ba-11eb-1293-0bae6bf3a318
# ╠═c9808aae-86b1-11eb-05e9-5fb82716a5a2
# ╠═1903eb40-86b2-11eb-11a2-91fdac60377c
# ╠═de8be7a0-86b7-11eb-1ca8-157bf909d2b9