# Schemes for nonlinear advection equation

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.

First of all, if you just want to solve the equation, NDSolve is enough:

L = 6; tfin = 1; vmax = 12; ρmax = 200;
With[{rho = rho[t, x]},
eq = D[rho, t] + vmax (1 - 2 rho/ρmax) D[rho, x] == 0;
ic = rho == Piecewise[{{ρmax, x <= L/2}}] /. t -> 0;
bc = {rho == ρmax /. x -> 0}];

mol[n_Integer, o_: "Pseudospectral"] := {"MethodOfLines",
"SpatialDiscretization" -> {"TensorProductGrid", "MaxPoints" -> n,
"MinPoints" -> n, "DifferenceOrder" -> o}};
sol = NDSolveValue[{eq, ic, bc, WhenEvent[rho[t, L] > ρmax/100, "StopIntegration"]},
rho, {t, 0, tfin}, {x, 0, L}, Method -> mol[100, 4]];
{t0, tend} = sol["Domain"][[1]];

Manipulate[Plot[sol[t, x], {x, 0, L}, PlotRange -> {-10, ρmax + 10}], {t, 0, tend}]


Then, if you need those FDM schemes for some reason, the following are the reasons why they "don't work".

# Issues in 1st upwind scheme

1. You need to introduce MachinePrecision number to make the calculation not that slow.

2. The first upwind scheme simply can't handle discontinuous initial data properly, please notice $f$ is always $0$ with the given initial condition (i.c.). If the i.c. is smooth enough, it'll have effect.

The following is the fixed code, with a smooth i.c. ρmax/2 Cos[2 Pi/(2 L) x] + ρmax/2:

Clear["*"]
L = 6; tfin = 0.1; vmax = 12; ρmax = 200;
NN = 100; Δx = N@(L/NN); Δt = Δx/vmax;
tfinal = Floor@(tfin/Δt);
Print["Δx = ", Δx, "   Δt = ", \
Δt, "   tfinal = ", tfinal, "   real final time = ",
N[tfinal Δt]]
(* A easier to undestand coding for g: *)
g[tjn_] = D[vmax tjn (1 - tjn/ρmax), tjn];
T0[x_] = ρmax/2 Cos[2 Pi/(2 L) x] + ρmax/2;
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] = T[j,
n - 1] - Δt/Δx Piecewise[{{
f[j, n - 1] - f[j - 1, n - 1], g@T[j, n - 1] > 0}},
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 -> "T[x,t], t=" <> ToString[Δt n // N]], {n, 0, tfinal,
Floor[tfinal/10]}]


# Issues in 2nd upwind scheme

1. Again, you need to introduce MachinePrecision number e.g. modifying Δx to Δx = L/NN // N; to make the calculation not that slow.

2. Your code for the difference formula contains simple mistake. When coding formula it's better to make the formula in the code as similar to the original as possible, in the code below I've included such a code piece, please read it carefully.

3. The time step size turns out to be too large, I've chosen Δt = Δx/vmax/3 instead.

4. There seems to be a typo in the original paper, I believe $F_{i+\frac{1}{2}}^n=v_{max} \rho_i^n \left( 1- \frac{\rho_{i+\frac{1}{2}}^n}{\rho_{max}}\right)$ should be $F_{i+\frac{1}{2}}^n=v_{max} \rho_{i+\frac{1}{2}}^n \left( 1- \frac{\rho_{i+\frac{1}{2}}^n}{\rho_{max}}\right)$. With this modification, the numeric solution is indeed more accurate.

The following is the fixed code:

