Stable rational approximations for parabolic equation methods

TL;DR

The paper addresses the challenge of stable, accurate rational approximations for parabolic equation operators in wave propagation, particularly in fluid-elastic waveguides where evanescent modes must be suppressed without sacrificing propagating modes. It applies the adaptive AAA algorithm to obtain barycentric rational approximations of operators like $f(q)=\sqrt{1+q}$ and $f(q)=\exp\{i\sigma(-1+\sqrt{1+q})\}$, yielding spectra-mapping that remains stable across both evanescent and propagating regions and enabling efficient split-step Padé marching. Compared with standard approaches (constraint equations and rotated-branch-cut Padé), AAA demonstrates comparable or superior accuracy to reference solutions and robust stability without parameter tuning, even in challenging geometries with thin elastic overlays. The findings support broader, more reliable PE-based simulations in complex fluid-elastic media and suggest AAA as a practical tool for achieving large-range steps and improved computational efficiency in parabolic-wave models.

Abstract

Modern parabolic equation (PE) methods for wave propagation rely on application of a variety of fractional-powered differential operators. Rational approximations of these operators need to properly map their spectra onto the complex plane, accurately handling propagating modes while annihilating evanescent ones. Standard approaches for stable and accurate rational approximations include rotating the branch cut of the operators or imposing stability constraint equations, and have yielded accurate results for wave propagation in a variety of fluid, elastic, and fluid-elastic waveguides. The stability constraint method, however, does not yield operators that are stable for all fluid-elastic waveguides, and a recent study of waveguides comprised of a thin elastic layer overlaying a thick fluid layer revealed instabilities in the approximations derived from rotated operators. In this paper, we demonstrate the applicability of a different rational approximation method, the recently-developed adaptive Antoulas-Anderson (AAA) algorithm, to simulations of wave propagation using the fluid-elastic parabolic equation. We find that simulations using operators approximated using the AAA algorithm provide excellent agreement with reference solutions, with errors in transmission loss comparable to, and often less than, that of simulations using the rotated operator method. In addition, we find that the AAA algorithm allows for the application of the split-step Padé method to fluid-elastic waveguides, which yields a large gain in computational efficiency.

Figures (19)

• Figure 1: The mapping of the real line segment $(-10,10)$ onto the complex plane by the square root operator $\sqrt{1+q}$. Shown (on both panels) are the the analytic result (solid black), and approximations using (a) the stability constraint method with a perturbation off the real line(dash-dash-dotted orange) and (b) the AAA algorithm (dashed red), the ROP method with $\theta=80^\circ$ (dashed-dot-dotted purple), $\theta=90^\circ$ (dash-dotted green), $\theta=100^\circ$ (dotted blue). All approximated operators are of order $N=12$.
• Figure 2: The mapping of the real line segment onto the complex plane by the square root operator $\sqrt{1+q}$ zoomed in around the origin. Shown are the analytic result (solid black), and approximations using the AAA algorithm (dashed red), the ROP method with $\theta=80^\circ$ (dashed-dot-dotted purple), $\theta=90^\circ$ (dash-dotted green), $\theta=100^\circ$ (dotted blue). All approximated operators are of order $N=12$.
• Figure 3: The error of the approximations $g(q)$ relative to the analytic function $f(q)=(1+q)^{1/2}$ for $q$ in the range $(-10,10)$. Shown on are the approximations using the stability constraint method (dash-dash-dotted orange), the AAA algorithm (dashed red), and the ROP method with $\theta=80^\circ$ (dashed-dot-dotted purple), $\theta=90^\circ$ (dash-dotted green), $\theta=100^\circ$ (dotted blue). All approximated operators are of order $N=12$.
• Figure 4: The mapping of the real line segment $(-10,10)$ onto the complex plane by the exponential operator $\exp\{i\sigma(-1+\sqrt{1+q})\}$ with (a),(c),(e) $\sigma=5$ and (b),(d),(f) $\sigma=10$. Shown (on all panels) are the the analytic result (solid black), and approximations using (a),(b) the stability constraint method with $g(-3) = 0$ (dash-dash-dotted orange), and (c),(d) the ROP method with $\theta=20^\circ$ (dash-dot-dotted purple), $\theta=40^\circ$ (dash-dotted green), and $\theta=60^\circ$ (dotted blue), and (e),(f) the AAA algorithm (dashed red). All approximated operators are of order $N=12$.
• Figure 5: The error of the approximations $g(q)$ relative to the analytic function $f(q)=\exp\{i\sigma(-1+\sqrt{1+q})\}$ with (a) $\sigma=5$ and (b) $\sigma=10$, for $q$ in the range $(-10,10)$. Shown on both panels are the approximations using the stability constraint method with $g(-3) = 0$ (dash-dash-dotted orange), the AAA algorithm (dashed red), and the ROP method with $\theta=20^\circ$ (dashed-dot-dotted purple), $\theta=40^\circ$ (dash-dotted green), $\theta=60^\circ$ (dotted blue). All approximated operators are of order $N=12$.