solve set of pde with stiffness

Question

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["DifferentialEquations`NDSolveProblems`"];
Needs["DifferentialEquations`NDSolveUtilities`"];
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!

Answer 1

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