Two-dimensional generalization of the Muller root-finding algorithm and its applications

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.

💡 Research Summary

**
The paper introduces a novel numerical method for solving systems of two nonlinear transcendental equations in two complex variables by extending the classic one‑dimensional Müller root‑finding algorithm to a two‑dimensional setting. The authors begin by reviewing the limitations of widely used methods: Newton’s method requires the explicit Jacobian and can fail to converge if the initial guess is poor; Broyden’s quasi‑Newton scheme avoids explicit differentiation but still relies on a linear approximation of the system and may suffer from slow convergence or divergence when the underlying functions are highly non‑linear or poorly conditioned.

The proposed “2‑D Müller” algorithm constructs a quadratic approximation of the vector‑valued function using three distinct points in the complex (z₁, z₂) space. By fitting a second‑order polynomial surface to the function values at these points, the method obtains an analytic expression for the complex roots of the quadratic surrogate. These roots become the next iterates, and the process repeats until a prescribed tolerance on the function norm is achieved. Crucially, the algorithm never computes derivatives; it only needs function evaluations, which makes it especially attractive for problems involving special functions whose derivatives are either unavailable or numerically unstable.

A rigorous convergence analysis is provided. The authors show that, provided the discriminant of the quadratic surrogate does not vanish in the region of interest, the iteration defines a locally contractive map. They also discuss branch selection rules to avoid jumping between spurious complex roots, and they derive error‑propagation formulas that demonstrate quadratic convergence in the vicinity of a simple root, comparable to Newton’s method but without the Jacobian cost.

To assess performance, the authors benchmark the 2‑D Müller method against Newton and Broyden on three families of test problems: (i) simple polynomial systems, (ii) transcendental systems involving exponential, sine, and cosine functions, and (iii) systems that contain the confluent Heun function (HeunC). In the first two families all three methods converge, with the 2‑D Müller method showing comparable iteration counts and CPU times. In the Heun‑function tests, however, Newton frequently diverges because numerical differentiation of HeunC is ill‑conditioned, and Broyden’s Jacobian approximations become unstable. The 2‑D Müller method, relying solely on HeunC evaluations, converges robustly and typically requires 30 % fewer iterations than Broyden.

The most compelling application presented is the computation of quasi‑normal modes (QNMs) of black holes, a problem that naturally leads to complex eigenvalue equations involving Heun functions. For the Schwarzschild case, the Regge–Wheeler equation reduces to a confluent Heun equation; the QNM frequencies are the complex roots of a transcendental condition that couples HeunC and its derivative. Using the 2‑D Müller algorithm, the authors locate these roots directly. Their results match published QNM frequencies obtained by Leaver’s continued‑fraction method and other high‑precision techniques to within 10⁻⁶ relative error. They further extend the approach to the Kerr black hole, where the Teukolsky master equation yields a more intricate Heun‑type condition that depends on the rotation parameter a. Again, the algorithm efficiently maps out the a‑dependent QNM spectrum, demonstrating flexibility in handling multi‑parameter complex root problems.

The paper concludes with a discussion of limitations and future directions. While the quadratic surrogate provides rapid local convergence, the method can become trapped in a local basin if the initial guess is far from any root, especially when multiple roots lie close together and the discriminant approaches zero. The authors suggest multi‑start strategies, adaptive selection of the three interpolation points, and possible extensions to higher‑order Müller schemes (cubic, quartic) to improve robustness. They also note that their current implementation is CPU‑based (C++/Fortran) and that GPU acceleration could enable large‑scale parameter scans, such as systematic QNM surveys across a wide range of black‑hole spins and charges.

In summary, the work delivers a practical, derivative‑free root‑finding technique that is particularly well‑suited for complex systems involving special functions like HeunC. By demonstrating comparable or superior performance to Newton and Broyden on challenging test cases and by successfully applying the method to physically important black‑hole QNM calculations, the authors provide a valuable addition to the numerical toolbox of both applied mathematicians and theoretical physicists.