# Looking for intuition as to what numerical issue makes this solution of NDSolve blow up unphysically

I am solving the spherically symmetric wave equation in 3 dimensions on $$r\in[0,2\pi]$$ $$(\partial_t^2 - \frac{1}{r^2} \partial_r r^2 \partial_r) \theta(t,r)=0$$ with initial conditions $$\theta(0,r)=e^{-(r-3)^2}$$ and $$\partial_t\theta(0,r)=0$$, imposing boundary conditions $$\partial_r\theta(t,0)=0, \qquad \partial_t(r\theta)(t,6)+\partial_r(6\theta)(t,6)=0.$$ The first is necessary for smoothness at the origin and the second (necessary to solve on a finite domain) should make sure that no information comes form outside the domain.

I include my code below. As you will see, the solution behaves reasonable at first and then starts growing for no reason. I would like to gain some intuition as to what numerically is the problem and how to extend the reasonable behaviour time.

tend = 100.;
rstart = 1/10^100;
rend = 6;

eq = D[θ[t, r], {t, 2}] - 1/r^2 D[r^2 D[θ[t, r], r], r] == 0;
ic = {θ[0, r] == Exp[-(r - 3)^2], (D[θ[t, r], t] /. t -> 0) == 0};
bc = {(D[θ[t, r], r] /. r -> rstart) ==
0, (D[r θ[t, r], t] /. r -> rend) + (D[r θ[t, r], r] /. r -> rend) == 0};

sol = NDSolve[Flatten@{eq, ic, bc}, θ, {t, 0, tend}, {r, rstart, rend}]
Animate[Plot[{θ[t, r] /. sol}, {r, rstart, rend}, PlotRange -> {{rstart, rend}, {-2, 2}}], {t, 0, tend}]


Edit: I'm running version 12.1. Below I include an snapshot of the animation when the solution starts to drift downwards

• Does “DifferenceOrder” -> “Pseudospectral” in the spatial discretization help? Jan 16, 2022 at 19:32
• “As you will see, the solution behaves reasonable at first and then starts growing for no reason.” I don't seem to see it in v12.3.1, can you add your version info and a picture showing the growing? Jan 17, 2022 at 7:05
• @xzczd I've added a picture of a timeslice when the solutions starts to go off. My version is 12.1. I hope that helps Jan 17, 2022 at 14:37
• Here's the result of v12.3.1: i.stack.imgur.com/8STbu.png What if you add Method -> {"MethodOfLines", "SpatialDiscretization" -> {"TensorProductGrid", "MaxPoints" -> 45, "MinPoints" -> 45, "DifferenceOrder" -> 4}, Method -> "Adams"} to NDSolve in v12.1? Also, you may need "DifferentiateBoundaryConditions" -> {True, "ScaleFactor" -> 100}. (To learn more about this option, you can read the following post: mathematica.stackexchange.com/a/127411/1871 ) Jan 17, 2022 at 14:47
• Just tested in v9.0.1, I observed similar phenomenon, and the option combination in my last comment fixes the problem in v9.0.1. A simpler solution is to add MaxStepSize -> {0.05, Automatic}. Do these options work in v12.1? Jan 17, 2022 at 15:01

MethodOfLines with FiniteElement gives a more plausible solution:

bc is changed to NeumannValue (first bc is automatically fullfilled)

teta = NDSolveValue [Flatten@{-((2 r \!$$\*SuperscriptBox[\(\[Theta]$$, TagBox[RowBox[{"(",RowBox[{"0", ",", "1"}], ")"}],Derivative],MultilineFunction->None]\)[t, r] +r^2\!$$\*SuperscriptBox[\(\[Theta]$$,TagBox[RowBox[{"(",RowBox[{"0", ",", "2"}], ")"}],Derivative],MultilineFunction->None]\)[t, r])/r^2) +\!$$\*SuperscriptBox[\(\[Theta]$$,TagBox[
RowBox[{"(", RowBox[{"2", ",", "0"}], ")"}],Derivative],MultilineFunction->None]\)[t, r] ==NeumannValue[-Derivative[1, 0][\[Theta]][t, r], r== 6] , ic },
\[Theta], {t, 0, tend}, {r, rstart, rend},Method -> {"MethodOfLines", "TemporalVariable" -> t,"SpatialDiscretization" -> {"FiniteElement" ,"MeshOptions" -> { "MaxCellMeasure" -> 0.1} }}]

Plot3D[teta[t, r], {t, 0, tend/5}, {r, rstart,rend}, PlotPoints -> 100, PlotRange -> All]


OK, since OP confirms the following solution also works in v12.1, let me turn it into an answer: I observe similar behavior in v9.0.1, and adding

Method -> {"MethodOfLines",
"SpatialDiscretization" -> {"TensorProductGrid", "MaxPoints" -> 45, "MinPoints" -> 45,
"DifferenceOrder" -> 4},
"DifferentiateBoundaryConditions" -> {True, "ScaleFactor" -> 100}, Method -> "Adams"}


is a way to resolve the problem. To understand why DifferentiateBoundaryConditions is added, you can read this post.

Putting my comment into an answer: "DifferenceOrder" -> "Pseudospectral" seems to work.

sol = NDSolve[
Flatten@{eq, ic, bc}, θ, {t, 0, tend}, {r, rstart, rend},
Method -> {"MethodOfLines",
"SpatialDiscretization" -> {"TensorProductGrid",
"DifferenceOrder" -> "Pseudospectral"}}]

Plot3D[θ[t, r] /. sol, {t, 0, tend/5}, {r, rstart, rend},
PlotPoints -> 100, PlotRange -> All]


It throws a warning:

NDSolve::ibcinc: Warning: boundary and initial conditions are inconsistent.