6
$\begingroup$

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.

Improve this question
$\endgroup$
12
  • 3
    $\begingroup$ reference.wolfram.com/language/FEMDocumentation/tutorial/… $\endgroup$
    – Henrik Schumacher
    Apr 18, 2019 at 17:17
  • 2
    $\begingroup$ Please express the equation in Mathematica format. $\endgroup$
    – bbgodfrey
    Apr 19, 2019 at 12:39
  • 1
    $\begingroup$ What is $\chi$? $\endgroup$
    – Chris K
    Apr 21, 2019 at 2:55
  • $\begingroup$ 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. $\endgroup$
    – Henrik Schumacher
    Apr 21, 2019 at 13:51
  • $\begingroup$ @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. $\endgroup$
    – Riku
    Apr 21, 2019 at 13:59

2 Answers 2

Reset to default
9
+25
$\begingroup$

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}]

fig1

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}]

Figure 2

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}]

Figure 3

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}];

fig4

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 fig5

Improve this answer
$\endgroup$
16
  • $\begingroup$ What version did you use? If it is Version 12, did you not get a message? $\endgroup$
    – user21
    Apr 22, 2019 at 4:14
  • $\begingroup$ @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. $\endgroup$
    – Alex Trounev
    Apr 22, 2019 at 10:44
  • 1
    $\begingroup$ Good, I think the message is correct. It's good that NDSolve warns about this since V12. $\endgroup$
    – user21
    Apr 22, 2019 at 11:10
  • $\begingroup$ This is very interesting. Why do you think the erosion phenomenon happens? And why doesn't anything happen if the initial datum is constant? $\endgroup$
    – Riku
    Apr 22, 2019 at 18:28
  • $\begingroup$ Also, I've asked a related question here: mathematica.stackexchange.com/questions/195776/… $\endgroup$
    – Riku
    Apr 22, 2019 at 18:32
2
$\begingroup$

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

Improve this answer
$\endgroup$
11
  • $\begingroup$ Could you write up an implementation of that technique in this case? $\endgroup$
    – Riku
    Apr 23, 2019 at 13:55
  • $\begingroup$ @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. $\endgroup$
    – Vsevolod A.
    Apr 23, 2019 at 14:33
  • $\begingroup$ @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 $\endgroup$
    – Alex Trounev
    Apr 25, 2019 at 10:35
  • $\begingroup$ @AlexTrounev how does non uniqueness affect my suggestion? $\endgroup$
    – Vsevolod A.
    Apr 25, 2019 at 14:52
  • $\begingroup$ @VsevolodA. Perhaps there can not be used the inverse flow map in a case of the Sobolev vector field. Can you show how it works? Use my first example. $\endgroup$
    – Alex Trounev
    Apr 25, 2019 at 18:25

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.