Coupled non-linear differential equations

Question

I have a system of coupled nonlinear differential equations to solve: $$ \frac{\partial m(x,t)}{\partial t}+v(x,t)\frac{\partial m(x,t)}{\partial x}=-\gamma \frac{\partial^2 v(x,t)}{\partial x^2}, \\ \frac{\partial v(x,t)}{\partial t}+v(x,t)\frac{\partial v(x,t)}{\partial x}=-\frac{m(x,t)}{2}.$$

$\gamma$ is a parameter which is strictly positive. The initial conditions are $v(x,0)=0$ and $m(x,0)=0.5*L*\sin(x/L)$ where $\pi L$ is the extension of the grid where I'm working on. Due to the symmetry of the problem, I can as well say that both $m$ and $v$ are vanishing at the boundaries of my grid, i.e. $m(0,t)=m(\pi L,t)=v(0,t)=v(\pi L,t)=0$. That is all I need to solve my problem. After trying to solve this system tt looks that there is a convergence issue around t = 0.788. The code goes as follows:

BackwardEuler = {"FixedStep", Method -> {"ImplicitRungeKutta", 
     "Coefficients" -> "ImplicitRungeKuttaRadauIIACoefficients", 
     "DifferenceOrder" -> 1, "ImplicitSolver" -> {"FixedPoint", 
         AccuracyGoal -> MachinePrecision, 
         PrecisionGoal -> MachinePrecision, 
         "IterationSafetyFactor" -> 1}}};

Tmax = 1.0; L = 60; gamma = 5000; A = 1.0;
xmin = 0.0; xmax = π L; Subscript[ρ, 0] = 1/(1 + A);

mi[x_] := 8 π Subscript[ρ, 0] A L Sin[x/L];

eq1nonstandard = D[m[x, t], t] + v[x, t] D[m[x, t], x] + gamma D[v[x, t], x, x];
eq2nonstandard = D[v[x, t], t] + v[x, t] D[v[x, t], x] + m[x, t]/2;
Vnonstandard = First[v /. NDSolve[{eq1nonstandard == 0, eq2nonstandard == 0,
     v[x, 0] == 0, m[x, 0] == mi[x], v[xmin, t] == 0.0, v[xmax, t] == 0, 
     m[xmin, t] == 0, m[xmax, t] == mi[xmax]}, {v}, {x, xmin, xmax},
     {t, 0, Tmax}, Method -> BackwardEuler, StartingStepSize -> 1/10000]]

But Mathematica complains

Repeated convergence test failure at t == 0.7885000000000001`; unable to continue. >>

Then, I tried to solve the same system of equations in Python using a forward in time/ backward in space finite difference method (explicit method) with a very small spatial and time step. Still, at some point the solution cease to exist. I believe that this is due to the fact that the system is stiff, because if I put $\gamma=0.0$, then (even though the solution diverges as expected from the analytical solution) I get something.

Any suggestions ?

Answer 1

I'm so excited now! I might found a solution for a certain type of PDE related problem! The key point is choosing a odd "DifferenceOrder"!

Let's define a auxiliary function first:

mol[n_, o_] := {"MethodOfLines", 
  "SpatialDiscretization" -> {"TensorProductGrid", "MaxPoints" -> n, 
    "MinPoints" -> n, "DifferenceOrder" -> o}}

Then try this:

Tmax = 1; L = 60; gamma = 5000; A = 1;
xmin = 0; xmax = π L; Subscript[ρ, 0] = 1/(1 + A);
mi[x_] := 8 π Subscript[ρ, 0] A L Sin[x/L];

eq1nonstandard = D[m[x, t], t] + v[x, t] D[m[x, t], x] + gamma D[v[x, t], x, x]; 
eq2nonstandard = D[v[x, t], t] + v[x, t] D[v[x, t], x] + 1/2 m[x, t]; 

Vnonstandard = 
 First[v /. NDSolve[{eq1nonstandard == 0, eq2nonstandard == 0, 
     v[x, 0] == 0, m[x, 0] == mi[x], v[xmin, t] == 0, v[xmax, t] == 0,
      m[xmin, t] == 0, m[xmax, t] == mi[xmax]}, {v}, {x, xmin, xmax}, 
     {t, 0, Tmax}, Method -> mol[25, 9]]]

(* The following line allows you to plot the result easily 
   when NDSolve stops at the half-way. *)
{{xl, xr}, {tl, tr}} = Vnonstandard["Domain"];
Plot3D[Vnonstandard[x, t], {x, xl, xr}, {t, tl, tr}]

Some observation:

1. The number of grid points can't be too large, I guess it's because something similar to this happens.

2. The bigger the "DifferenceOrder" is, the better. This is the result under mol[12, 3]:

and this is under mol[16, 5]:

3. Under some "DifferenceOrder", Mathematica can choose suitable number of grid points automatically. For example mol[Automatic, 9].

Nevertheless, I'm not sure why this Method works, this is just a rare victory among my numerous failures when trying to solve the PDE related problem in this site by trial and error.