# Using NDSolve to solve a system of coupled PDEs

I am trying to solve the Gross-Neveu model in one dimension for a specific soliton initial condition. I am trying

X[x_, t_, V_] := (x - V t)/Sqrt[1 - V^2];
T[x_, t_, V_] := (t - V x)/Sqrt[1 - V^2];
eta[x_, t_, V_, Ω_] := Sqrt[1 - Ω^2] X[x, t, V];
alpha[Ω_] := ArcCos[Ω];
beta[V_] := ArcTanh[V];
u0[x_, t_, V_, Ω_] := (Cosh[2 eta[x, t, V, Ω]] - Ω)/(2 (Ω Cosh[2 eta[x, t, V, Ω]] - 1)) Sqrt[1 - Ω^2]/Cosh[eta[x, t, V, Ω] - I/2 (alpha[Ω] - Pi)] Exp[I Ω T[x, t, V] + beta[V]/2];
v0[x_, t_, V_, Ω_] := -(Cosh[2 eta[x, t, V, Ω]] - Ω)/(2 (Ω Cosh[2 eta[x, t, V, Ω]] - 1)) Sqrt[1 - Ω^2]/Cosh[eta[x, t, V, Ω] + I/2 (alpha[Ω] - Pi)] Exp[I Ω T[x, t, V] - beta[V]/2];
V = 0;
Ω = -0.1;
L = 20;
tmax = 20;
gnm  = {I (D[u[x, t], {t, 1}] + D[u[x, t], {x, 1}]) + v[x, t] +
2 v[x, t] (u[x, t] Conjugate[v[x, t]] +
v[x, t] Conjugate[u[x, t]]) == 0,
I (D[v[x, t], {t, 1}] - D[v[x, t], {x, 1}]) + u[x, t] +
2 u[x, t] (u[x, t] Conjugate[v[x, t]] +
v[x, t] Conjugate[u[x, t]]) == 0,
u[x, 0] == u0[x, 0, V, Ω],
v[x, 0] == v0[x, 0, V, Ω],
Derivative[1, 0][u][L, t] == Derivative[1, 0][u][-L, t],
Derivative[1, 0][v][L, t] == Derivative[1, 0][v][-L, t]
};
sol = NDSolve[gnm, {u, v}, {x, -L, L}, {t, 0, tmax},
Method -> {"MethodOfLines", "SpatialDiscretization" -> {"TensorProductGrid"}}] ​


I am getting the following error though

NDSolve::eerr: Warning: scaled local spatial error estimate of 447.1246272034737 at t = 20. in the direction of independent variable x is much greater than the prescribed error tolerance. Grid spacing with 179 points may be too large to achieve the desired accuracy or precision. A singularity may have formed or a smaller grid spacing can be specified using the MaxStepSize or MinPoints method options.


I think I am making a silly mistake but I just can't see it.

• Hello, welcome to Mathematica.SE! Then please show us the complete code sample, without the definition of V, L, etc. we can't reproduce your problem. BTW, InterpolatingFunction can't be transferred with a simple copy&paste even inside Mathematica. And you may want to have a look at this tool for your code formatting. – xzczd Oct 16 '14 at 5:08
• It might take me a while to get the tool working. Thank you for your speedy reply. The complete code sample is now there. – Robert Moerman Oct 16 '14 at 7:21

You need the magic of "Pseudospectral":

sol = With[{nxy = 250},
NDSolve[gnm, {u, v}, {x, -L, L}, {t, 0, tmax},
Method -> {"MethodOfLines",
"SpatialDiscretization" -> {"TensorProductGrid",
"MaxPoints" -> nxy, "MinPoints" -> nxy, "DifferenceOrder" -> "Pseudospectral"}}]];

Animate[Plot[Through[{Re, Im}@#] & /@ {u[x, t], v[x, t]} /. sol
// Evaluate, {x, -L, L}, PlotRange -> 1.5], {t, 0, tmax}] • Thanks for the help. Although it works, I expected the initial conditions to preserve their shape in time since $u0[x_, t_, V_, Ω_]$ and $v0[x_, t_, V_, Ω_]$ are analytic travelling soliton solutions. Do you have any thoughts on this? – Robert Moerman Oct 16 '14 at 11:34
• @RobertMoerman You mean the current plot doesn't match the physical prediction and the plot should be two traveling waves from the origin? If so, My tools at hand can't solve this problem. I just noticed that, strictly speaking, the i.c. and b.c. aren't consistent, can this be a reason for the undesired behavior? BTW, you don't need to accept an answer so fast, it'll be better to wait for 24 hours or more so others may give better answers. And, I just checked the wiki of Gross-Neveu model but got nothing valid, maybe you can talk about this a little more in your question. – xzczd Oct 16 '14 at 12:25