Clear["*"]
L = 6; tfin = 0.1; vmax = 12; ρmax = 200;
NN = 200; Δx = L/NN // N; Δt = Δx/vmax/3;
tfinal = Floor[tfin/Δt];
Print["Δx = ", Δx, "   Δt = ", \
Δt, "   tfinal = ", tfinal, "   real final time = ",
N[tfinal*Δt]]
T0[x_] := Piecewise[{{ρmax, x <= L/2}, {0, x > L/2}}];
T[j_, 0] := T0[x] /. x -> j*Δx;
(*Your code for the formula,fixed.*)
(*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*A1[j,n](1-A1[j,n]/ρmax);
B2[j_,n_]:=B2[j,n]=vmax*A2[j,n](1-A2[j,n]/ρmax);
T[j_,n_]:=T[j,n]=T[j,n-1]-(Δt (B1[j,n]-B2[j,n]))/Δx;*)
(*My code for the formula:*)
F[i_, -1/2, n_] := F[i - 1, 1/2, n]
F[i_, 1/2, n_] := vmax T[i, 1/2, n] (1 - T[i, 1/2, n]/ρmax)
T[i_, 1/2, n_] := (T[i, n] + T[i + 1, n])/2
T[j_, n_] :=
T[j, n] = T[j,
n - 1] - (Δt (F[j, 1/2, n - 1] - F[j, -1/2, n - 1]))/Δx
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,
Floor[tfinal/10]}]


BTW the following is my implmentation for the scheme:

vmax = 12; rhomax = 200;
L = 6; tend = 0.1;
xpoints = 100;
dx = L/xpoints // N;
dt = dx/vmax/3;
tpoints = tend/dt // Floor;

Clear[rho, F]
F[n_, i_, -1/2] := F[n, i - 1, 1/2]
F[n_, i_, 1/2] := u@rho[n, i, 1/2] rho[n, i, 1/2]
rho[n_, i_, 1/2] := (rho[n, i] + rho[n, i + 1])/2
u[rho_] := vmax (1 - rho/rhomax)
rho[0, i_] = With[{x = i dx}, Piecewise[{{rhomax, x < L/2}}]];
rho[n_, 0] = rhomax;
rho[n_, xpoints] = 0;
(rho[n_Integer, i_] := rho[n, i] = #) &[
Solve[(rho[n + 1, i] - rho[n, i])/dt + (F[n, i, 1/2] - F[n, i, -1/2])/dx == 0,
rho[n + 1, i]][[1, 1, -1]] /. n -> n - 1];
dat = Table[rho[n, i], {n, 0, tpoints}, {i, 0, xpoints}];
dat // ListPlot3D
ListLinePlot[dat[[-1]], DataRange -> {0, L}]


The resulting graph is the same so I'd like to omit it here.

# Issue in Lax-Wendroff scheme

The only issue is you didn't introduce MachinePrecision number. Fixed code:

Clear["*"]
L = 6; tfin = 0.07; vmax = 12; ρmax = 200;
NN = 100; Δx = L/NN // N; Δ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 <= L/2}, {0, x > L/2}}];
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/10}]


• Brilliant, thank you. – MikelBa Jan 16 '18 at 15:54
• In the First Upwind scheme I did not define the $g$ function because $g(\rho )=f'(\rho )=v_{max} \left( 1- \frac{2 \rho}{\rho_{max}} \right)$ and this means that $g(\rho) \begin{cases} \geq 0 & \text{if } 0\leq \rho \leq \rho_{max}/2 \\ <0 & \text{if } \rho_{max}/2 < \rho\leq \rho_{max} \end{cases}$ – MikelBa Jan 16 '18 at 16:09
• @MikelBa Oops, I didn't notice that's the code for $g$, corrected, thanks for pointing out. – xzczd Jan 16 '18 at 16:19
• In your own implementation of the second Upwind scheme you use the command (rho[n_Integer, i_] := rho[n, i] = #) &[Solve[(rho[n + 1, i] - rho[n, i])/dt + (F[n, i, 1/2] - F[n, i, -1/2])/dx == 0, rho[n + 1, i]][[1, 1, -1]] /. n -> n - 1]; What does [[1, 1, -1]] /. n -> n - 1] do? – MikelBa Jan 16 '18 at 16:21
• @MikelBa 1. Part([[]]) doesn't only work on List({}) , for example, a[b[c, d]][[1, 2]] outputs d, in this case, I've just extracted the solution. (You can execute Solve[(rho[n + 1, i] - rho[n, i])/dt + (F[n, i, 1/2] - F[n, i, -1/2])/dx == 0, rho[n + 1, i]][[1, 1, -1]] and observe. ) 2. I've Solved for rho[n + 1, i] but what we need is rho[n, i], so I replaced all of the n with n-1`. – xzczd Jan 16 '18 at 16:38