# Numerical methods to solve a continuity equation

What numerical methods can be used to study the initial value problem for the continuity equation where $$u = u(t, x)$$

$$u_t + \nabla\cdot(\boldsymbol b u) = 0, \qquad t \in [0,T], \quad x=(x_1,x_2) \in \mathbb{R}^2$$ where $$\boldsymbol b = (0,\chi_{\{x_1 \le x_2\}})$$ and how can the solution be plotted using Mathematica?

• Here $$\chi_{x_1 \le x_2}$$ is the characteristic function of the set $$\{(x_1,x_2):x_1\le x_2\}$$. And to fix ideas, we may take $$u(0,\cdot) = 1$$.

As a first step towards a complete solution, we ca take $$\chi_{x_1 \le x_2}$$ as the characteristic function of the set $$\{(x_1,x_2):0 \le x_1\le x_2 \le 1\}$$.

A related (more theoretical) question is on MathOverflow.

• Please express the equation in Mathematica format. – bbgodfrey Apr 19 at 12:39
• What is $\chi$? – Chris K Apr 21 at 2:55
• Is it really necessary that $\chi$ as unbounded support? If $\chi$ were the characteristic function of a compact set, e.g. $\{x \in \mathbb{R}^2 | 0 \leq x_1 \leq x_2 \leq 1 \}$, this would be much easier to do. – Henrik Schumacher Apr 21 at 13:51
• @HenrikSchumacher I'd like to see the case of unbounded support. But as a first step, I'd be happy to see the compact support case too. – Riku Apr 21 at 13:59

To illustrate the problem, I will give an example that differs from the one proposed by Riku. But in this case, numerical instability is better seen. The result is similar to erosion. Perhaps geologists will like this.

b = {1, HeavisideTheta[x - y]}; L = 4; reg =
DiscretizeRegion[Rectangle[{-L, -L}, {L, L}], MaxCellMeasure -> .01];
eq = D[u[t, x, y], t] + Div[b*u[t, x, y], {x, y}] == 0;

ic = u[0, x, y] == Exp[-x^2 - y^2];
bc = {u[t, L, y] == 0, u[t, -L, y] == 0, u[t, x, L] == 0,
u[t, x, -L] == 0};
sol = NDSolveValue[{eq, ic, bc}, u, {x, y} \[Element] reg, {t, 0, 1}];

Table[Plot3D[sol[t, x, y], {x, y} \[Element] reg, Mesh -> None,
PlotRange -> All, PlotLabel -> Row[{"t = ", t}],
AxesLabel -> {"x", "y", ""}], {t, 0, 1, .1}] This code has a message

NDSolveValue::femcscd: The PDE is convection dominated and the result may not be stable. Adding artificial diffusion may help


This is confirmed by the data shown in Figure 1. Use the Sobolev vector field, as suggested by Riku. The following code has no messages, but the solution of the problem shown in Fig. 2 demonstrates similar erosion, as in Fig. 1.

b = {0, HeavisideTheta[x - y]}; L = 4; reg =
DiscretizeRegion[Rectangle[{-L, -L}, {L, L}], MaxCellMeasure -> .01];
eq = D[u[t, x, y], t] + Div[b*u[t, x, y], {x, y}] == 0;

ic = u[0, x, y] == Exp[-x^2 - y^2];
bc = {u[t, L, y] == 0, u[t, -L, y] == 0, u[t, x, L] == 0,
u[t, x, -L] == 0};
sol = NDSolveValue[{eq, ic, bc}, u, {x, y} \[Element] reg, {t, 0, 1}];

Table[Plot3D[sol[t, x, y], {x, y} \[Element] reg, Mesh -> None,
PlotRange -> All, PlotLabel -> Row[{"t = ", t}],
AxesLabel -> {"x", "y", ""}], {t, 0, 1, .1}] Finally, we use the initial data, as suggested by Riku. The following code has no messages, and the solution shown in Figure 3 has no features.

