differential equations - Solution diverges in periodic PDE

Problem introduced in 11.0.1 and persisting through 11.3

Mathematica version 11 introduces PeriodicBoundaryCondition which is very useful in solving periodic PDE systems. I'm considering using it to solve a 1D time-dependent Schrodinger equation (1D periodic in space + time). But as a first test, I find that the norm of the solution I get is diverging as a function of time, which seems to be incorrect.

Consider a periodic potential

V[x_] := -0.2 (Cos[(π x)/5] + 1)

The eigenstates of this periodic potential can be calculated using NDEigensystem

{vals, funs} = 
 Transpose@
  SortBy[Transpose[
    NDEigensystem[{-(1/2) u''[x] + V[x] u[x], 
      PeriodicBoundaryCondition[u[x], x == -5, 
       TranslationTransform[{10}]]}, 
     u[x], {x} ∈ Line[{{-5}, {5}}], 3, 
     Method -> {"SpatialDiscretization" -> {"FiniteElement", 

         "MeshOptions" -> {"MaxCellMeasure" -> 0.01}}}]], First];

And these are the first five eigenstates

Plot[funs, {x, -5, 5}]

Now as a test, I take the second eigenstates and do a free time-propagation

ufun = NDSolveValue[{I D[u[t, x], t] == -(1/2) D[u[t, x], {x, 2}] + 
     V[x] u[t, x], u[0, x] == funs[[2]],
   PeriodicBoundaryCondition[u[t, x], x == -5, 

    TranslationTransform[{10}]]
   }, u, {t, 0, 5}, {x} ∈ Line[{{-5}, {5}}], 
  Method -> {"FiniteElement", 
    "MeshOptions" -> {"MaxCellMeasure" -> 0.01}}, MaxStepSize -> 0.01];//AbsoluteTiming
(*{13.1789, Null}*)

Since there is no external interaction with the system (there is no term like f[t]*u[t,x] in the equation), the solution should be the same as the initial condition, except for a phase difference. And the norm of the solution should be independent of time. However, for this example, the norm seems to diverge

ListPlot[Table[
  NIntegrate[Abs[ufun[t, x]]^2, {x, -5, 5}], {t, 0, 3, .1}], 
 DataRange -> {0, 3}, PlotRange -> All, Mesh -> Full, 

 FrameLabel -> {"time", "Norm"}]

So why does the numerical solution diverge? I tried to make MaxStepSize and "MaxCellMeasure" smaller, but it doesn't seem to help.

Answer

The results presented in the question suggest an inconsistency between NDEigensystem and NDSolveValue using the new PeriodicBoundaryCondition. This inconsistency can be localized by plotting ufun at various times.

Table[Plot[Evaluate[ReIm[ufun[t, x]]], {x, -5, 5}], {t, 0, 1, .2}]

Evidently, an error is occurring at the boundaries and propagating in. Moreover, spatial derivatives of the solution visibly are not periodic at the boundary, even though the solution itself is. In contrast, derivatives of funs[[2]] do appear to be periodic, if a bit noisy away from the boundary.

Plot[(-(1/2) D[funs[[2]] , {x, 2}] + V[x] funs[[2]]) /. x -> z, {z, -5, 5}]

(The noise can be reduced by decreasing "MaxCellMeasure". Nonetheless, using {"MaxCellMeasure" -> 0.001} in both functions, although painfully slow, reproduces the spurious growth shown in the second plot of the question.) Thus it appears that a bug in NDSolveValue has been introduced in Version 11.

Addendum

Plot[Im[ufun[0, x]], {x, -5, 5}]

I would have expected Im[ufun] to be zero at t == 0, as Im[funs[[2]]] is.

By the way, WorkingPrecision is not an allowed option for NDEigensystem. I hope it will be accommodated in future versions.