# 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.

• I updated my answer. Please have a look and let me know what you think about it. – user21 Jun 28 '19 at 6:55

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.

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.

• Thanks for the answer. The new PeriodicBoundaryCondition looks so promising, it's a pity that it doesn't work well. – xslittlegrass Aug 26 '16 at 2:30
• @xslittlegrass I imagine it will in due course. Do you wish to flag this as a bug and report it to Wolfram, Inc., or would you prefer that I do. Your choice. – bbgodfrey Aug 26 '16 at 3:31
• Concerning WorkingPrecision for NDEigensystem this is not on the agenda for quite a while. This would mean to extend the FEM code to non machine precision numbers. This is not a small undertaking and I think there are more important new features to look at first. If disagree, I'd be curious to hear your arguments for it and I may be able re-prioritize this request. – user21 Sep 2 '16 at 21:21
• @user21 Thank you for the unexpected bounty. You are most generous. – bbgodfrey Sep 2 '16 at 21:26
• When I started looking into this, you had already done quite a bit of analysis that I found valuable to get started. So, you are most welcome. – user21 Sep 2 '16 at 21:59

Update

I spend a lot of time to trying to fix this and made several attempt to do so. The idea for a fix is to use "ConstraintMethod" -> "Project" shown in the old part of the post and make that an option on the NDSolve level. The projection method finds a subspace and does the time integration in that subspace and the projects the solution back. While this sounds easy this is very complicated to implement. To give you an idea, this means that if the DE had a WhenEvent then the current solution from the subspace needs to be projected to the original space, customer code runs and the result needs to be projected back to the subspace. Many more details on that level need to thought of and implemented in an efficient manner. At one point I started to wonder, why is this fix so complicated and I stepped back and asked myself why does the current method not handle this as expected. This update is about why the default method behaves as it does; and why that actually is the correct thing in the context of the finite element method.

Consider a time-dependent equation discretized with the finite element method. An initial condition u and implicit Neumann zero boundary conditions on both sides are specified:

ufun = NDSolveValue[{D[u[t, x], t] - D[u[t, x], {x, 2}] == 0,
u[0, x] == Sin[x]}, u, {t, 0, 1}, {x, -\[Pi], \[Pi]},
Method -> {"MethodOfLines",
"SpatialDiscretization" -> {"FiniteElement"}}];


Visualize the solution at various times:

frames = Table[
Plot[ufun[t, x], {x, -\[Pi], \[Pi]}, PlotRange -> {-1, 1}], {t, 0,
1, 0.1}];
ListAnimate[frames, SaveDefinitions -> True]


Note that as soon as the time integration starts at both boundaries, implicit Neumann 0 boundary conditions are satisfied.

Following is the same equation and initial condition as previously and an additional periodic boundary condition that has its source on the left and its target on the right:

ufun = NDSolveValue[{D[u[t, x], t] - D[u[t, x], {x, 2}] == 0,
u[0, x] == Sin[x],
PeriodicBoundaryCondition[u[t, x], x == \[Pi],
Function[X, X - 2 \[Pi]]]}, u, {t, 0, 1}, {x, -\[Pi], \[Pi]},
Method -> {"MethodOfLines",
"SpatialDiscretization" -> {"FiniteElement"}}]


Visualize the solution at various times:

frames = Table[
Plot[ufun[t, x], {x, -\[Pi], \[Pi]}, PlotRange -> {-1, 1}], {t, 0,
1, 0.1}];
ListAnimate[frames, SaveDefinitions -> True]


Note how the value of the implicit Neumann 0 boundary condition at the left is projected to the right. If you had a DirichletCondition on the left, one would expect it to be projected to the right too, right?

So this is the expected behavior for the finite element method. The tensor product grid method behaves differently, as that method does not have implicit boundary conditions:

ufunTPG =
NDSolveValue[{D[u[t, x], t] - D[u[t, x], {x, 2}] == 0,
u[0, x] == Sin[x], u[t, -\[Pi]] == u[t, \[Pi]]},
u, {t, 0, 1}, {x, -\[Pi], \[Pi]},
Method -> {"MethodOfLines",
"SpatialDiscretization" -> {"TensorProductGrid"}}]


Visualize the tensor product grid solution at various times:

frames = Table[
Plot[ufunTPG[t, x], {x, -\[Pi], \[Pi]}, PlotRange -> {-1, 1}], {t,
0, 1, 0.1}];
ListAnimate[frames, SaveDefinitions -> True]


A similar behavior can be achieved with the finite element method by specifying a DirichletCondition on the left and a PeriodicBoundaryCondition:

ufunFEM =
NDSolveValue[{D[u[t, x], t] - D[u[t, x], {x, 2}] == 0,
u[0, x] == Sin[x],
PeriodicBoundaryCondition[u[t, x], x == \[Pi],
Function[X, X - 2 \[Pi]]],
DirichletCondition[u[t, x] == 0, x == -\[Pi]]},
u, {t, 0, 1}, {x, -\[Pi], \[Pi]},
Method -> {"MethodOfLines",
"SpatialDiscretization" -> {"FiniteElement"}}]


