Inaccuracy with this solution of this PDE from Mathematica NDSolve

Question

I have a problem with this script that solves a nonlinear PDE. Please note that the PDE is defined for the complex function \Psi of x and t whose module is called \phi and argument called \chi.

ClearAll["Global`*"];
M = 1; (* Mass scale*)
A = M*10^-3; (* Gaussian Amplitude*)
L = 10/M; (* Length scale*)
Lb = 100*L; (* Boundary distance*)
tmax = 30*L; (* Maximal time*)
\[Chi]0[x_] := A*Exp[-x^2/(2*L^2)]; (* IC for light field*)
\[Phi]0[x_] := Sqrt[M^2 - 0 + (\[Chi]0'[x])^2/M^2]; (* IC for heavy field*)
sol = NDSolve[{
    -D[\[CapitalPsi][x, t], t, t] + D[\[CapitalPsi][x, t], x, x] + 
      2*((Abs[\[CapitalPsi][x, t]])^2 - M^2 ) \[CapitalPsi][x, t] == 0,
    \[CapitalPsi][x, 0] == \[Phi]0[x]*Exp[I*\[Chi]0[x]],
    (D[\[CapitalPsi][x, t], t] /. t -> 0) == 0,
    D[\[CapitalPsi][-Lb, t], t] == D[\[CapitalPsi][Lb, t], t],
    \[CapitalPsi][-Lb, t] ==  \[CapitalPsi][Lb, t]
    },
   {\[CapitalPsi]}, {x, -Lb, Lb}, {t, 0, tmax}];

y[x_, t_] = \[CapitalPsi][x, t] /. sol[[1]];
\[Chi][x_, t_] = Arg[y[x, t]];
\[Phi][x_, t_] = Abs[y[x, t]];
Plot[{\[Chi]0[x], \[Chi][x, 0]}, {x, -Lb, Lb}, PlotRange -> All]

The last line helps me to compare the function that defines the initial condition for \chi to the evaluation of the evaluation of the argument of \Psi at time t=0 (these two are supposed to match). Unfortunately, the plot gives a different result for the two functions (the real initial condition in blue and the solution in orange):

I don't understand why the solution does not respect the initial condition. Can anybody help me please?

Many thanks.

Answer 1

It looks like the solver has trouble converging. If you use smaller step sizes or other measures to increase numerical stability it should help. For example with MaxStepSize -> 0.5:

ClearAll["Global*"]; M = 1;(*Mass scale*)A = M*10^-3;(*Gaussian Amplitude*)L = 10/M;(*Length scale*)Lb = 100*L;(*Boundary distance*)tmax = 30*L;(*Maximal time*)\[Chi]0[x_] := A*Exp[-x^2/(2*L^2)];(*IC for light field*)\[Phi]0[x_] := Sqrt[M^2 - 0 + (\[Chi]0'[x])^2/M^2];(*IC for heavy field*)sol = NDSolve[{-D[\[CapitalPsi][x, t], t, t] + D[\[CapitalPsi][x, t], x, x] + 2*((Abs[\[CapitalPsi][x, t]])^2 - M^2) \[CapitalPsi][x, t] == 0, \[CapitalPsi][x, 0] == \[Phi]0[x]* Exp[I*\[Chi]0[x]], (D[\[CapitalPsi][x, t], t] /. t -> 0) == 0, D[\[CapitalPsi][-Lb, t], t] == D[\[CapitalPsi][Lb, t], t], \[CapitalPsi][-Lb, t] == \[CapitalPsi][Lb, t]}, {\[CapitalPsi]}, {x, -Lb, Lb}, {t, 0, tmax}, MaxStepSize -> 0.5];

y[x_, t_] = \[CapitalPsi][x, t] /. sol[[1]]; \[Chi][x_, t_] = Arg[y[x, t]]; \[Phi][x_, t_] = Abs[y[x, t]]; Plot[{\[Chi]0[x], \[Chi][x, 0]}, {x, -Lb, Lb}, PlotRange -> All]