The parabolic equation (PE) is an effective method for analyzing the propagation of electromagnetic waves, sound waves, and light waves in large areas situated in complex environments. Common PE analysis techniques include split-step parabolic equation (SSPE) and finite element method-based parabolic equation (FEM-PE) [1]. SSPE, which is based on the fast Fourier transform (FFT), offers the advantage of fast analysis but is accompanied by the disadvantage of complicated boundary condition processing. In contrast, FEM-PE is a matrix-based analysis method that has a more complex formulation process than SSPE but can easily handle boundary conditions.

In conventional PE based on the extended analysis region, Hamming or Hann windows are employed as the transparent boundary condition (TBC) to suppress non-physical reflections of radio waves at the top boundary, resulting in an increase in computational burden [2].

In this letter, we propose a novel TBC algorithm based on the matrix pencil method (MPM) that suppresses non-physical reflections in FEM-PE analysis and does not require additional analysis region for TBC processing. To verify the feasibility of the proposed algorithm, its analytical results and computational complexity were compared with those of a conventional SSPE.

Fig. 1 presents the two-dimensional (2D) projected propagation region. Notably, we adopted the time convention e^−iωt throughout the analysis. In 2D radio wave propagation analysis, when dealing with very large areas, analysis is conducted using PE, as shown in Eq. (1), instead of the Helmholtz equation:

where k₀ = 2π/λ₀ denotes the free-space wave number, u is the electromagnetic field, and λ₀ is the free-space wavelength. Furthermore, n = n(x) represents the refractive index, where z and x denote the range and height coordinates, respectively. Notably, in Eq. (1), a one-dimensional FEM analysis was applied to solve the second-order differential equation for x of the PE. Along these lines, if first-order linear polynomials are considered the basis functions for each element, u can be expressed as follows:

where,

In Eq. (3), x1e and x2e are the coordinates of the first and second nodes of the e-th element.

Based on Eq. (2), Galerkin’s method and the Crank-Nicolson method can be applied to calculate the vertical field profile at the next range step as follows [3]:

where,

In Eq. (4), matrix A refers to the system matrix, which has a tridiagonal matrix structure. The boundary condition of the ground or sea surface, assumed to be the perfect electric conductor (PEC) in this letter, was implemented at x = 0 and the TBC was applied at x = X.

To apply the TBC at x = X, the first-order radiation boundary equation was applied as follows:

In Eq. (9), the value of k_x varies for each range step due to the influence of the beam direction and the atmospheric refractive index. Notably, in this letter, the value for k_x estimated at the l-th range step was used to apply Eq. (9) to the (l + 1)-th range step. To estimate k_x, the MPM was applied to the field u at the l-th range step. The process for estimating k_x based on the MPM is presented below.

In the l-th range step, a Hankel matrix was constructed by sampling S numbers of u from the (N_x – 1) node on the x-axis, with S having a considerably smaller value than N_x, as follows:

Singular value decomposition was applied to Eq. (10) to find that H_S = U∑V_H. Using the orthonormal matrix U, eigenvalues ξ_p = e^ik_x,pΔx of U1† U2 were evaluated [4], where “†” refers to the Moore-Penrose pseudo-inverse symbol, and U₁ and U₂ denote the matrices obtained by deleting the last and first rows of the U matrix, respectively. Ultimately, the required k_x,p in Eq. (9) was calculated as kx,p=ln ξpiΔx [4]. Since only the dominant k_x,p was required in Eq. (9), only case k_x,1, with p = 1, was considered. Therefore, using k_x,1 at the (l + 1)-th range step, the central difference method was applied to Eq. (9) as follows:

Finally, Eq. (11) was applied as the TBC in Eq. (4), where matrix A refers to a tridiagonal band matrix with a bandwidth of 3. The direct method achieved a computational complexity of O(9N_x) at each range step. In contrast, the computational complexity of an FFT-based SSPE is O(8N_xlog(4N_x), indicating that the proposed FEM-PE exhibits better computational performance [5].

To validate the proposed FEM-PE, a comparison of its results with those of a conventional SSPE with a Hamming window was conducted. The analysis region was set considering a range of 0–100 km and a height of 0–1 km, with 1001 nodes along both the range and height axes. At the initial range of z = 0, the vertical profile of u¹ was assumed to follow a Gaussian pattern with a half power beam width of 1° and elevation angle of −1°. The source frequency was 300 MHz, and the altitude was 600 m, which was set to horizontal polarization. The ground surface was assumed to be an irregular PEC terrain with a maximum height of 150 m. Notably, atmospheric refractivity n(x) was replaced by modified refractivity M(x) to account for the radio horizon effect. The modified refractivity was bilinear, with M(0) = 360 M-units and dM(x)/dx being −200 M-units/km up to an altitude of 300 m and +400 M-units/km above 300 m, resulting in surface ducting. Furthermore, S = 6 was considered for estimating k_x,1, with the matrix H in Eq. (10) having a size of 4 × 3. Consequently, the computational load of the MPM process was negligible.

Fig. 2(a) and 2(b) depict the results of the SSPE analysis and the proposed FEM-PE method. The results of the proposed algorithm are similar to those attained by Hamming window-based SSPE, showing the suppression of non-physical reflected waves at the upper boundary. Furthermore, the normalized root-mean-square error (RMSE) of the two methods was 0.018, indicating a high degree of agreement. Employing a system equipped with a 13th Gen Intel Core i7-13700KF CPU and 128 GB RAM, it was observed that the SSPE method took 0.155 seconds to complete, whereas the proposed method required only 0.072 seconds, indicating that the latter was approximately 2.15 times faster than the former. This further emphasizes that although the proposed FEM-PE necessitates MPM calculation at each range step, it has a negligible effect on the amount of calculation required.

In this letter, we propose a novel TBC using the MPM algorithm as an alternative to Hamming or Hann window techniques for suppressing non-physical reflections at the upper boundary in PE analysis. The proposed FEM-PE with TBC/MPM effectively suppressed non-physical reflections at the upper boundary and demonstrated superior computational efficiency owing to the characteristics of its tridiagonal banded matrix. Therefore, the proposed FEM-PE can be used in various PE analysis applications involving complex media and boundary conditions.