b = {0, HeavisideTheta[x - y]}; L = 4; reg =
DiscretizeRegion[Rectangle[{-L, -L}, {L, L}], MaxCellMeasure -> .01];
eq = D[u[t, x, y], t] + Div[b*u[t, x, y], {x, y}] == 0;

ic = u[0, x, y] == 1;
bc = {u[t, L, y] == 1, u[t, -L, y] == 1, u[t, x, L] == 1,
u[t, x, -L] == 1};
sol = NDSolveValue[{eq, ic, bc}, u, {x, y} \[Element] reg, {t, 0, 1}];

Table[Plot3D[sol[t, x, y], {x, y} \[Element] reg, Mesh -> None,
PlotRange -> All, PlotLabel -> Row[{"t = ", t}],
AxesLabel -> {"x", "y", ""}], {t, 0, 1, .1}] I add a numerical example for the solution that Vsevolod A. proposed. Equation has the form $$u_t+\nabla .(u\vec {b})=0$$ with $$b=(1,\sigma (x-y) )$$ and $$\sigma (s)=$$2*HeavisideTheta[s]-1.

b = {1, 2*(-1/2 + HeavisideTheta[x - y])}; L = 4; reg =
DiscretizeRegion[Rectangle[{-L, -L}, {L, L}], MaxCellMeasure -> .01];
eq = D[u[t, x, y], t] + Div[b*u[t, x, y], {x, y}] == 0;

ic = u[0, x, y] == Exp[-x^2 - y^2];
bc = {u[t, L, y] == 0, u[t, -L, y] == 0, u[t, x, L] == 0,
u[t, x, -L] == 0};
sol = NDSolveValue[{eq, ic, bc}, u, {x, y} \[Element] reg, {t, 0, 1}]; Here we see numerical instability. Now we have to study how this instability arises from the solution in the form of the ODE system. First, we consider how the solution changes as k changes in approximate expression $$\sigma =\frac {k(x-y)}{\sqrt {1+k^2(x-y)^2}}$$

plot[p_] :=
Block[{q = p}, b = {1, q*(x - y)/Sqrt[1 + (q*(x - y))^2]}; L = 4;
reg = DiscretizeRegion[Rectangle[{-L, -L}, {L, L}],
MaxCellMeasure -> .01];
eq = D[u[t, x, y], t] + Div[b*u[t, x, y], {x, y}] == 0;
ic = u[0, x, y] == Exp[-x^2 - y^2];
bc = {u[t, L, y] == 0, u[t, -L, y] == 0, u[t, x, L] == 0,
u[t, x, -L] == 0};
sol = NDSolveValue[{eq, ic, bc},
u, {x, y} \[Element] reg, {t, 0, 1}, Method -> {
"PDEDiscretization" -> {"MethodOfLines",
"SpatialDiscretization" -> {"FiniteElement",
"MeshOptions" -> {"MaxCellMeasure" -> 0.01}}}}];
Plot3D[sol[1, x, y], {x, y} \[Element] reg, Mesh -> None,
PlotRange -> All, PlotLabel -> Row[{"k = ", p}],
AxesLabel -> {"x", "y", ""}]] // Quiet
Table[plot[k], {k, 1, 11, 2}]


Here are the results for t = 1 and different k • What version did you use? If it is Version 12, did you not get a message? – user21 Apr 22 at 4:14
• @user21 of course I got a message NDSolveValue::femcscd: The PDE is convection dominated and the result may not be stable. Adding artificial diffusion may help. And this is shown in Figure 1. – Alex Trounev Apr 22 at 10:44
• Good, I think the message is correct. It's good that NDSolve warns about this since V12. – user21 Apr 22 at 11:10
• This is very interesting. Why do you think the erosion phenomenon happens? And why doesn't anything happen if the initial datum is constant? – Riku Apr 22 at 18:28
• Also, I've asked a related question here: mathematica.stackexchange.com/questions/195776/… – Riku Apr 22 at 18:32

