NDSolve 1D Heat equation with NeumannValue poorly satisfies boundary condition

Question

Context

I am trying to measure the slope of the stellar cuspide around our MilkyWay's central black hole using the observed distribution of S stars orbiting around it.

Attempt

Close to the Documentation test case, I am trying to solve 1D PDE corresponding to heat diffusion (bound in a box) of a initially Gaussian Distribution a follows:

with

 s = 0.1; a = 0.2; tmax = 5;

the solution is found via

usol = NDSolveValue[
   eqn = {D[u[x, t], t] - 1/5/tmax^2 D[D[u[x, t], x], x] == 
      NeumannValue[0, True],DirichletCondition[
      u[x, t] == (bound[x_] = PDF[NormalDistribution[a, s], x]) // 
       Evaluate, t == 0]
     }, u, {x, -1, 1}, {t, 0, tmax}, AccuracyGoal -> 20, 
   PrecisionGoal -> 20
   ];

I get an as warning

while the answer looks reasonable:

ContourPlot[Evaluate[usol[x, t]], {x, -1, 1}, {t, 0, tmax}, 
 PlotRange -> All, PlotLegends -> Automatic, 
  AspectRatio -> tmax/2]

If I check how well the boundary condition is satisfied:

Plot[usol[x, 0] - bound[x], {x, -1, 1}, PlotRange -> All]

which is not great but ok. If I then change the integration time:

     tmax = 5;

Then the same diagnostic yields

while for tmax=10 it gets even worse:

Question

How can make sure the boundary condition is satisfied whatever the time interval over which I want to integrate?

Attempt to circumvent the problem.

I have tried using explicitly FEM as follows

Needs["NDSolve`FEM`"];
reg = Rectangle[{-1, 0}, {1, tmax}];
reg = ToElementMesh[reg, "MaxBoundaryCellMeasure" -> 0.025, 
   "MeshElementType" -> TriangleElement];
usol2 = NDSolveValue[eqn, u, {x, t} \[Element] reg];
ContourPlot[Evaluate[usol2[x, t]], {x, t} \[Element] reg,  
PlotRange -> All, PlotLegends -> Automatic, 
 AspectRatio -> tmax/2]

Note the noisy contours, while the two solutions differ somewhat:

Plot3D[Evaluate[usol2[x, t] - usol[x, t]], {x, 0, 1}, {t, 0, tmax}, 
 PlotRange -> All]

The fact that I cannot get the solver to provide me with a good asymptotic solution (at late time) is a problem because we need to use this late time limit to constrain the galactic center's cuspide.

Answer 1

Your code doesn't give a good enough result, because NDSolve has chosen pure FiniteElement method to solve the problem i.e. FiniteElement method has been used for discretization in both $t$ and $x$ direction, while the problem is an initial value problem (IVP) in $t$ direction, and FiniteElement method isn't designed for IVP. (I remember user21 has mentioned this in several places, for example here. )

FiniteElement method is triggered in $t$ direction because DirichletCondition has been used to set the initial condition. Then FiniteElement becomes the only choice in $x$ direction, because NDSolve cannot combine TensorProductGrid and FiniteElement at the moment, AFAIK. This topic has been discussed in detail here.

The following is the fixed code:

usol = NDSolveValue[ eqn = {D[u[x, t], t] - 1/5/tmax^2 D[D[u[x, t], x], x] == 0, 
                     u[x, 0] == (bound[x_] = PDF[NormalDistribution[a, s], x])}, 
                     u, {x, -1, 1}, {t, 0, tmax}, 
                       Method -> {MethodOfLines, 
                         SpatialDiscretization -> {FiniteElement, 
                           MeshOptions -> MaxCellMeasure -> 0.01}}];

In this code piece FiniteElement is used only in $x$ direction, and MaxCellMeasure -> 0.01 is added to make the mesh denser. NeumannValue[0, True] is omitted, because zero Neumann value is the default setting of FiniteElement method. This is mentioned in Details section of document of NeumannValue:

When no boundary condition is specified on a part of the boundary $∂Ω$, then the flux term $∇·(-c ∇u-α u+γ)+…$ over that part is taken to be $f=f+0=f+\text{NeumannValue}[0,…]$, so not specifying a boundary condition at all is equivalent to specifying a Neumann 0 condition.

Actually, as mentioned here, even if NeumannValue[0, whatever] is added to the code, it will be simply taken out at parser level.