Traffic example with time-dependent external velocity

Import the necessary packages.

using ConservationLawsParticles
using RecursiveArrayTools, DifferentialEquations, Plots

Define the model

• The external velocity is constant and directed to the right.
• The interaction is repulsive at small scales and attractive at long range.
• The mobility models total congestion when the density reaches $1$.

$$ \begin{aligned} V(x) &= [1 + 0.2\sin(t)] \cdot [1 + 0.2\cos(3x+t)] \\ W(x) &= \frac1{(|x|+1)^2} - \frac1{|x|+1} + 0.02x^2 \\ \mathop{\mathrm{mob}}(\rho) &= (1-\rho)_+^2 \end{aligned} $$

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")

Define the ODE problem and the initial conditions

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],)

Solve the ODE

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!()

Plot the density

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")