# 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 – xzczd Apr 19 '19 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. – apt45 Apr 19 '19 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 – apt45 Apr 19 '19 at 12:28
• @xzczd How could I have understood the value of AccuracyGoal in order to catch the solution?? This is a mystery..... – apt45 Apr 19 '19 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 – xzczd Apr 19 '19 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. – apt45 Apr 24 '19 at 12:45
• @apt45 Can't reproduce in v11.3, which version are you in? What code did you use to generate the plot? – xzczd Apr 24 '19 at 12:59
• I also use v11.3. Here you can download the notebook dropbox.com/s/gal9yav4caw06ec/example.nb?dl=0 – apt45 Apr 24 '19 at 13:26
• Try this Plot[iTfTimeToPlot[3.6, r], {r, 10^-10, 1}, PlotRange -> All] I get this plot ibb.co/DzQBYY1 – apt45 Apr 24 '19 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/… – xzczd Apr 24 '19 at 13:37