Coupled non-linear differential equations

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",
"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 ?

• Hi ! This is a strictly Mathematica.SE, so asking about Python is off-topic. – Sektor Mar 20 '15 at 17:48
• Hi! Ok, sorry. Anyway, the python part of the question is not relevant. It was just to put forward that maybe the system is stiff. – Ricoleto Mar 20 '15 at 17:54
• Seems that this hasn't been done yet: Hello, welcome to Mathematica.SE! I suggest the following: 1) As you receive help, try to give it too, by answering questions in your area of expertise. 2) Read the faq! 3) When you see good questions and answers, vote them up by clicking the gray triangles, because the credibility of the system is based on the reputation gained by users sharing their knowledge. Also, please remember to accept the answer, if any, that solves your problem, by clicking the checkmark sign! – xzczd Mar 21 '15 at 7:55

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

1. 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.

• Wow! That really looks amazing. Thanks! How did you figure out that method?! Nonetheless, it still looks that there is an issue if you increase the time Tmax to higher values such as Tmax=10. Don't you have that problem? – Ricoleto Mar 21 '15 at 23:59
• @Ricoleto As mentioned above, it's just a result of trial and error. Sadly it seems not to be able to handle Tmax=10`, and I haven't think out a corresponding solution yet. I think it might help if you add some background information about your equations to your question. BTW, I just noticed that the b.c. in your code is different from what you claimed in the text i.e. $m(0,t)=m(\pi L,t)=v(0,t)=v(\pi L,t)$, is it intended? – xzczd Mar 22 '15 at 11:33
• Basically, $m(x,t)$ is the integral of some density, while $v(x,t)$ corresponds to the velocity. If you forget about the $\gamma$-term, those equations are just the conservation equation and Euler equation in disguise. About my boundary equations, sorry I just forgot to say that all of them are equal to 0 (that's why they are all equal). – Ricoleto Mar 22 '15 at 11:39