Has any work been done about *numerical methods* for the continuity equation 
$$
\partial\_t \rho(x,t) + \operatorname{div} (b(x,t) \rho(x,t)) = 0, \qquad t \in [0,T], \quad x \in \mathbb R^N, 
$$
where $b \in L^1\_tW^{1,p}\_x$?




---


A related equation has been asked on [Mathematica StackExchange](https://mathematica.stackexchange.com/questions/195547/numerical-methods-to-solve-a-continuity-equation).



Numerics for continuity equation with Sobolev vector field
While the standard numerical methods (at least some of them) work even in the Sobolev regularity setting, the analysis of convergence is far from being trivial (due to non-smoothness of $b$). It is addressed for instance in the papers 


1. A. Schlichting, C. Seis: *[Convergence rates for upwind schemes with rough coefficients](https://arxiv.org/abs/1606.09156)* SIAM J. Numer. Anal., 55(2), 2017.
2. A. Schlichting, C. Seis: *[Analysis of the implicit upwind finite volume scheme with rough coefficients](https://arxiv.org/abs/1702.02534)*, Numer. Math. 139(1), 2018.


