# solve set of pde with stiffness

I am trying to solve set of PDEs as shown below. Because a[x_, t_] := aL*Exp[-(x - t)^2/L^2] Sin[x - t] when x < -xmax and t<0, I put a[x,0] and D/Dt(a[x,0]) into ICs.

There appear several error messages, one of which is:

NDSolveValue::eerr: Warning: scaled local spatial error estimate of 208.66558109964893at t = 39.666155299563634 in the direction of independent variable x is much greater than the prescribed error tolerance. Grid spacing with 135 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 changed value of precision goal and accuracy goal, but it didn't solve the problem. Thus, I have three equations:

1. Is the ICs correct for a[x,t]?
2. How to avoid the errors?
3. When initional condition of n and ni are changed to n[x,0]==0, ni[x,0]==0, there is another error message, hinting stiffness at around t=4. How to solve it?
Needs["DifferentialEquationsNDSolveProblems"];
Needs["DifferentialEquationsNDSolveUtilities"];
L = 10; b = 10; \[Mu] = 1/1800; aL = 1;
\[Gamma][x_, t_] := ((1 + a[x, t]^2)/(1 - u[x, t]^2))^(1/2);
\[Gamma]i[x_, t_] := ((1 + \[Mu]*a[x, t]^2)/(1 - ui[x, t]^2))^(1/2);
equ = {
D[a[x, t], {x, 2}] - D[a[x, t], {t, 2}] - (4/b^2) a[x, t] -
n[x, t]/\[Gamma][x, t]*a[x, t] == 0,
D[\[Gamma][x, t] u[x, t], t] == -eField[x, t] -
D[\[Gamma][x, t], x],
D[\[Gamma]i[x, t] ui[x, t], t] == -\[Mu]*eField[x, t] -
D[\[Gamma]i[x, t], x],
D[eField[x, t], t] == n[x, t] u[x, t] - ni[x, t] ui[x, t],
D[n[x, t], t] + D[n[x, t]*u[x, t], x] == 0 ,
D[ni[x, t], t] + D[ni[x, t]*ui[x, t], x] == 0
};
ic = {
a[x, 0] == aL*E^(-x^2/L^2) Sin[x],
Derivative[0, 1][a][x,
0] == (-E^(-(x^2/L^2)) *Cos[x] +  E^(-(x^2/L^2)) *2x/L^2 *Sin[x])*aL,
n[x, 0] == 10, ni[x, 0] == 10,
u[x, 0] == 0,  ui[x, 0] == 0, eField[x, 0] == 0
};
xmax = 40; tmax = 100;
{asol, nsol, nisol, usol, uisol, eFsol} =
NDSolveValue[{equ, ic}, {a, n, ni, u, ui, eField}, {x, -xmax,
xmax}, {t, 0, tmax}]


Thank you!

• As mentioned by Alex below, you need to add b.c. to avoid the bcart warning. Please notice bcart warning is a serious problem (much more serious than the eerr mentioned in the question). For more info check this post: mathematica.stackexchange.com/q/73961/1871 May 21, 2020 at 1:22
• Thank you @xzczd. But how to define BCs correctly? For example, it works when bc = {a[-xmin, t] == aL*E^(-xmin^2/L^2) Sin[-xmin], a[xmax, t] == aL*E^(-xmax^2/L^2) Sin[xmax], n[-xmin, t] == n0, ni[-xmin, t] == n0, u[-xmin, t] == 0, ui[-xmin, t] == 0, eField[-xmin, t] == 0};, but not works when bc = {a[-xmin, t] == aL*E^(-xmin^2/L^2) Sin[-xmin], Derivative[1, 0][a][-xmin, t] == aL*Exp[-xmin^2/L^2] (Cos[xmin] - 2 xmin/L^2* Sin[xmin]), n[-xmin, t] == n0, ni[-xmin, t] == n0, u[-xmin, t] == 0, ui[-xmin, t] == 0, eField[-xmin, t] == 0}; May 21, 2020 at 7:49
• The topic is rather complicated, but for your specific problem, add Method -> {MethodOfLines, TemporalVariable -> t}. May 21, 2020 at 8:19

We need to add boundary conditions and Method to solve this system as follows

