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.

  • $\begingroup$ It produces a solution for me (the warning is just a warning). Is the solution acceptable? $\endgroup$
    – Michael E2
    Dec 18, 2018 at 16:25
  • $\begingroup$ I was suggesting you could check... $\endgroup$
    – Michael E2
    Dec 18, 2018 at 17:56
  • $\begingroup$ The implication `[CurlyPi][t, 1] == 0 [DoubleRightArrow] [CurlyPhi][t, 1] == const``is wrong! $\endgroup$ Dec 18, 2018 at 19:03
  • $\begingroup$ $\varpi(t,1)=0\Rightarrow\partial_t\varphi(t,1)=0\Rightarrow \varphi(t,1)=\varphi(0,1)=\cos(\pi)=-1$ $\endgroup$
    – dkstack
    Dec 19, 2018 at 7:50
  • $\begingroup$ The problem is that "artificial boundary conditions" may modify the problem and then offer an acceptable solution but for the wrong problem. $\endgroup$
    – dkstack
    Dec 19, 2018 at 10:24

1 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$:

bc@2 = Numerator@Together[Subtract @@ eq] == 0 /. r -> 0 /. h_[t, 0] :> h[t, e] // 
(* 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}}

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}]

Mathematica graphics

  • $\begingroup$ Use of intermediate variable $\varpi$ is a matter of technique, not necessity. The technical problem is that Mathematica should understand that boundary condition for $\varpi$ implies a boundary condition $\varphi$. So something is wrong with my code. $\endgroup$
    – dkstack
    Dec 18, 2018 at 17:24
  • 1
    $\begingroup$ @dkstack Oops, I made a mistake in translating the b.c. at r == rmax. Now it's corrected. Notice the key point is b.c. at $r=0$ needs to be given explicitly when TensorProductGrid is chosen. $\endgroup$
    – xzczd
    Dec 19, 2018 at 8:05

Your Answer

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