Visualize the difference between the finite element and tensor product grid solutions at various times:

frames = Table[
Plot[ufunFEM[t, x] - ufunTPG[t, x], {x, -\[Pi], \[Pi]},
PlotRange -> {-5 10^-4, 5 10^-4}], {t, 0, 1, 0.1}];
ListAnimate[frames, SaveDefinitions -> True]


With all of this I start to think that the issue mention in the question is not a bug but is in fact the correct behavior. I admit that I was baffled by this for a long time. The workarounds presented in the old part of this post are still valid. I have updated the possible issues section of the reference page of PeriodicBoundaryCondition to explain this behavior. I'd like to apologize to @xslittlegrass that it took so long to come up with an explanation. Sorry about that. I still learn about the finite element method too.

Old Post

Let me reorder the equations a bit:

V[x_] := -0.2 (Cos[(\[Pi] x)/5] + 1)
pde = I D[u[t, x], t] + (1/2) D[u[t, x], {x, 2}] - V[x] u[t, x] == 0;

{vals, funs} =
Transpose@
SortBy[Transpose[
NDEigensystem[{-(1/2) u''[x] + V[x] u[x],
PeriodicBoundaryCondition[u[x], x == -5,
TranslationTransform[{10}]]},
u[x], {x} \[Element] Line[{{-5}, {5}}], 3,
Method -> {"SpatialDiscretization" -> {"FiniteElement",
"MeshOptions" -> {"MaxCellMeasure" -> 0.01}}}]], First];
(*Plot[funs[[2]], {x, -5, 5}]*)
init = funs[[2]];


A workaround is the following:

ufunTG = NDSolveValue[{pde, u[0, x] == init, u[t, -5] == u[t, 5]},
u, {t, 0, 5}, {x, -5, 5},
Method -> {"MethodOfLines",
"SpatialDiscretization" -> {"TensorProductGrid"}}];
Table[Plot[Evaluate[ReIm[ufunTG[t, x]]], {x, -5, 5}], {t, 0, 1, .2}]


And:

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


Here are two more possible ideas: One you could enforce the boundary to behave properly as you expect by either specifying:

eqn = {pde, u[0, x] == init,
PeriodicBoundaryCondition[u[t, x], x == -5,
TranslationTransform[{10}]]
, DirichletCondition[u[t, x] == 0, x == 5]};


or

eqn = {I D[u[t, x], t] + (1/2) D[u[t, x], {x, 2}] - V[x] u[t, x] ==
NeumannValue[-10^6 I u[t, x], x == 5], u[0, x] == init,
PeriodicBoundaryCondition[u[t, x], x == -5,
TranslationTransform[{10}]]};


If you need to use the FEM as spatial discretization method then the low level functions can help around the issue:

Needs["NDSolveFEM"]
{dpde, dbc, vd, sd, md} = ProcessPDEEquations @@ fun;
mesh = md["ElementMesh"];
initvals = Flatten[init /. x -> mesh["Coordinates"]];
ifiv = ElementMeshInterpolation[{mesh}, initvals];


We then extract the system matrices:

{lm, sm, dm, mm} = dpde["SystemMatrices"];
deplBCs =
DeployBoundaryConditions[{lm, sm, dm, mm}, dbc,
"ConstraintMethod" -> "Project"]


We construct the intial values and the derivatives there of

pm = "ProjectionMatrix" /. deplBCs["ConstraintMethodData"];
vals = Flatten[sm.(initvals.pm) - lm];
xp = LinearSolve[-dm, vals];
ifxp = ElementMeshInterpolation[{mesh},
NDSolveFEMProcessDiscretizedSolutions[xp, deplBCs]];


Next, the time integration is done:

dif = NDSolveValue[{dm.D[u[t], t] + sm.u[t] == lm,
u[0] == initvals.pm}, u, {t, 0, 5}];


Visualize the solutions:

Table[Plot[
Evaluate[ReIm[
ElementMeshInterpolation[{mesh},
NDSolveFEMProcessDiscretizedSolutions[dif[t], deplBCs]][
x]]], {x, -5, 5}], {t, 0, 1, 0.2}]


Table[Plot[
Evaluate[ReIm[ufunTG[t, x]] -
ReIm[ElementMeshInterpolation[{mesh},
NDSolveFEMProcessDiscretizedSolutions[dif[t], deplBCs]][
x]]], {x, -5, 5}, PlotRange -> All], {t, 0, 1, 0.2}]


ListPlot[Table[
NIntegrate[
Abs[ElementMeshInterpolation[{mesh},
NDSolveFEMProcessDiscretizedSolutions[dif[t], deplBCs]][
x]]^2, {x, -5, 5}], {t, 0, 1, 0.2}], DataRange -> {0, 3},
PlotRange -> All, Mesh -> Full, FrameLabel -> {"time", "Norm"}]