# NDSolve 1D Heat equation with NeumannValue poorly satisfies boundary condition

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

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["NDSolveFEM"];
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.

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.