2
$\begingroup$

I'm trying to solve the spherically symmetric wave equation $$0 = (\partial_t^2 - \partial_r^2 + 1)\phi(t,r)\,,$$ where $\phi(t,0) = 0$. Without doing anything fancy, we can solve this equation "in a vacuum" by simply putting it in a large enough box. For example

rmax = 40;
ti = 0;
tf = 40;
ω = 1/2;
θfn = 
  NDSolveValue[{ 
    Derivative[0, 2][θ][r, t] - 
      Derivative[2, 0][θ][r, t] + θ[r, t] == 
     0, θ[0, t] == 0, 
    Derivative[1, 0][θ][rmax, t] == 0, θ[rmax, t] == 
     0, θ[r, ti] == 0, 
    Derivative[0, 1][θ][r, ti] == 
     r  A  Sech[r ω]}, θ, {r, 0, rmax}, {t, ti, tf}, 
   MaxStepSize -> 1/10];
Manipulate[
 Plot[θfn[r, t]/r, {r, 0, rmax}, PlotRange -> {-10, 10}], {t, 
  0, tf}]

However, if I want to run this code for a very long time, then the size of the grid grows like $t_f^2$. It would be desirable if I could put absorbing boundary conditions in this box, however, I run into trouble with the grid refinement. For example, the following code runs into trouble

A = 8;
rmax = 20;
ti = 0;
tf = 80;
ω = 1/2;
θfn = 
  NDSolveValue[{Derivative[0, 2][θ][r, t] - 
      Derivative[2, 0][θ][r, t] + θ[r, t] + 
      NeumannValue[Derivative[0, 1][θ][r, t], r == rmax] == 
     0, θ[0, t] == 0, θ[r, ti] == 0, 
    Derivative[0, 1][θ][r, ti] == 
     r  A  Sech[r ω] (1 - Erf[r/6])}, θ, {r, 0, 
    rmax}, {t, ti, tf}, MaxStepSize -> 1/100];
Manipulate[
 Plot[θfn[r, t]/r, {r, 0, rmax}, PlotRange -> {-1, 1}], {t, 0, 
  tf}]

As you can see, standing waves form, which should not exist. I suspected that this might be a problem with the absorbing boundary conditions I implemented using the NeumannValue function. However, even as we take the size of the box large enough for the reflection at the boundary to never happen, the standing waves still form. This is clearly a result of a grid that is too-coarse. However, the maximum step size is even smaller than in the simple example I gave above which does have the correct behavior. Therefore, I suspect that including the NeumannValue function is somehow causing the grid to be unrefined, leading to large errors.

Any suggestions to circumvent this would be appreciated!!!

$\endgroup$
6
  • $\begingroup$ "I'm trying to solve the spherically symmetric wave equation " Are you sure you're solving the correct equation?: Laplacian[f[r], {r, th, phi}, "Spherical"] $\endgroup$ – xzczd Aug 18 '20 at 3:06
  • $\begingroup$ Yeah, you can just replace $f\to g/r$ and then you return to the 1D wave equation. The only difference is the boundary condition at zero must now be $g(0) = 0$. $\endgroup$ – David Aug 18 '20 at 3:10
  • $\begingroup$ But Laplacian[g[r]/r, {r, th, phi}, "Spherical"] // Simplify gives g''[r]/r. $\endgroup$ – xzczd Aug 18 '20 at 3:17
  • $\begingroup$ Exactly - so the substitution gives an overall factor of r. Regularity at the origin just requires $g[0]=0$ $\endgroup$ – David Aug 18 '20 at 3:37
  • 1
    $\begingroup$ If Derivative[1, 0][θ][rmax, t] == 0, θ[rmax, t] == 0, θ[r, ti] == 0 is droped the rest works fine for me and fullfills this boundary condition too. With Your input I get warnings for inconsistency, should be scalar function of the spatial variables (twice), boundary and initial condition inconsistent. $\endgroup$ – Steffen Jaeschke Aug 18 '20 at 8:08
1
$\begingroup$

If

Derivative[1, 0][θ][rmax, t] == 0, θ[rmax, t] ==       0, θ[r, ti] == 0 

is dropped the rest works fine for me and fulfills this boundary condition too. With Your input, I get warnings for inconsistency, should be a scalar function of the spatial variables (twice), boundary and initial condition inconsistent.

From this experience, it follows easily that the rate of convergence for

Derivative[0, 1][θ][r, ti] == 
     r  A  Sech[r ω] (1 - Erf[r/6]) 

has to be enhanced.

There are several options to do that. For example with Piecewise or simply check whether to multiplicate this function with the appropriate factors stemming from the Laplacian again? This seems a study book example and therefore well designed and calculated.

The solutions from the view at the results in Manipulate have different orders of linearity in the origin and therefore the corresponding behavior for very large values. The Sech and Erf do not really converge to zero rapidly enough to compensate for all.

A cutoff for

Derivative[0, 1][θ][r, ti] ==  Piecewise[{{
     r  A  Sech[r ω] (1 - Erf[r/6]),0<=r<=rmax},{0,r>rmax}}]

For a sufficient rmax.

With[{rmax = 11.7, \[Omega] = 1/2, A = 1}, 
 Plot[Piecewise[{{r A Sech[r \[Omega]] (1 - Erf[r/6]), 
     0 <= r <= rmax}, {0, r > rmax}}], {r, 0, 15}]]

Plot of the radial disturbation boundary function

All for the first part of the question.

$\endgroup$
1
  • 1
    $\begingroup$ (-1) It's is the second sample that OP is having trouble. Once again, please read the question carefully before posting an answer. $\endgroup$ – xzczd Aug 18 '20 at 12:18

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.