NDSolve: method of lines: unexpected error: insufficient boundary conditions

Question

As a test of my ability to master the Method of Lines I tried to NDSolve the PDE's system $$\partial_t \varphi = \varpi\qquad\partial_t\varpi=\frac{1}{r}\partial^2_r(r\varphi)-\varphi$$ in the area $$(t,r)\in[0,1]\times[0,1] $$ along with initial-boundary conditions $$ \varphi(0,r)=\cos(\pi r)\qquad \varpi(0,r)=\sin(\pi r)\qquad \varpi(t,1)=0 $$

The problem is mathematically well posed:

1.Boundary condition $\varpi(t,1)=0$ implies $\varphi(t,1)=const.$

2.Initial data at $r=0$ should evolve in time just as at any other point inside the sphere of radius $R=1$.

However my code

    rmax = 1; tmax = 1; e = 10^-21;
s = NDSolve[{D[φ[t, r], t] == p[t, r], 
    D[p[t, r], t] == D[r*φ[t, r], {r, 2}]/r - φ[t, r], 
    p[0, r] == Sin[ π*r/rmax], φ[0, r] == Cos[ π*r/rmax], 
    p[t, rmax] == 0}, {p, φ}, {t, 0, tmax}, {r, e, rmax}, 
   Method -> {"MethodOfLines", 
     "SpatialDiscretization" -> {"TensorProductGrid", 
       "MinPoints" -> 80, "MaxPoints" -> 100, 
       "DifferenceOrder" -> "Pseudospectral"}}];

results in the following error:

NDSolve::bcart: Warning: an insufficient number of boundary conditions have been specified for the direction of independent variable r. Artificial boundary effects may be present in the solution

How can I make Method of lines and Mathematica understand that boundary conditions are just fine?

PS: The equations above are the Klein-Gordon equations for scalar field of mass $m=1$ in hamiltonian form and assuming spherical symmetry.

Answer 1

1. You don't need to introduce the intermediate variable $ \varpi$, NDSolve is clever enough to handle 2nd order derivative respects to $t$.

2. Though not necessary when analytically solving the problem, the boundary condition (b.c.) at $r=0$ is necessary for numeric solver, that's the reason why bcart warning pops up.

3. Don't trust the result when bcart pops up, for more information check this post.

4. Strictly speaking, there's an implicit b.c. called finiteness condition for this problem i.e. the solution is bounded in the domain, and the b.c. at $r=0$ can be deduced based on this condition.

OK, now let's fix your code. First simplify the code a little:

(* const is chosen casually because it's not given in the question. *)
const = -1;
With[{φ = φ[t, r]},
 eq = D[φ, t, t] == D[r φ, r, r]/r - φ;
 ic = {φ == Cos[π*r/rmax], D[φ, t] == Sin[π*r/rmax]} /. 
   t -> 0;
 bc@1 = φ == const /. r -> rmax;]

Then let's deduce the b.c. at $r=0$:

Clear@e;
bc@2 = Numerator@Together[Subtract @@ eq] == 0 /. r -> 0 /. h_[t, 0] :> h[t, e] // 
  Simplify
(* Derivative[0, 1][φ][t, e] == 0 *)

The last step is solving the problem, and we have 2 possible choices for spatial discretization. One is the default "TensorProductGrid" method:

mol[n:_Integer|{_Integer..}, o_:"Pseudospectral"] := {"MethodOfLines", 
  "SpatialDiscretization" -> {"TensorProductGrid", "MaxPoints" -> n, 
    "MinPoints" -> n, "DifferenceOrder" -> o}}
mol[tf:False|True,sf_:Automatic]:={"MethodOfLines",
"DifferentiateBoundaryConditions"->{tf,"ScaleFactor"->sf}}

rmax = 1; tmax = 1; e = 10^-21;

sol = NDSolveValue[{eq, ic, bc@1, bc@2}, φ, {t, 0, tmax}, {r, e, rmax}, 
  Method -> Union[mol[100, 4], mol[True, 1]]]

Actually mol[True, 1] happens to be unnecessary for your specific b.c., but this option is often necessary when ibcinc warning pops up, for more information please check this post.

The other is the "FiniteElement" method:

solfem = NDSolveValue[{eq, ic, bc@1}, φ, {t, 0, tmax}, {r, 0, rmax}, 
  Method -> {"MethodOfLines", 
    "SpatialDiscretization" -> {"FiniteElement", 
      "MeshOptions" -> MaxCellMeasure -> 0.01}}]

Notice "FiniteElement" shows several advantages compared to "TensorProductGrid" in this case:

1. We can use $r=0$ as the left boundary because "FiniteElement" is able to handle the removable singularity there.

2. bc@2 can be omitted because 0 Neumann value is the default setting of FiniteElement. For more information please check the Details of document of NeumannValue.

3. b.c. is imposed in a different way so we don't need to deal with the "DifferentiateBoundaryConditions" option.

Finally an illustration of the solution:

Plot3D[solfem[t, r], {t, 0, tmax}, {r, e, rmax}]