# Boundary conditions are not satisfied with NDSolve

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: *)

• 1. The b.c. at $r=M$ isn't included in the code. Though NDSolve doesn't spit out bcart, please notice this will still cause severe problem: mathematica.stackexchange.com/questions/141364/… 2. Have you read this post?: mathematica.stackexchange.com/a/127411/1871 Commented Apr 19, 2019 at 12:24
• Isn't included because I think it is redundant given the dynamics of the equation. The static solution already satisfies it and I expect this to hold. However, I have tried adding this boundary condition and nothing changes. Commented Apr 19, 2019 at 12:27
• @xzczd I have seen the problem is solved by adding AccuracyGoal -> 6 when solving the time-dependent solution. Now, it gives me the solution but it returns the message Warning: scaled local spatial error estimate is too large. However, the solution looks nice, but its merely an accident Commented Apr 19, 2019 at 12:28
• @xzczd How could I have understood the value of AccuracyGoal in order to catch the solution?? This is a mystery..... Commented Apr 19, 2019 at 12:51
• 1. The missing of b.c. at $r=M$ doesn't seem to cause problem in this case, but it's no more than an accident, please add all the necessary b.c. to the code. 2. Please read the post linked above about ibcinc warning. 3. AFAIK, struggling with AccuracyGoal and PrecisionGoal is almost always the wrong way to go when dealing with PDE, for more info please read this post: mathematica.stackexchange.com/q/118249/1871 Commented Apr 19, 2019 at 12:55

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?

• Thank you. However, the solution behaves correctly at values $r> 1$, while for small values $r \sim 0$ but still $r>10^{-10}$ the solution oscillates, look at the red square here ibb.co/M1FDNdj. This is not good.... this problem is related to the interpolation error, right? I am going to read the link you referred to. Commented Apr 24, 2019 at 12:45
• @apt45 Can't reproduce in v11.3, which version are you in? What code did you use to generate the plot? Commented Apr 24, 2019 at 12:59
• I also use v11.3. Here you can download the notebook dropbox.com/s/gal9yav4caw06ec/example.nb?dl=0 Commented Apr 24, 2019 at 13:26
• Try this Plot[iTfTimeToPlot[3.6, r], {r, 10^-10, 1}, PlotRange -> All] I get this plot ibb.co/DzQBYY1 Commented Apr 24, 2019 at 13:28
• @apt45 1. Please avoid sharing code with external download links. 2. Can't reproduce the result at t==3.6, either. 3. Are you facing the problem mentioned here?: mathematica.stackexchange.com/questions/191358/… Commented Apr 24, 2019 at 13:37