This is a first order PDE which can be reduced to a system of ODE (page 9) and integrated numerically.

I'm going to solve this:

$$\frac{\partial u}{\partial t}+(\frac{\partial}{\partial x},\frac{\partial}{\partial y})\cdot(u,\sigma(x-y))=0$$.

Here $$\sigma$$ is some kind of smooth function with parameter than can give close-enough theta function. A substitution:

$$\tilde{x}=x-y$$

$$\tilde{y}=y$$

Gives

$$\frac{\partial u}{\partial t}+(\frac{\partial}{\partial \tilde{x}},\frac{\partial}{\partial \tilde{y}}-\frac{\partial}{\partial \tilde{x}})\cdot(u,\sigma(\tilde{x}))$$.

Or (without tildes)

$$\frac{\partial u}{\partial t}+\frac{\partial u}{\partial x}(1-\sigma(x))+\frac{\partial u}{\partial y}\sigma(x)=u\sigma'(x)$$

Which results in 4 ordinary differential equations.

$$\frac{dt}{dr}=1$$

$$\frac{du}{dr}=u\sigma'(x)$$

$$\frac{dy}{dr}=\sigma(x)$$

$$\frac{dx}{dr}=1-\sigma(x)$$

Starting from the last:

$$\frac{dx}{1-\sigma(x)}=dr$$

Choosing $$\sigma(x)=\frac{\kappa x}{\sqrt{1+\kappa^2 x^2}}$$

Gives the solution:

$$\frac{2}{3 \kappa } \left(\kappa ^3 x^3+\kappa ^2 x^2 \sqrt{\kappa ^2 x^2+1}+\sqrt{\kappa ^2 x^2+1}+3 \kappa x\right)=r+C_{1}$$

Which mathematica can solve for $$x(r)$$, giving 4 different roots.

The first differential equation gives $$t=r$$ which gives $$x(t)$$.

The second and the third ones are expressed as integrals from analytical expresions:

$$\frac{dy}{dt}=\sigma(x(t))$$

$$y=y_{0}+\int_{t_{0}}^{t}\sigma(x(t))dt$$

$$\frac{du}{dt}=u\sigma'(x(t))$$

$$u=u_{0}\exp(\int_{t_{0}}^{t}\sigma'(x(t))dt$$

Where $$u_{0}=F(y_{0},C_{1})$$ is an arbitrary function determined by initial conditions.

The solution is thus:

$$u=\exp(\int_{t_{0}}^{t}\sigma'(x(t))dt)\cdot F(y-\int_{t_{0}}^{t}\sigma(x(t))dt,\frac{2}{3 \kappa } \left(\kappa ^3 x^3+\kappa ^2 x^2 \sqrt{\kappa ^2 x^2+1}+\sqrt{\kappa ^2 x^2+1}+3 \kappa x\right)-t)$$

DON'T FORGET THAT TO SUBSTITUTE VARIABLES BACK

• Could you write up an implementation of that technique in this case? – Riku Apr 23 at 13:55
• @Riku I suggest you to describe this method for your particular case (three ordinary differential equations with c1=(...), c2=(...) and add that to your question - answers are pretty much guaranteed at this point. – Vsevolod A. Apr 23 at 14:33
• @VsevolodA. This is not the case. There is a Sobolev vector field, so the "failure of uniqueness" is possible, as proved in the article arxiv.org/pdf/1712.03867.pdf – Alex Trounev Apr 25 at 10:35
• @AlexTrounev how does non uniqueness affect my suggestion? – Vsevolod A. Apr 25 at 14:52
• @VsevolodA. Perhaps there can not be used the inverse ﬂow map in a case of the Sobolev vector field. Can you show how it works? Use my first example. – Alex Trounev Apr 25 at 18:25