Numerical method for Coupled PDE (system of conservation laws) I have the following PDE's:
$\frac{{\partial \rho }}{{\partial t}} =  - \frac{{\partial \rho }}{{\partial x}}u - \rho \frac{{\partial u}}{{\partial x}}$
$\frac{{\partial u}}{{\partial t}} =  - u\frac{{\partial u}}{{\partial x}} - \frac{\partial }{{\partial x}}\left[ {\frac{1}{{\sqrt \rho  }}\frac{{{\partial ^2}(\sqrt \rho  )}}{{\partial {x^2}}}} \right]$
where
$\rho (x,t)\,\,\,,\,\,\,u(x,t)$
the boundary conditions are given for: 
$\rho (x,t = 0),\,\,\,\,u(x,t = 0)$
I want to solve the equations numerically and see how $\rho$ and $u$ evolve with time.
I'm familiar with the finite difference and finite valume methods, but since the two equations are coupled I don't know how to start.
Can anyone suggest a practical numerical method for this kind of problem ?
Thank you !
 A: The system of equations rewrites as a system of conservation laws with relaxation
$$
\frac{\partial}{\partial t} \boldsymbol{U} + \frac{\partial}{\partial x} \boldsymbol{f}(\boldsymbol{U}) = \boldsymbol{r}(\boldsymbol{U}) \, ,
$$
where $\boldsymbol{U} = (\rho,u)^\top$, $\boldsymbol{f}(\boldsymbol{U}) = (\rho u, \frac{1}{2}u^2)^\top$ and $\boldsymbol{r}(\boldsymbol{U}) = (0, -\frac{\partial}{\partial x}[\frac{1}{\sqrt{\rho}} \frac{\partial^2\sqrt{\rho}}{\partial x^2}])^\top$.
An efficient strategy to solve numerically this kind of equations is splitting (e.g. Godunov splitting or Strang splitting), which consists in solving successively the propagative part
$$
\frac{\partial}{\partial t} \boldsymbol{U} + \frac{\partial}{\partial x} \boldsymbol{f}(\boldsymbol{U}) = \boldsymbol{0}
$$
by a well-suited finite-volume method, and the relaxation part
$$
\frac{\partial}{\partial t} \boldsymbol{U} = \boldsymbol{r}(\boldsymbol{U})
$$
at each time-step $t \leftarrow t+\Delta t$. Otherwise, a straightforward explicit discretization
$$
\frac{\boldsymbol{U}_i^{n+1}-\boldsymbol{U}_i^{n}}{\Delta t} +  \frac{\hat{\boldsymbol{f}}_{i+1/2}^n - \hat{\boldsymbol{f}}_{i-1/2}^n}{\Delta x} = \hat{\boldsymbol{r}}_i^n
$$
of the full system may be used, where $\boldsymbol{U}_i^{n} \approx \boldsymbol{U}(i\Delta x, n\Delta t)$ denotes the numerical solution, $\hat{\boldsymbol{f}}_{i+1/2}^n$ is the numerical flux from a finite-volume method (e.g. Lax-Friedrichs or Godunov), and $\hat{\boldsymbol{r}}_i^n \approx \boldsymbol{r}(\boldsymbol{U}(i\Delta x, n\Delta t))$. Thus, one obtains the time-stepping formula
$$
\boldsymbol{U}_i^{n+1} = \boldsymbol{U}_i^{n} - \frac{\Delta t}{\Delta x} \left( \hat{\boldsymbol{f}}_{i+1/2}^n - \hat{\boldsymbol{f}}_{i-1/2}^n\right) + \Delta t \, \hat{\boldsymbol{r}}_i^n \, .
$$
The latter approach is often not as efficient as splitting, and may have stability issues. For further reading, a reference book dedicated to this kind of systems is LeVeque, 2002.
