# Solve a wave equation with boundary conditions using NDSolve?

Posted 1 year ago
642 Views
|
0 Replies
|
0 Total Likes
|
 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, 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*)