6
$\begingroup$

Consider the Burgers equation $$\partial_t u + \partial_x\left(u^2/2 \right) = 0, \quad u(0, x) = u_0(x).$$

eq = D[u[t, x], t] + D[u[t,x]^2/2, x] == 0

How can I use Mathematica to get the (unique entropy) solution that has the Cantor staircase function as initial data $u_0 = F_C$?

And how can I plot this solution?

Improve this question
$\endgroup$
7
  • 1
    $\begingroup$ I don't think that this is possible becaus $\partial_x(u^2/2)$ is a singular measure (neither absolutely continuous with respect to the Lebesgue measure nor a countable sum of Dirac measures). But you may try to use NDSolve with the first few iterates of the piecewise-linear functions that converge to the Cantor staircase function. In the (for me) best imaginable universe, the solutions might converge in a very weak sense to your desired solution so that you might get an impression of how it looks like. $\endgroup$
    – Henrik Schumacher
    Apr 29, 2019 at 12:08
  • 1
    $\begingroup$ Where are you getting stuck -- solving the Burgers equation? setting up the Cantor's staircase? Have you tried anything? $\endgroup$
    – Chris K
    Apr 29, 2019 at 12:19
  • $\begingroup$ Actually, since the finite iterates $F_n$ that are usually used to approximate $F_C$ are piecewise-linear and monotonically increasing, it should be possible to solve Burger's equation analytically in a piewise polynomial way. Due to shocks, the polynomial patches have to be truncated also in time and I expect that first time $T_n$ for which the a shock emerges arises with $u_0 = F_n$ converges to $T_n \to 0$ for $n \to \infty$. So, things will become really complicated. $\endgroup$
    – Henrik Schumacher
    Apr 29, 2019 at 12:21
  • 2
    $\begingroup$ @Zyl If you have NDSolve working for other initial conditions, post your code here so people don't have to reinvent the wheel. $\endgroup$
    – Chris K
    Apr 29, 2019 at 19:16
  • 1
    $\begingroup$ A Burger equation is part of the nonlinear FEM verification tests: Convection—FEM-NL-Transient-1D-Single-Convection-0001 or in the help system under FEMDocumentation/tutorial/NonlinearFiniteElementVerificationTests#\ 1129285755 $\endgroup$
    – user21
    Apr 30, 2019 at 4:21

1 Answer 1

Reset to default
9
$\begingroup$

Here is how to obtain the "$n$-th approximating function of CantorStaircase.

ClearAll[f];
f[0] = x \[Function] x;
f[n_Integer?Positive] := f[n] = x \[Function] Evaluate[PiecewiseExpand[
      Piecewise[{
        {0, x <= 0},
        {1/2 f[n - 1][3 x], 0 <= x <= 1/3},
        {1/2, 1/3 <= x <= 2/3},
        {1/2 + 1/2 f[n - 1][3 x - 2], 2/3 <= x <= 1},
        {1, 1 <= x}
        }]]];

Dislcaimer: As nonlinear FEM is involved in the following, this will work only with _Mathematica_ version 12 or later.

Next step is to solve the PDE numerically. This is a hyperbolic PDE of firsts order which is notorious for its shocks. So we use some diffusion to obtain a stable numerical scheme (so ϵ converging to 0, we will obtain the viscousity sol). The second, very important trick is to general a spacial discretization that has the nondifferentiable points of the initial conditions contained in its vertex list.

Needs["NDSolve`FEM`"]
sign = 1;
ϵ = 0.0001;
T = 1;
n = 5;
nElements = 10 3^n;
Ω = ToElementMesh[
   "Coordinates" -> Partition[Subdivide[-1., 2., nElements], 1],
   "MeshElements" -> {LineElement[Partition[Range[nElements + 1], 2, 1]]},
   "MeshOrder" -> 1
   ];
sol = NDSolveValue[
  {
   D[u[t, x], t] - ϵ D[u[t, x], x, x] + D[u[t, x]^2/2, x] == 0,
   u[0, x] == sign f[n][x]
   },
  u, {t, 0, T}, x ∈ Ω
  ];

Plot3D[sol[t, x], {x, -1/2, 3/2}, {t, 0, T},
 AxesLabel -> {"x", "t"},
 PlotRange -> sign {-0.1, 1.1},
 NormalsFunction -> None,
 PlotPoints -> {200, 200},
 ViewPoint -> {- sign 1.3, -2.4, 2.}
 ]

enter image description here

We can nicely see the rarefaction waves. Shock waves can only be observed after setting sign to a negative value; here is the plot for sign = -1:

enter image description here

Improve this answer
$\endgroup$
6
  • $\begingroup$ @Zyl sign flips the initial condition from $F_C$ to $-F_C$. Rarefaction waves are the "opposite" of shocks. Please, inform yourself on Burger's equation for more details. $\endgroup$
    – Henrik Schumacher
    Apr 29, 2019 at 21:14
  • $\begingroup$ Thank you. In addition to the previous three questions ( 1. Do you get anything by setting the initial data as "CantorStaircase[x]"? 2. If not, is it possible to have a table with multiple plots (for different and bigger values of n? 3. Can you get from Mathematica an explicit expression of the solution at time 1? It kind of looks like a step function.), a fourth one occurred: 4. What happens for longer times, for example, up to $t=3$? $\endgroup$
    – user64094
    Apr 29, 2019 at 21:25
  • $\begingroup$ 1. Just try it. 2. You're free to try it yourself. 3. Not with NDSolve. As I told you before, it is principally possible to solve the ODE exactly for f[n] as initial condition. But that's simply to much work for me. Try it yourself! 4. Again, you are free to try it yourself. $\endgroup$
    – Henrik Schumacher
    Apr 29, 2019 at 21:46
  • $\begingroup$ I'd gladly try myself, but I don't have a newer version of Mathematica to make the nonlinear FEM run. $\endgroup$
    – user64094
    Apr 29, 2019 at 21:50
  • $\begingroup$ Maybe you have access to Mathematica online? That runs already under version 12. $\endgroup$
    – Henrik Schumacher
    Apr 29, 2019 at 21:54

Your Answer

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