Boundary conditions are not satisfied with NDSolve

Question

I want to study the time evolution of a small perturbation around the static solution of the following Wave Equation

$$ -\partial_t^2 v(t,r) + \partial_r^2v(t,r) + \frac{2}{r}\partial_r v(t,r) = \frac{\partial V(v)}{\partial v}(t,r) $$

for some expression of the potential $V(v)$ that is written in the code below. The coordinates $t,r$ run over $[0,+\infty]$.

By definition, the static solution $\hat{v}(r)$ is time-independent and I require the following initial/boundary conditions

$$ \partial_r \hat{v}(r)|_{r=0} =0\,,\qquad \hat{v}(r \rightarrow +\infty) = 0. $$

Obviously, to perform numerical computations the limit $r\rightarrow+\infty$ is replaced by $r=M$ where $M\gg \ell$ where $\ell$ is the characteristic length of the problem; it turns out to be $\ell \sim 2$ for the static solution.

I want to perturb this solution at $t=0$ and see how it evolves with time. So, now I am interested in the time-dependent solution which satisfies

$$ v(t=0,r) = \hat{v}(r)\,\qquad \partial_t v(r,t)|_{t=0}=\delta \cdot 10^{-2}\,,\\ \partial_r v(t,r)|_{r=0}=0\,, \qquad v(r=M) = 0. $$

where $\delta\ll 1$.

• While the numerical static solution satisfies $\hat{v}'(r=0)=0$, the time-dependent solution I got does not. I don't understand why. For a specific example with $\delta = 0.001$, I get $\partial_r v(t,r) \sim -0.000701892$ irrespectively of the value of the time variable t. In particular, it looks the initial condition $v(t=0,r) = \hat{v}(r)$ is not satisfied. Is this normal?
• Moreover, I get the error , why? Are my boundary conditions really inconsistent?

This is my code.

V[v_] = (-1 + (1/8 (-9 + Sqrt[145]) - v)^2)^2 + 3 (1/8 (-9 + Sqrt[145]) - v)^3;

sol[rmax_, \[Delta]_] := Last@Last@ Last@NDSolve[{+D[v[r], {r, 2}] + 2/r D[v[r], {r, 1}] - (D[V[v], v] /. v -> v[r]) == 0, (D[v[r], r] /. r -> SetPrecision[10^-10, 100]) == 0, v[SetPrecision[10^-10, 100]] == SetPrecision[\[Delta], 100]}, v, {r, 10^-10, rmax}, WorkingPrecision -> 50, Method -> "Extrapolation"]

iTf = sol[30, 1.506400187591933106770472351];
Plot[{iTf[r]}, {r, 0, 30}, PlotRange -> All, Frame -> True]

iTfTime = v /. ParametricNDSolve[{-D[v[t, r], {t, 2}] + D[v[t, r], {r, 2}] + 2/r D[v[t, r], {r, 1}] - (D[V[v], v] /. v -> v[t, r]) == 0, v[0, r] == iTf[r], ((D[v[t, r], t]) /. t -> 0) == +\[Delta] 10^-2, (D[v[t, r], r] /. r -> 10^-10) == 0,v[t,30]==0}, v, {t, 0, 40}, {r, 10^-10, 30}, {\[Delta]}, WorkingPrecision -> MachinePrecision, Method -> {"MethodOfLines", "TemporalVariable" -> t, "SpatialDiscretization" -> {"TensorProductGrid", "MinPoints" -> 200}}, PrecisionGoal -> 15]

iTfTimeToPlot0 = iTfTime[0.001];

(*Checking boundary conditions in generic points*)
((D[iTfTimeToPlot0[t, r], t] /. t -> 0) /. r -> RandomReal[]) == +0.001 10^-2
(*Output: True*)

((D[iTfTimeToPlot0[t, r], r] /. r -> 10^-10) /. t -> RandomReal[]) == 0
(*Output: False*)

Update

I have tried adding the following method

Method -> {"MethodOfLines", 
"DifferentiateBoundaryConditions" -> {True, "ScaleFactor" -> 1}}

but still the two solutions ($v(t=0,r)$ and $\hat{v}(r)$ differs for small values of $r$, instead they should coincide given the boundary condition)

iTfTime = v /. ParametricNDSolve[{-D[v[t, r], {t, 2}] + D[v[t, r], {r, 2}] + 2/r D[v[t, r], {r, 1}] - (D[V[v], v] /. v -> v[t, r]) == 0, v[0, r] == iTf[r], ((D[v[t, r], t]) /. t -> 0) == +\[Delta] 10^-2, (D[v[t, r], r] /. r -> 10^-10) == 0, v[t, 30] == 0}, v, {t, 0, 40}, {r, 10^-10, 30}, {\[Delta]},  WorkingPrecision -> MachinePrecision, Method -> {"MethodOfLines", "DifferentiateBoundaryConditions" -> {True, "ScaleFactor" -> 1}}]

iTfTimeToPlot = iTfTime[0.001]
Plot[{iTfTimeToPlot[0, r], iTf[r]}, {r, 10^-10, 0.003}, PlotRange -> All]

(*Output: *)

Answer 1

I've been hesitant about whether I should mark this question as duplicate or not for quite a while, but finally decide to post an answer to elaborate because the relationship between OP's question and

2D Heat equation: inconsistent boundary and initial conditions

may be not immediately clear.

In short, it's all because the spatial grid is not dense enough, so the interpolation error shows up, even if you've added "MinPoints" -> 200, given the requirement for precision is so demanding. Adding "ScaleFactor" -> 1 (as discussed in detail in this post) is a relatively cheap way to relief this.

Let's increase the spatial grid points to 1500 and recheck:

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

r0 = 10^-10; rmax = 30; tmax = 40;

iTfTime = v /. 
   ParametricNDSolve[{-D[v[t, r], {t, 2}] + D[v[t, r], {r, 2}] + 
       2/r D[v[t, r], {r, 1}] - (D[V[v], v] /. v -> v[t, r]) == 0,

     v[0, r] == iTf[r],
     ((D[v[t, r], t]) /. t -> 0) == +δ 10^-2,
     (D[v[t, r], r] /. r -> r0) == 0,
     v[t, rmax] == 0}, v, {t, 0, tmax}, {r, r0, rmax}, {δ}, 
    Method -> Union[mol[1500, 2], mol[True, 1]]];

iTfTimeToPlot = iTfTime[0.001];
Plot[{iTfTimeToPlot[0, r], iTf[r]}, {r, r0, 0.003}, PlotRange -> All, 
 PlotStyle -> {Automatic, {Red, Dashed}}]

D[iTfTimeToPlot[0, r], r] /. r -> r0
(* 0.0000370957 *)

D[iTfTimeToPlot[20, r], r] /. r -> r0
(* 7.07301*10^-11 *)

Plot[D[iTfTimeToPlot[t, r], r] /. r -> r0 // Evaluate, {t, 0, tmax}, PlotRange -> All]

Manipulate[Plot[iTfTimeToPlot[t, r], {r, r0, rmax/4}, PlotRange -> {-1/2, 3/2}], {t, 0, 
  tmax/4}]

GIF made with 9.0.1 because 11.3 isn't at hand, but the result just looks the same.

The 0.0000370957 at $t=0$ is, as already mentioned, interpolation error:

iTfinter = Interpolation@Array[{#, iTf@#} &, 1500, {r0, rmax}];
iTfinter'[r0]
(* 0.0000370957281084964292812983780513760449138607 *)

BTW here's another post about interpolation error:

How to avoid this kind of numerical error caused by extreme parameters when using NDSolve?