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.