Import the necessary packages.
using ConservationLawsParticles
using RecursiveArrayTools, DifferentialEquations, Plots
@time_independent V(x) = 1.
@time_independent W(x) = 2(exp(abs(x)/4) + exp(-2abs(x))) - 4
@time_independent W′(x) = 2sign(x) * (exp(abs(x)/4)/4 - 2exp(-2abs(x)))
mob(ρ) = max(1 - ρ, 0)
smodel = SampledModel((V,), ((W′,),), (mob,))
imodel = IntegratedModel((V,), ((W,),), (mob,))
plot(W, -5, 5, title="Interaction", label="W")
plot!(W′, -5, 5, 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., 5.)
sprob = ODEProblem(velocities_gen!, x0, tspan, smodel)
iprob = ODEProblem(velocities_gen!, x0, tspan, imodel)
;
abstol, reltol = 1e-7, 1e-7
@time ssol = solve(sprob, BS5(); abstol=abstol, reltol=reltol)
@time isol = solve(iprob, BS5(); abstol=abstol, reltol=reltol)
length.((ssol, isol))
0.015444 seconds (1.96 k allocations: 646.656 KiB) 0.011731 seconds (3.02 k allocations: 1.132 MiB)
(50, 50)
plot(title="Traffic trajectories, sampled (blue) and integrated (red)", legend=false)
plot!(ssol; color=:blue)
plot!(isol; color=:red)
plot!(xlabel="time", ylabel="space")
savefig("plots/trajectories.png")
plot!()
Let us first compute a more refined solution.
n = 101
x0 = ArrayPartition(vcat(range(-1, -.5, length=n), range(.5, 1, length=n)))
sprob = ODEProblem(velocities_gen!, x0, tspan, smodel)
iprob = ODEProblem(velocities_gen!, x0, tspan, imodel)
abstol, reltol = 1e-7, 1e-7
@time ssol = solve(sprob, BS5(); abstol=abstol, reltol=reltol)
@time isol = solve(iprob, BS5(); abstol=abstol, reltol=reltol)
length.((ssol, isol))
0.507735 seconds (5.36 k allocations: 5.627 MiB) 0.486930 seconds (8.30 k allocations: 10.697 MiB)
(137, 136)
Now we can plot an animation of the density.
anim = @animate for t in range(tspan..., step=1/24)
plot_density(isol(t); legend=false, xrange=(-1, 5), yrange=(0, 1),
title="Traffic density", xlabel="position", ylabel="density")
end
gif(anim, "plots/traffic.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(isol(t)); legend=false, xrange=(-2, 2), yrange=(0, 1),
title="Traffic density (centered)", xlabel="position", ylabel="density")
end
gif(anim, "plots/traffic_centered.gif")