Import the necessary packages.
using ConservationLawsParticles
using RecursiveArrayTools, DifferentialEquations, Plots
V(t, x) = (1 + 0.2sin(t)) * (1 + 0.2cos(3x + t))
@time_independent W(x) = 1/(abs(x)+1)^2 - 1/(abs(x)+1) + .05x^2
mob(ρ) = max(1 - ρ, 0)^2
model = IntegratedModel((V,), ((W,),), (mob,))
plot(W, -4, 4, title="Interaction", label="W")
The initial condition approximates $\rho_0 = 1_{[-1,-1/2]} + 1_{[1/2,1]}$.
n = 25
x0 = ArrayPartition(vcat(range(-1, -.5, length=n), range(.5, 1, length=n)))
tspan = (0., 20.)
prob = ODEProblem(velocities_gen!, x0, tspan, model)
ODEProblem with uType ArrayPartition{Float64, Tuple{Vector{Float64}}} and tType Float64. In-place: true timespan: (0.0, 20.0) u0: ([-1.0, -0.9791666666666666, -0.9583333333333334, -0.9375, -0.9166666666666666, -0.8958333333333334, -0.875, -0.8541666666666666, -0.8333333333333334, -0.8125 … 0.8125, 0.8333333333333334, 0.8541666666666666, 0.875, 0.8958333333333334, 0.9166666666666666, 0.9375, 0.9583333333333334, 0.9791666666666666, 1.0],)
abstol, reltol = 1e-7, 1e-7
@time sol = solve(prob, BS5(); abstol=abstol, reltol=reltol)
length(sol)
0.007236 seconds (7.70 k allocations: 2.900 MiB)
130
plot(title="Traffic trajectories", legend=false)
plot!(sol; color=:blue)
plot!(xlabel="time", ylabel="space")
savefig("plots/trajectories2.png")
plot!()
Let us first compute a more refined solution.
n = 201
x0 = ArrayPartition(vcat(range(-1, -.5, length=n), range(.5, 1, length=n)))
tspan = (0, 20.)
prob = ODEProblem(velocities_gen!, x0, tspan, model)
abstol, reltol = .5e-7, .5e-7
@time sol = solve(prob, BS5(); abstol=abstol, reltol=reltol)
length(sol)
1.773375 seconds (43.16 k allocations: 104.561 MiB)
721
Now we can plot an animation of the density.
anim = @animate for t in range(tspan..., step=1/24)
plot_density(sol(t); legend=false, xrange=(-1, 15), yrange=(0, 1),
title="Traffic density", xlabel="position", ylabel="density")
end
gif(anim, "plots/traffic2.gif")
We can translate the particles so that the barycenter stays at the origin.
barycenter(x::AbstractVector) = (sum(x) - (x[1] + x[end]) / 2) / (length(x) - 1)
recenter(x::AbstractVector) = x .- barycenter(x)
recenter (generic function with 1 method)
anim = @animate for t in range(tspan..., step=1/24)
plot_density(recenter(sol(t)); legend=false, xrange=(-3, 4), yrange=(0, 1),
title="Traffic density (centered)", xlabel="position", ylabel="density")
end
gif(anim, "plots/traffic_centered2.gif")