search

BlackHolePerturbationToolkit / Teukolsky

A Mathematica package for computing solutions to the Teukolsky equation.
http://bhptoolkit.org/Teukolsky/
MIT License
19 stars 3 forks source link

FindRoot Error #51

Closed susannabarsanti closed 4 months ago

susannabarsanti commented 4 months ago

For specific input parameters, the TeukolskyPointParticleMode function doesn't work. It returns a "FindRoot" error and never completes the computation.

Three examples are: 

orbit = KerrGeoOrbit[0.9`30, 10, 0, 1];
TeukolskyPointParticleMode[0, 8, 6, 0, 0, orbit];
orbit = KerrGeoOrbit[0.9`300, 10, 0, 1];
TeukolskyPointParticleMode[0, 8, 6, 0, 0, orbit];
orbit = KerrGeoOrbit[0.9`300,4.4775`300, 0, 1];
TeukolskyPointParticleMode[-2, 8, 2, 0, 0, orbit];

The error message is " FindRoot::cvmit: Failed to converge to the requested accuracy or precision within 100 iterations, " which appears twice.

barrywardell commented 4 months ago

There appear to be two things going on here:

  1. The messages are generated by the SpinWeightedSpheroidalHarmonics package when it tries to compute the eigenvalue. As far as I can see the eigenvalue is computed fairly accurately so the errors can probably be ignored (the main place high precision on the eigenvalue may be required is in computing the renormalized angular momentum, but that doesn't seem to be the case here). Commit 607d306 of SpinWeightedSpheroidalHarmonics adds an extra Message to make clear where the FindRoot message is coming from.
  2. It takes a long time/never computes the Teukolsky solution. The calculation is getting stuck while computing the spheroidal harmonic, line 134 of TeukolskyMode.m. Commit 3c73799 of SpinWeightedSpheroidalHarmonics adds a MaxTerms option which will prevent the calculation from getting stuck and instead produce a warning message.

With the two fixes above, your test cases now run fine for me. However, unless you really need high accuracy an even more efficient solution in this case is to just work at machine precision:

orbit = KerrGeoOrbit[0.9`30, 10, 0, 1];
TeukolskyPointParticleMode[0, 8, 6, 0, 0, orbit];

produces an answer that is accurate to about 8 digits in the latest release (Teukolsky 1.0.4 and SpinWeightedSpheroidalHarmonics 1.0.0) without requiring either of the above fixes (previously it was necessary to work at high precision, but recent versions of the Teukolsky package have improved a lot so that it's not necessary most of the time).