Can someone walk me through solving a PDE numerically?

[cmkluza]

I've recently been making some posts around the web and on this forum attempting to figure out how to use a PDE that models traffic flow in concrete examples. I realize that I have to solve this PDE in order to use it, but I'm sort of lost on how exactly one solves it. The PDE is as follows:

$$\frac{\partial \rho}{\partial t} + \frac{\partial v_{max}(\rho - \frac{\rho^2}{\rho_{max}})}{\partial x} = 0$$
where $\rho$ is density of vehicles, $v$ is velocity, $\rho_{max}$ and $v_{max}$ are max density and velocity respectively, $t$ is time, and $x$ is distance.

That large bit on the top of the second partial is just velocity as a function of density in $q(\rho) = \rho v(\rho)$, where $q$ is flow (number of vehicles per unit time).

From what I've been told, solving this is going to have to be solved numerically, or at least an analytical solution would be way over my head. I've been told that Euler's method could be used here, but I'm not seeing how. I feel like I need to define $\rho(x, t)$ first in order to do anything as far as solving this goes, but no one else has hinted that that is necessary, nor have I ever seen an expression for $\rho(x, t)$ in the different sources I've read through.

Can anyone explain (as simply as possible) how I can go about solving this numerically? Or perhaps there are mathematical utilities I can use to solve this? I do have access to Wolfram's Mathematica program at the moment, which I'm lead to believe can be helpful in solving this, though I'm not overly familiar in Wolfram's language so even with this I don't know how to solve this PDE.

Any help will be greatly appreciated!

 

[BvU]

Your PDE reminds me of the advection equation. It appears in time-dependent particle size distributions when modelling crystallization processes .
I came across the Nagel-Schreckenberg model that looks a lot like what you describe.

Solution requires an intial value $\rho(x,0)$. If that's a constant, the solution is pretty boring.

For a start you can subdivide the road in N sections and impose a $\rho(x,0)$ -- that allows you to calculate $\rho(x,\Delta t)$.

It gets more sophisticated if you let the cars enter with a distribution of speeds, but then you have an extra dimension to deal with.

 

[BvU]

I had fun doing this in Excel
36 sections $\rho(x,0) = 1$ except 7-10: 3 and 21-28: 2 . $v_{\rm \;max} = 0.3, \ \ \rho_{\rm \;max} = 6$ .

I noticed that a block of 5 cars/section doesn't move when in between 1 car/section ( $1-{1\over 6} = 5-{25\over 6}$ ), so there are some limitations to your expression (or I took the steps too big... )

But otherwise you see the 3 cars/section block progress slower than the 2 cars/section , which is nice.

 

[Hesch]

Can anyone explain (as simply as possible) how I can go about solving this numerically?

I'm not quite sure, but think it can be done this way:

Regard the parameter set as a complex number, z, so that z = x + jt.

In the z-plane you can "plot" a landscape ( mountains and valleys ), as for various values of z. You have to reflect all the heights of the landscape in the plane where height = 0.
Now you place a ball somewhere in this landscape, and the ball will roll downwards until level zero is reached. You can let the ball roll, using Newtons method:

