3
$\begingroup$

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.

enter image description here

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

enter image description here

while the answer looks reasonable:

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

enter image description here

If I check how well the boundary condition is satisfied:

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

enter image description here

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

 tmax = 5;

Then the same diagnostic yields

enter image description here

while for tmax=10 it gets even worse:

enter image description here

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]

enter image description here

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]

enter image description here

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.

Improve this question
$\endgroup$
0

1 Answer 1

Reset to default
5
$\begingroup$

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.

Improve this answer
$\endgroup$

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.