We propose a new algorithm for solving a system of two nonlinear transcendental equations with two complex variables based on the Muller algorithm. The two-dimensional Muller algorithm is tested on systems of different type and is found to work comparably to Newton's method and Broyden's method in many cases. The new algorithm is particularly useful in systems featuring the Heun functions whose complexity may make the already known algorithms not efficient enough or not working at all. In those specific cases, the new algorithm gives distinctly better results than the other two methods. As an example for its application in physics, the new algorithm was used to find the quasi-normal modes (QNM) of Schwarzschild black hole described by the Regge-Wheeler equation. The numerical results obtained by our method are compared with the already published QNM frequencies and are found to coincide to a great extent with them. Also discussed are the QNM of the Kerr black hole, described by the Teukolsky Master equation.
1 Overview Solving a system of two complex-valued nonlinear transcendental equations numerically is a task with varying difficulty, depending on the non-linearity of the system, the types of functions involved and the dimension of the space determined by the system. There are many well-known iterative root-finding algorithms, but most of them are specialized and optimized to work with a narrow set of functions -for example polynomials or functions with real-valued roots. The two most heavily relied upon one-dimensional algorithms -the secant method and Newton's method (or the Newton-Raphson method, [1,2,3,4]) can work with a wide set of complex valued functions under proper conditions (see [4]), but they have their weak sides. Newton's method requires the evaluation of the function and its first derivative at each iteration. This increases the computational cost of the algorithm and it makes the algorithm unusable when the procedure evaluating the derivative of the function has numerical problems (for example see the discussion for the Heun functions below) or when derivative becomes zero or changes sign. The secant method avoids this limitation, but in the general case, it has lower order of convergence (∼ 1.618) compared to that of Newton's method (= 2) and the convergence of both of them is strongly dependent on the initial guess.
These problems of the algorithms are inherited by their multi-dimensional versions such as the generalized Newton-Raphson method ( [3]) and the generalized secant method (Broyden’s method, [5]). Although those problems can have varying severity, there are systems in which those algorithms cannot be used effectively. There are also some novel approaches (see [6], [7]), but when they rely on the same one-dimensional algorithms, they are likely to share their weaknesses as well. Clearly there is a need for new algorithms that will enlarge the class of functions we are able to work with efficiently.
A great challenge in front of root-finding algorithms in modern physics can be found in systems including the Heun functions. The Heun functions are unique particular local solutions of a second-order linear ordinary differential equation from the Heun type [8,9,10,11] which in the general case have 4 regular singular points. Two or more of those regular singularities can coalesce into an irregular singularity leading to confluent differential equations and their solutions: confluent Heun function, biconfluent Heun function, double confluent Heun function and triconfluent Heun function. The Heun functions generalize the hypergeometric function (and also include the Lame function, Mathieu function and the spheroidal wave functions [10,11]) and their wide applications in physics ( [11]) was summarized recently in [12]. From that paper and the cited therein, it is clear that the Heun functions will be encountered more and more in modern physics from quantum mechanics to astrophysics, hence we need adequate numerical algorithms able to deal with them.
The work with the Heun functions, however, is more than complicated. While there are analytical works on the Heun functions, they were largely neglected until recently and therefore the theory is far from complete. The only software package currently able to work with them is maple. Although, the routines that evaluate them in that package work well in the general case, there are some peculiarities -there are values of the parameters where the routines break down leading to infinities or to numerical errors. The situation with the derivatives of the Heun functions is even worst -in some cases, they simply do not work or their precision is lower than that of the Heun function itself. Also, in some cases there are no convenient power-series representations and then the Heun functions are evaluated in maple using numerical integration. Therefore the procedure goes slowly in the complex domain (compared to the hypergeometric function) which means that the convergence of the root-finding algorithm is essential. Despite all the challenges in the numerical work with the Heun functions, they offer many opportunities to modern physics. For example, they occur in the problem of quasi-normal modes (QNM) of rotating and non-rotating black holes. In this case, one has to solve a two-dimensional connected spectral problem with two complex equations in each of which one encounters the confluent Heun functions. The analytical theory of the confluent Heun function is more developed than that of the other types of Heun functions, but still many unknowns remain. Again, the evaluation of the derivative of the confluent Heun function is problematic in maple and the behavior of the function is hard to predict. In that situation Newton’s method cannot be used as a root-finding algorithm. Broyden’s algorithm works well in most cases, but it is slowly convergent even close to a root. It is clear, then, that we need a novel algorithm, that will offer quicker converg
This content is AI-processed based on open access ArXiv data.