z_n+1 = z_n - k*(z_n/z'_n)

I have inserted this "k" to stabilize the iteration. A typical value could be k = 0.1

Say you have two points: z_n and z_n+1, z'_n+1 ≈ ( h( z_n+1) - h(z_n) ) / ( z_n+1 - z_n ).
h(z) = reflected height.

 

[JJacquelin]

This is not a reply to the question about numerical calculus.
In fact the analytical solving is not so hard :
Thanks to the characteristics method of solving PDEs, the general solution expressed on implicite form is :
$\rho=F\left( x-v_{max} \left(\rho-\frac{\rho^2}{\rho_{max} } \right) t \right)$
in which $F$ is any derivable function, to be determined according to the initial condition.
If the initial condition is $\rho(x,0)=f(x)$ with a given function $f(x)$ , then $\rho(x,0)=F\left( x\right)=f(x)$
Thus, the solution expressed on implicite form is : $\rho=f\left( x-v_{max} \left(\rho-\frac{\rho^2}{\rho_{max} } \right) t \right)$
One can explicitly expressed $t$ as a function of$\rho$ :
$t=\frac{x-f^{-1}(\rho)}{v_{max} \left(\rho-\frac{\rho^2}{\rho_{max} } \right)}$
where $f^{-1}$ is the inverse function of $f$.

 

[BvU]

Coming back on this v function: is it really $v = v_{\rm max} (1-{\rho \over \rho_{\rm max} } )$ ?

 

[cmkluza]

Coming back on this v function: is it really $v = v_{\rm max} (1-{\rho \over \rho_{\rm max} } )$ ?

Yep, that's what velocity as a function of density is (as best as I know). Multiplying that by density ($\rho v$) will get you flow ($q$) which is what the second PD in the equation is of, flow with respect to distance. Perhaps it would've been helpful for me to explicitly state that in the OP?

 

[cmkluza]

This is not a reply to the question about numerical calculus.
In fact the analytical solving is not so hard :
Thanks to the characteristics method of solving PDEs, the general solution expressed on implicite form is :
$\rho=F\left( x-v_{max} \left(\rho-\frac{\rho^2}{\rho_{max} } \right) t \right)$
in which $F$ is any derivable function, to be determined according to the initial condition.
If the initial condition is $\rho(x,0)=f(x)$ with a given function $f(x)$ , then $\rho(x,0)=F\left( x\right)=f(x)$
Thus, the solution expressed on implicite form is : $\rho=f\left( x-v_{max} \left(\rho-\frac{\rho^2}{\rho_{max} } \right) t \right)$
One can explicitly expressed $t$ as a function of$\rho$ :
$t=\frac{x-f^{-1}(\rho)}{v_{max} \left(\rho-\frac{\rho^2}{\rho_{max} } \right)}$
where $f^{-1}$ is the inverse function of $f$.

Thanks for your response! I do like analytical methods, but was led to believe that they would be very difficult for this equation. I'm going to try to learn about this characteristics method now. Could you possibly explain more on how I might find the function $f(x)$ in the solution you showed me?

I had fun doing this in Excel
36 sections ρ(x,0)=1\rho(x,0) = 1 except 7-10: 3 and 21-28: 2 . vmax=0.3, ρmax=6v_{\rm \;max} = 0.3, \ \ \rho_{\rm \;max} = 6 .

I noticed that a block of 5 cars/section doesn't move when in between 1 car/section ( 1−16=5−2561-{1\over 6} = 5-{25\over 6} ), so there are some limitations to your expression (or I took the steps too big... )

But otherwise you see the 3 cars/section block progress slower than the 2 cars/section , which is nice.

Thanks a bunch for your help, @BvU, I can't believe you did all that in Excel! I'm sorry if this is obvious, but I'm still not quite wrapping my head around what's going on in that graph. I can see the sections on the x-axis, and the density on the y-axis, but I'm not quite following where those curves come from. Could you possibly show me the input for that graph in Excel?

 

[BvU]

Sure. There's three blocks. Columns are time, rows are distance.

First block: $\rho(x,t)$

Second: $\rho\, v = \rho \,v_{\rm max} (1-{\rho \over \rho_{\rm max} } )$

Third: $\Delta (\rho\, v) \qquad$ (I use $\Delta x = 1$)​

Values in column i of the third block give the change from column i to i + 1 in the first block.

 

[jackmell]

There is a discussion of traffic flow and a PDE describing it in the textbook "Basic Partial Differential Equations" by David Bleecker and George Csordas. Should be available in any university library and may help you better understand the dynamics which involves a cusp catastrophe (a wreck).

 

[cmkluza]

Sure. There's three blocks. Columns are time, rows are distance.

First block: $\rho(x,t)$

Second: $\rho\, v = \rho \,v_{\rm max} (1-{\rho \over \rho_{\rm max} } )$

Third: $\Delta (\rho\, v) \qquad$ (I use $\Delta x = 1$)​

Values in column i of the third block give the change from column i to i + 1 in the first block.

Sorry to bother you even more about this, but I'm still working on understanding all of this. I think I've finally got a good understanding of your table, but I'm wondering if you could explain how you picked your values for initial density at those times ($\rho(x,0)$ = 1, except $x$ = 7-10, 3 and $x$ = 21-28, 2), as well as the maximum density and velocity. I tried inputting some different values for maximum density and volume (225 $\frac{\rm vehicles}{\rm mile}$ and 25 $\frac{\rm miles}{\rm hour}$), but the values became too big very quickly resulting in a lot of errors.

Thanks for all your help thus far!

 

[BvU]

initial density: pretty random choice.

different numbers: then you want to take smaller time steps. So instead of delta t = 1 as I did, you take delta t = 0.01 or something.

 

[depizixuri]

An easy way to calculate numerically is

$$\frac{\Delta\rho}{ \Delta t } = - \frac{\Delta v_{max}(\rho - \frac{\rho^2}{\rho_{max}})}{\Delta x}$$

So, if you know all the values at a time t, and because $\Delta\rho = \rho_{t+1} - \rho_{t}$, then you can calculate $\rho_{t+1}$ this way:$$\rho_{t+1} = \rho_{t} - { \Delta t } * \frac{v_{max}(\rho_t - \frac{\rho_t^2}{\rho_{max}})}{\Delta x}$$

But any numerical solutions will only work if you use a $\Delta t$ restricted to
$v_{max}(\rho_t - \frac{\rho_t^2}{\rho_{max}})< \frac {\Delta x}{\Delta t}$

In other words, if you "speed" crosses your grid on a single time step, then numerical errors will grow exponentially, and your calculation will be ruined. So, if your $\Delta x$ is constant, you need to check that $\Delta t$ is small enough on each timestep.