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


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:


 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


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

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.


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


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.