I am working in the traffic flow problem using the Lighthill-Whitham-Richards model together with the Greenshields equation. The equation of that model is this:
$$ \frac{\partial\rho}{\partial t}+v_{max} \left( 1- \frac{2\rho}{\rho_{max}}\right) \frac{\partial\rho}{\partial x}=0. $$
The particular case in which I am having problems with the simulation is the case when a traffic light turns green after being red. For this case, the initial condition is
$$\rho(x,0)= \begin{cases} \rho_{max} & \text{for } x \leq 0 \\ 0 & \text{for } x>0 \end{cases}$$
I have find some conservative method for the problem in the book Modelización de problemas de flujo vehicular (page 80, is in spanish).
UPWIND
$$U_i^{j+1}=U_i^j- \begin{cases} \frac{\Delta t}{\Delta x}\left[ f_i^j -f_{i-1}^j \right] & \text{for } g(U_i^j)>0 \\ \frac{\Delta t}{\Delta x}\left[ f_{i+1}^j -f_i^j \right] & \text{for } g(U_i^j)<0 \end{cases}$$
where $f_i^j=f(U_i^j)=v_{max} U_i^j \left( 1- \frac{U_i^j}{\rho_{max}}\right)$ and $g(u)=f'(u)$.
I use the code from below to do the simulation and the result is that the solution is same as initial condition for any time and that is wrong (solution is posed in in Habermans Book, sec 72, pages 342-350).
L = 2; tfin = 1; k = 0.1; vmax = 12; ρmax = 200;
NN = 200; Δx = L/NN; Δt = Δx/vmax;
tfinal = Ceiling[tfin/Δt, 10];
Print["Δx = ", Δx, " Δt = ",
Δt, " tfinal = ", tfinal, " real final time = ",
N[tfinal*Δt]]
T0[x_] := Piecewise[{{ρmax, x <= 1}, {0, x > 1}}]
T[j_, 0] := T0[x] /. x -> j*Δx;
f[j_, n_] := f[j, n] = vmax*T[j, n]*(1 - T[j, n]/ρmax);
T[j_, n_] := T[j, n] = If[T[j, n - 1] <= ρmax/2, T[j, n - 1] -
(Δt/Δx)*(f[j, n - 1] - f[j - 1, n - 1]), T[j, n
- 1] - (Δt/Δx)*(f[j + 1, n - 1] - f[j, n -
1])];
T[0, n_] := T[0, n] = ρmax;
T[NN, n_] := T[NN, n] = 0;
Table[ListPlot[Table[{Δx*j, T[j, n]}, {j, 0, NN}], Joined ->
True, AxesLabel -> {"x", " "}, PlotLabel -> StringJoin["T[x,t], t=",
ToString[Δt*n]]], {n, 0, tfinal, tfinal/20}]
Other version of Upwind is the one from the article The Conservative Upwind Scheme For Simple Traffic Flow Model (pages 5-6, and yo can see how the solution must to be in page 7):
$$ \rho_i^{n+1}=\rho_i^n - \frac{\Delta t}{\Delta x} \left[ F_{i+\frac{1}{2}}^n - F_{i-\frac{1}{2}}^n \right], $$
where $F_{i+\frac{1}{2}}^n=v_{max} \rho_i^n \left( 1- \frac{\rho_{i+\frac{1}{2}}^n}{\rho_{max}}\right)$ and $\rho_{i+\frac{1}{2}}^n=\frac{\rho_i^n+\rho_{i+1}^n}{2}$.
I use the code from below to do the simulation in this case, but It has been running for more than one hour and I don´t get nothing. I don´t know what the problem can be.
L = 2; tfin = 0.6; k = 0.1; vmax = 12; ρmax = 200;
NN = 200; Δx = L/NN; Δt = Δx/vmax;
tfinal = Ceiling[tfin/Δt, 10];
Print["Δx = ", Δx, " Δt = ",
Δt, " tfinal = ", tfinal, " real final time = ",
N[tfinal*Δt]]
T0[x_] := Piecewise[{{ρmax, x <= 1}, {0, x > 1}}];
T[j_, 0] := T0[x] /. x -> j*Δx;
A1[j_, n_] := A1[j, n] = (T[j, n - 1] + T[j + 1, n - 1])/2;
A2[j_, n_] := A2[j, n] = (T[j - 1, n - 1] + T[j, n - 1])/2;
B1[j_, n_] := B1[j, n] = vmax*(1 - A1[j, n]/ρmax);
B2[j_, n_] := B2[j, n] = vmax*(1 - A2[j, n]/ρmax);
T[j_, n_] := T[j, n] = T[j, n - 1]*(1 - (Δt/Δx)*
(B1[j, n] - B2[j, n]));
T[0, n_] := T[0, n] = ρmax;
T[NN, n_] := T[NN, n] = 0;
Table[ListPlot[Table[{Δx*j, T[j, n]}, {j, 0, NN}], Joined ->
True, AxesLabel -> {"x", " "}, PlotLabel -> StringJoin["T[x,t], t=",
ToString[Δt*n]]], {n, 0, tfinal, tfinal/10}]
Lax-Wendroff
$$ U_{i+\frac{1}{2}}^{j+\frac{1}{2}}=\frac{1}{2}\left[ U_i^j+U_{i+1}^j-\frac{\Delta t}{\Delta x} \left( f_{i+1}^j - f_i^j \right) \right] $$
$$ U_{i-\frac{1}{2}}^{j+\frac{1}{2}}=\frac{1}{2}\left[ U_i^j+U_{i-1}^j-\frac{\Delta t}{\Delta x} \left( f_{i}^j - f_{i-1}^j \right) \right] $$
$$ U_{i}^{j+1}=U_i^j-\frac{\Delta t}{\Delta x} \left[ f_{i+\frac{1}{2}}^{j+\frac{1}{2}} - f_{i-\frac{1}{2}}^{j+\frac{1}{2}} \right] $$
I use the code from below to do the simulation and in this case, it doesn't show nothing again.
L = 2; tfin = 0.2; k = 0.1; vmax = 10; ρmax = 200;
NN = 200; Δx = L/NN; Δt = Δx/vmax;
tfinal = Ceiling[tfin/Δt, 10];
Print["Δx = ", Δx, " Δt = ",
Δt, " tfinal = ", tfinal, " real final time = ",
N[tfinal*Δt]]
T0[x_] := Piecewise[{{ρmax, x <= 1}, {0, x > 1}}];
T[j_, 0] := T0[x] /. x -> j*Δx;
f[j_, n_] := f[j, n] = T[j, n - 1]*vmax*(1 - T[j, n - 1]/ρmax);
A[j_, n_] := A[j, n] = (1/2)*(T[j, n - 1] + T[j + 1, n - 1] -
(Δt/Δx)*(f[j + 1, n] - f[j, n]));
B[j_, n_] := B[j, n] = (1/2)*(T[j, n - 1] + T[j - 1, n - 1] -
(Δt/Δx)*(f[j, n] - f[j - 1, n]));
f1[j_, n_] := f1[j, n] = A[j, n]*vmax*(1 - A[j, n]/ρmax);
f2[j_, n_] := f2[j, n] = B[j, n]*vmax*(1 - B[j, n]/ρmax);
T[j_, n_] := T[j, n] = T[j, n - 1] - (Δt/Δx)*
(f1[j, n] - f2[j, n]);
T[0, n_] := T[0, n] = ρmax;
T[NN, n_] := T[NN, n] = 0;
Table[ListPlot[Table[{Δx*j, T[j, n]}, {j, 0, NN}], Joined ->
True, AxesLabel -> {"x", " "}, PlotLabel -> StringJoin["T[x,t], t=",
ToString[Δt*n]]], {n, 0, tfinal, tfinal/20}]
Questions
- Why are my simulations not working properly? Is it because the scheme (stability, convergence problems) or the problem is in the code?
- Does any one know another more appropriated Finite Difference method for resolving this problem? I know I can use FEM or FVM (Godunov) methods, but I am interested in FDM.
