Message Boards Message Boards

Solve a wave equation with boundary conditions using NDSolve?

Posted 5 years ago

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 enter image description here, 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*)
POSTED BY: mathPhys User
Reply to this discussion
Community posts can be styled and formatted using the Markdown syntax.
Reply Preview
Attachments
Remove
or Discard

Group Abstract Group Abstract