L = 10; b = 10; \[Mu] = 1/1800; aL = 1; xmin = xmax = 40; tmax = 100;
\[Gamma][x_, t_] := ((1 + a[x, t]^2)/(1 - u[x, t]^2))^(1/2);
\[Gamma]i[x_, t_] := ((1 + \[Mu]*a[x, t]^2)/(1 - ui[x, t]^2))^(1/2);
equ = {D[a[x, t], {x, 2}] - D[a[x, t], {t, 2}] - (4/b^2) a[x, t] -
n[x, t]/\[Gamma][x, t]*a[x, t] == 0,
D[\[Gamma][x, t] u[x, t], t] == -eField[x, t] -
D[\[Gamma][x, t], x],
D[\[Gamma]i[x, t] ui[x, t], t] == -\[Mu]*eField[x, t] -
D[\[Gamma]i[x, t], x],
D[eField[x, t], t] == n[x, t] u[x, t] - ni[x, t] ui[x, t],
D[n[x, t], t] + D[n[x, t]*u[x, t], x] == 0,
D[ni[x, t], t] + D[ni[x, t]*ui[x, t], x] == 0};
ic = {a[x, 0] == aL*E^(-x^2/L^2) Sin[x],
Derivative[0, 1][a][x,
0] == (-E^(-(x^2/L^2)) Cos[x] + 1/L E^(-(x^2/L^2)) x Sin[x])*aL,
n[x, 0] == 10, ni[x, 0] == 10, u[x, 0] == 0, ui[x, 0] == 0,
eField[x, 0] == 0}; bc = {a[-xmin, t] ==
aL*E^(-xmin^2/L^2) Sin[xmin],
a[xmax, t] == aL*E^(-xmax^2/L^2) Sin[xmax], n[-xmin, t] == 10,
n[xmax, t] == 10, ni[-xmin, t] == 10, ni[xmax, t] == 10,
u[-xmin, t] == 0, u[xmax, t] == 0, ui[-xmin, t] == 0,
ui[xmax, t] == 0, eField[-xmin, t] == 0, eField[xmax, t] == 0};

{asol, nsol, nisol, usol, uisol, eFsol} =
NDSolveValue[{equ, ic, bc}, {a, n, ni, u, ui, eField}, {x, -xmin,
xmax}, {t, 0, tmax},
Method -> {"PDEDiscretization" -> {"MethodOfLines",
"SpatialDiscretization" -> {"TensorProductGrid",
"MinPoints" -> 81, "MaxPoints" -> 81,
"DifferenceOrder" -> "Pseudospectral"}}}];


Visualization of solution

var = {asol, nsol, nisol, usol, uisol, eFsol}; var0 = {a, n, ni, u,
ui, eField}; Table[
Plot3D[Re[var[[i]][x, t]], {x, -xmax, xmax}, {t, 0, tmax},
ColorFunction -> Hue, AxesLabel -> Automatic, PlotRange -> All,
Mesh -> None, PlotLabel -> var0[[i]]], {i, Length[var]}]


• Thank you, it works! When the initial conditions are corrected (my mistake), it still works. However, when some parameters are changed, error message arises again. For example, (1) aL=5; (2) n0=0 or n0=0.1; or (3) xmin=xmax=100, and tmax=40. May 20, 2020 at 15:02
• Since n,ni, u, ui, eField only has 1order of derivative for x, BCs are changed (together with modified ICs), which also works. ic = {a[x, 0] == aL*Exp[-x^2/L^2] Sin[x], Derivative[0, 1][a][x, 0] == aL*Exp[-x^2/L^2] (-Cos[x] + 2 x/L^2* Sin[x]), n[x, 0] == n0, ni[x, 0] == n0, u[x, 0] == 0, ui[x, 0] == 0, eField[x, 0] == 0}; bc = {a[-xmin, t] == aL*E^(-xmin^2/L^2) Sin[-xmin], a[xmax, t] == aL*E^(-xmax^2/L^2) Sin[xmax], n[-xmin, t] == n0, ni[-xmin, t] == n0, u[-xmin, t] == 0, ui[-xmin, t] == 0, eField[-xmin, t] == 0};. May 21, 2020 at 7:35
• However, when I changed BCs to below, it does not work bc = {a[-xmin, t] == aL*E^(-xmin^2/L^2) Sin[-xmin], Derivative[1, 0][a][-xmin, t] == aL*Exp[-xmin^2/L^2] (Cos[xmin] - 2 xmin/L^2* Sin[xmin]), n[-xmin, t] == n0, ni[-xmin, t] == n0, u[-xmin, t] == 0, ui[-xmin, t] == 0, eField[-xmin, t] == 0};. Furthermore, there are errors when some parameters are changed, for example, (1) aL=5; (2) n0=0 or n0=0.1; or (3) xmin=xmax=100, and tmax=40.@Alex Trounev May 21, 2020 at 7:52
• @sixpenny It is numerical model. We cant arbitrary change parameters. Do it slowly, step by step like you driving your car. May 21, 2020 at 11:44