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

enter image description here

  • $\begingroup$ Does “DifferenceOrder” -> “Pseudospectral” in the spatial discretization help? $\endgroup$
    – Michael E2
    Jan 16, 2022 at 19:32
  • $\begingroup$ “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? $\endgroup$
    – xzczd
    Jan 17, 2022 at 7:05
  • $\begingroup$ @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 $\endgroup$
    – Rudyard
    Jan 17, 2022 at 14:37
  • 2
    $\begingroup$ 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 ) $\endgroup$
    – xzczd
    Jan 17, 2022 at 14:47
  • $\begingroup$ 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? $\endgroup$
    – xzczd
    Jan 17, 2022 at 15:01

3 Answers 3


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] 

enter image description here


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]

enter image description here

It throws a warning:

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


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.