An Iterative High-Accuracy ADI Method for the 3D Parabolic Equation

Abstract

The alternating direction implicit parabolic equation (ADI-PE) method and the Crank–Nicolson parabolic equation (CN-PE) method have been widely used for solving the 3D parabolic equation (3D-PE) in radio wave propagation. The ADI-PE method is more computationally efficient than the CN-PE method. The accuracy of the ADI-PE method is improved by the higher-order Mitchell–Fairweather (MF)-ADI method. This paper presents an iterative high-accuracy (IHA)-ADI method for the 3D parabolic equation. A derivation of the proposed method is presented. The convergence and stability of the proposed method are estimated. Several numerical examples are considered to illustrate the advantages of the proposed method. The results of error analysis and a comparative study show that the proposed method is unconditionally stable and computationally efficient. The proposed method is more numerically accurate than the MF-ADI method.

This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

1. Introduction

The parabolic equation (PE) is one of the most widely used methods for the prediction of radio wave propagation [1–3]. It models the wave in a paraxial direction in a centered cone, and it can deal with the refraction and diffraction caused simultaneously by the inhomogeneous atmospheric structure and irregular topography [4–6]. PE is computationally efficient for long-range wave propagation [7]. It acts as a trade-off between the ray tracing (RT) method and the finite-differencetime-domain (FDTD) method. In fact, RT has limited accuracy, while FDTD has high computing resource consumption in some critical cases. Recently, PE has been used in the guided system for the accuracy calculation of the navigation parameter [8, 9]. Therefore, the numerical accuracy of PE is hard to overstate in the case of its application.

The split-step Fourier (SSF) method [10] and the finite difference approximation (FD) method are two methods that are used to solve PE. The FD method is more decent, due to the fact that it can model quite complex boundary conditions. Most of the PE studies mainly focus on the two-dimensional parabolic equation (2D-PE). Taking into account the solution of the FD method, some high-accuracy 2D-PE schemes [11–14] implemented using higher-order numerical approximations are proposed. Solving the three-dimensional parabolic equation (3D-PE) using the FD method is a more complex task. The Crank–Nicolson parabolic equation (CN-PE) is one of the popular FD methods for 3D-PE. The CN-PE method is second-order accurate and unconditionally stable in different spatial variables [1]. The Peaceman–Rachford (PR) alternating direction implicit (ADI) method [15] and the locally one-dimensional (LOD)-ADI method [16] are modifications of the CN-PE method. The PR-ADI and LOD-ADI methods for 3D-PE (ADI-PE method) are more computationally efficient than the CN-PE method. The ADI-PE is second-order accurate in the x, y, and z directions and is unconditionally stable, similar to the CN-PE. The accuracy of the ADI-PE method is improved by the Mitchell–Fairweather (MF)-ADI method [17], whose advantage is verified by error analysis and comparative study [18]. In recent years, compared with other 3D-PE methods, the MF-ADI method has proven to be computationally efficient and higher-order accurate. The author in [19] proposes an iterative method to improve the accuracy of the ADI-FDTD method. The convergence and stability analysis of this iterative method is presented in [20, 21].

In this paper, an iterative high-accuracy (IHA)-ADI method for the 3D parabolic equation is proposed. It is the first time to treat the important accuracy problem of the 3D parabolic equation using iteration. A deduction from the proposed method is presented. The convergence and stability of the proposed method are estimated. The numerical dispersion relations for the proposed method are also given. Moreover, the advantages of the proposed method are then demonstrated by a comparative study and error analysis. The results illustrate that the proposed IHA-ADI method approaches the efficiency of the MF-ADI method and is more numerically accurate than the MF-ADI method.

2. Iterative High-Accuracy ADI Method for 3D-PE

2.1. CN-PE Method and MF-ADI Method

The ADI-PE method is an approximate factorization of the CN-PE method. Firstly, we consider the scalar three-dimensional wave equation in Cartesian coordinates (x, y, z) [22], given by$$\begin{array}{}\text{(1)}& \frac{{\partial }^2\psi }{\partial x^2}+\frac{{\partial }^2\psi }{\partial y^2}+\frac{{\partial }^2\psi }{\partial z^2}+k^2{n}^2\psi =0\text{,}\end{array}$$where $\psi$ is a scalar component that represents the magnetic field or electric field in horizontal or vertical polarization, n is the refractive index (equals to 1, in the air), and k is the wavenumber. Assuming the x-axis is the direction of wave propagation, one of the solutions of (1) can be expressed as$$\begin{array}{}\text{(2)}& ux,y,z=\psi x,y,z{e}^{-j{k}_{0}x}\text{.}\end{array}$$

Considering the same assumption as in [1],$$\begin{array}{}\text{(3)}& \frac{{\partial }^2\psi }{\partial x^2}\ll k_0\frac{\partial \psi }{\partial x}\text{.}\end{array}$$

Substituting (2) into (1) results in$$\begin{array}{}\text{(4)}& \frac{\partial \psi }{\partial x}=\frac{1}{2j{k}_0}\frac{{\partial }^2\psi }{\partial y^2}+\frac{{\partial }^2\psi }{\partial z^2}\text{.}\end{array}$$

Equation (4) is the standard parabolic equation. By reducing the second-order derivative into a first-order derivative along the propagation direction and using the finite difference (FD) approximations for the second-order derivative, (4) can be simplified as in [1]$$\begin{array}{}\text{(5)}& 1-\frac{r_y}{4j{k}_0}\delta y-\frac{r_z}{4j{k}_0}\delta z{u}_{i,j}^{k+1}=1+\frac{r_y}{4j{k}_0}\delta y+\frac{r_z}{4j{k}_0}\delta z{u}_{i,j}^k\text{,}\end{array}$$where i, j, and k are the ordinal numbers of discrete points along y, z, and x directions.$$\begin{array}{c}\text{(6)}& \begin{array}{c}{r}_y=\frac{∆x}{∆y^2},\\ r_z=\frac{∆x}{∆z^2},\end{array}\begin{array}{c}{\delta }_y{u}_{i,j}^k=u_{i-1,j}^k-2u_{i,j}^k+u_{i+1,j}^k,\\ {\delta }_z{u}_{i,j}^k=u_{i,j-1}^k-2u_{i,j}^k+u_{i,j+1}^k.\end{array}\end{array}$$

Note that (5) is the CN-PE method. The update matrix for the CN-PE method is a sparse matrix. The solution of this method is time-consuming for electrically large problems.

Adding the error term $r_y/4j{k}_0{r}_z/4j{k}_0\delta y\delta z$ on the left side and the right side of (5), with the fourth-order approximation for the transverse plane coordinates, (5) can be expressed as$$\begin{array}{}\text{(7)}& 1-\frac{r_y}{4j{k}_0}-\frac{1}{12}\delta y1-\frac{r_z}{4j{k}_0}-\frac{1}{12}\delta z{u}_{i,j}^{k+1}=1+\frac{r_y}{4j{k}_0}+\frac{1}{12}\delta y1+\frac{r_z}{4j{k}_0}+\frac{1}{12}\delta z{u}_{i,j}^k\text{.}\end{array}$$

Equation (7) can be solved by two splitting steps given by$$\begin{array}{}\text{(8a)}& 1-\frac{r_y}{4j{k}_0}-\frac{1}{12}\delta y{u}_{i,j}^{k+1/2}=1+\frac{r_y}{4j{k}_0}+\frac{1}{12}\delta y{u}_{i,j}^k\text{,}\text{(8b)}& 1-\frac{r_z}{4j{k}_0}-\frac{1}{12}\delta z{u}_{i,j}^{k+1}=1+\frac{r_z}{4j{k}_0}+\frac{1}{12}\delta z{u}_{i,j}^{k+1/2}\text{.}\end{array}$$

Note that (9) is the MF-ADI method [17]. It is an unconditionally stable method and can be solved with high computational efficiency.

2.2. Iterative High-Accuracy ADI Method for 3D-PE

Adding $r_y/4j{k}_0{r}_z/4j{k}_0\delta y\delta z$ on the left side and the right side of (7), we obtain$$\begin{array}{}\text{(9)}& 1-\frac{r_y}{4j{k}_0}-\frac{1}{12}\delta y1-\frac{r_z}{4j{k}_0}-\frac{1}{12}\delta z{u}^{k+1}=1+\frac{r_y}{4j{k}_0}+\frac{1}{12}\delta y1+\frac{r_z}{4j{k}_0}+\frac{1}{12}\delta z{u}^k-\frac{r_y}{4j{k}_0}\frac{r_z}{4j{k}_0}\delta y\delta z{u}^{k+1}-u^k\text{.}\end{array}$$

We consider the following:$$\begin{array}{c}\text{(10)}& Q=1-\frac{r_y}{4j{k}_0}-\frac{1}{12}\delta y1-\frac{r_z}{4j{k}_0}-\frac{1}{12}\delta z\text{,}R=1+\frac{r_y}{4j{k}_0}+\frac{1}{12}\delta y1+\frac{r_z}{4j{k}_0}+\frac{1}{12}\delta z+\frac{r_y}{4j{k}_0}\frac{r_z}{4j{k}_0}\delta y\delta z{u}^k\text{,}\\ T=-\frac{r_y}{4j{k}_0}\frac{r_z}{4j{k}_0}\delta y\delta z\text{.}\end{array}$$

Then, (13) can be rewritten as$$\begin{array}{}\text{(11)}& Q{u}^{k+1}=T{u}^{k+1}+R\text{.}\end{array}$$

(11) is a linear system that can be iteratively solved in the following form:$$\begin{array}{}\text{(12)}& Q{u}_{n+1}^{k+1}=T{u}_n^{k+1}+R\text{,}\end{array}$$where $n$ denotes the $n$th iteration.

If the initial guess, $u_0^{k+1}$, is set to $u^k$, (7) can be recovered from (12). Equation (7) is the first iterative solution of (12). The scheme in (12) is straightforward. For the $n+1$th iteration, we obtain$$\begin{array}{}\text{(13)}& 1-\frac{r_y}{4j{k}_0}-\frac{1}{12}\delta y1-\frac{r_z}{4j{k}_0}-\frac{1}{12}\delta z{u}_{n+1}^{k+1}=1+\frac{r_y}{4j{k}_0}+\frac{1}{12}\delta y1+\frac{r_z}{4j{k}_0}+\frac{1}{12}\delta z{u}^k-\frac{r_y}{4j{k}_0}\frac{r_z}{4j{k}_0}\delta y\delta z{u}_n^{k+1}-u^k\text{.}\end{array}$$

Afterwards, (13) is split into two steps:$$\begin{array}{}\text{(14a)}& 1-\frac{r_y}{4j{k}_0}-\frac{1}{12}\delta y{u}_{n+1}^{tmp}=1+\frac{r_y}{4j{k}_0}+\frac{1}{12}\delta y1+\frac{r_z}{4j{k}_0}+\frac{1}{12}\delta z{u}^k+\frac{r_y}{4j{k}_0}\frac{r_z}{4j{k}_0}\delta y\delta z{u}_n^{k+1}-u^k\text{,}\text{(14b)}& \begin{array}{c}1-\frac{r_z}{4j{k}_0}-\frac{1}{12}\delta z{u}_{n+1}^{k+1}=u_{n+1}^{tmp}\end{array}\text{.}\end{array}$$

(17) is the proposed iterative high-accuracy ADI method for the 3D parabolic equation which is depicted in Figure 1.

[figure(s) omitted; refer to PDF]

2.3. Convergence and Stability Analysis

The convergence and stability of the iterated ADI-FDTD are analyzed in [20, 21]. In this section, the convergence and stability analysis of the proposed IHA-ADI method is presented.

The proposed method converges if the eigenvalues $\mu$ of the iteration matrix$$\begin{array}{c}\text{(15)}& S={K-B}^{-1}{J-A}^{-1}AB={K-B}^{-1}{J-A}^{-1}AB{K-B}^{-1}K-B\\ ={K-B}^{-1}\mathcal{C}A\mathcal{C}BK-B\text{,}\end{array}$$satisfy $\mu \le 1$, where $J=1+1/12\delta y$, $A=r_y/4j{k}_0\delta y$, $K=1+1/12\delta z$, $B=r_z/4j{k}_0\delta z$, $\mathcal{C}A={J-A}^{-1}A$, and $\mathcal{C}B={K-B}^{-1}B$. $S$ is similar to $\mathcal{C}A\mathcal{C}B$. A matrix is skew-Hermitian, and its eigenvalues $jw$ are located on the imaginary axis. The J matrix is a real positive definite symmetric matrix, and its eigenvalues satisfy ${\lambda }^J\ge 0$. The eigenvalues $\mathcal{C}j\omega$ of $\mathcal{C}A$ then satisfy$$\begin{array}{}\text{(16)}& \mathcal{C}j\omega =\frac{j\omega }{{\lambda }_J-j\omega }=\frac{\omega }{\sqrt{{\lambda }_J+{\omega }^2}}<1\text{.}\end{array}$$

The matrix $\mathcal{C}A\mathcal{C}{A}^{\ast }={J-A}^{-1}A{A}^{\ast }{J^{\ast }-A^{\ast }}^{-1}=-{J-A}^{-1}{A}^2{J+A}^{-1}=-{J-A^2}^{-1}{A}^2=-\mathcal{C}{A}^2$ has eigenvalues $0\le -\mathcal{C}{j\omega }^2={\omega }^2/{\lambda }_J^2+{\omega }^2<1$ so that$$\begin{array}{c}\text{(17)}& {\mathcal{C}A}^2=\underset{\lambda }{\max }\lambda \mathcal{C}A\mathcal{C}{A}^{\ast }=\underset{\omega }{\max }\frac{{\omega }^2}{1+{\omega }^2}<1\text{,}\mu \le \mathcal{C}A\mathcal{C}B\le \mathcal{C}A\mathcal{C}B<1\text{.}\end{array}$$

The CN-PE and the MF-ADI methods are unconditionally stable, which is estimated via von Neumann [23] by the amplification factor $G\le 1$. The stability of the proposed IHA-ADI method can be estimated via von Neumann, amplification factor $G$ satisfies.$$\begin{array}{c}\text{(18)}& G=\frac{1-R_y-1/12{\sin }^2{\sigma }_1∆y/21-R_z-1/12}{1+R_y+1/12{\sin }^2{\sigma }_1∆y/21+R_z+1/12}\frac{{\sin }^2{\sigma }_2∆z/2-R_y{R}_z{\sin }^2{\sigma }_1∆y/2{\sin }^2{\sigma }_2∆z/2{\sin }^2{\sigma }_2∆z/2}{{\sin }^2{\sigma }_2∆z/2-R_y{R}_z{\sin }^2{\sigma }_1∆y/2{\sin }^2{\sigma }_2∆z/2}\le 1\text{,}\end{array}$$where $R_y=r_y/4j{k}_0$ and $R_z=r_z/4j{k}_0$. Consequently, the proposed method converges and is unconditionally stable.

2.4. Error Analysis

The numerical dispersion relations of the CN-PE and MF-ADI methods are presented in [17]. In this section, the numerical dispersion relations for the proposed IHA-ADI method are presented.

An exact solution of the Helmholtz wave equation is given by$$\begin{array}{}\text{(19)}& \psi x,y,z=A{e}^{-j{k}_{x}x+k_{y}y+k_{z}z}\text{,}\end{array}$$where A is a constant coefficient, and $ux,y,z$ is expressed as$$\begin{array}{}\text{(20)}& ux,y,z=A{e}^{-j{k}_{x}x+k_{y}y+k_z-k_{0}z}\text{,}\end{array}$$where $k_0={k_x}^2+{k_y}^2+{k_z}^2$.

The normalized numerical wavenumber $\overline{k_0}$ can be extracted by substituting (20) into a different PE scheme. By substituting into different PE scheme, the normalized numerical wavenumber $\overline{k_0}$ can be extracted. The numerical wavevector at the x, y, and z directions defined by the angles of the spherical coordinate system is given by$$\begin{array}{c}\text{(21)}& \overline{k_x}=\overline{k_0}\cos \phi \sin \theta \text{,}\overline{k_y}=\overline{k_0}\sin \phi \sin \theta \text{,}\\ \overline{k_z}=\overline{k_0}\cos \theta \text{.}\end{array}$$

The numerical dispersion relations for the proposed iterative high-accuracy ADI method are given by$$\begin{array}{c}\text{(22a)}& e^{j{\overline{k}}_{\mathrm{0,1}}^n\lambda 1-\cos \theta \alpha \pi /N_{\lambda }^2}=1+j\alpha -\frac{1}{3}{\sin }^2\frac{{\overline{k}}_0\lambda \sin \phi \sin \theta }{2N_{\lambda }}\frac{1+j\alpha -1/3{\sin }^2{\overline{k}}_0\lambda \cos \phi \sin \theta /2N_{\lambda }}{1-j\alpha +1/3{\sin }^2{\overline{k}}_0\lambda \cos \phi \sin \theta /2N_{\lambda }}\\ +\frac{{\alpha }^2{\sin }^2{\overline{k}}_0\lambda \sin \phi \sin \theta /2N_{\lambda }\ast {\sin }^2{\overline{k}}_0\lambda \cos \phi \sin \theta /2N_{\lambda }\ast M^n-1}{1-j\alpha +1/3{\sin }^2{\overline{k}}_0\lambda \cos \phi \sin \theta /2N_{\lambda }}\text{,}\\ \text{(22b)}& e^{j{\overline{k}}_{\mathrm{0,2}}^n\lambda 1-\cos \theta \alpha \pi /N_{\lambda }^2}=\frac{1}{1-j\alpha +1/3{\sin }^2{\overline{k}}_0\lambda \sin \phi \sin \theta /2N_{\lambda }}\text{,}\end{array}$$where $\alpha =∆x/k_0∆z^2$, $∆z=∆y$, and $N_{\lambda }=\lambda /∆z=2\pi /k_0∆z$, $M^0=1$,$$\begin{array}{c}\text{(23)}& M^n=\frac{1+j\alpha -1/3{\sin }^2{\overline{k}}_0\lambda \sin \phi \sin \theta /2N_{\lambda }1+j\alpha -1/3{\sin }^2{\overline{k}}_0\lambda \cos \phi \sin \theta /2N_{\lambda }}{1-j\alpha +1/3{\sin }^2{\overline{k}}_0\lambda \cos \phi \sin \theta /2N_{\lambda }1-j\alpha +1/3{\sin }^2{\overline{k}}_0\lambda \sin \phi \sin \theta /2N_{\lambda }}\text{,}+\frac{{\alpha }^2{\sin }^2{\overline{k}}_0\lambda \sin \phi \sin \theta /2N_{\lambda }}{1-j\alpha +1/3{\sin }^2{\overline{k}}_0\lambda \cos \phi \sin \theta /2N_{\lambda }}\\ \frac{{\sin }^2{\overline{k}}_0\lambda \cos \phi \sin \theta /2N_{\lambda }{M}^{n-1}-1}{1-j\alpha +1/3{\sin }^2{\overline{k}}_0\lambda \sin \phi \sin \theta /2N_{\lambda }}\text{.}\end{array}$$

3. Numerical Results and Discussion

The relative error, which includes both the standard dispersion error of parabolic solvers and the error introduced by the PE approximation of the wave equation, is presented in this section. The relative error of the MF-ADI and CN-PE methods is verified in [18] that the MF-ADI method is more accurate. In this section, to verify the accuracy of the proposed algorithm, the proposed IHA-ADI method is compared with the CN-PE method and MF-ADI method. The numerical relative errors are analyzed as follows:

Figure 2 illustrates the relative error of the three methods with respect to $\theta$ for one iteration, where $\alpha =5$, $N_{\lambda }=4$, and $\phi =45$. It can be seen that the relative error increases with the increase in $\theta$. The relative error differences for the three methods are higher when $\theta$ increases. Among the three methods, the relative error of the proposed IHA-ADI method is the lowest. In the same error requirement, the proposed method has a larger propagation angle than other methods.

[figure(s) omitted; refer to PDF]

Figure 3 illustrates the relative error with respect to $\alpha$ when $N_{\lambda }=4$, $\theta =6$, and $\phi =45$. It can be seen that the relative error of the proposed IHA-ADI method is lower than the CN-PE method and the MF-ADI method. The result is obvious with an appropriate $\alpha$. It can be deduced from section II that $\alpha$ is proportional to $∆x$ (the step size along wave-propagation). The last adding term of the proposed iterative method in (13) is proportional to $∆x$. After iteration, the proposed scheme can modify the solution of the MF-ADI method and brings a significant advantage with an appropriate $\alpha$.

[figure(s) omitted; refer to PDF]

Figure 4 illustrates the relative error function of the azimuth angle, where $\theta =6$, $N_{\lambda }=4$, and $\alpha =5$. It can be observed that the wave-propagation along the diagonal ($\phi \approx {45}^{°}$ or ${135}^{°}$) results in the minimum error. This result is mainly caused by the numerical dispersion error of FD approximation. It is also encountered in other methods such as MF-ADI and FDTD.

[figure(s) omitted; refer to PDF]

To better illustrate the effect of the number of grid points per wavelength $N_{\lambda }$ and the zenith angle $\theta$ on the obtained results, Figure 5 shows that the relative errors of the three methods for $\theta =2^{°}$ or $6^{°}$, $\alpha =5$, and $\phi =45$. It can be seen that the relative error of the three methods decreases with the increase in $N_{\lambda }$, and finally converges to the same error floor, which represents the approximation error of PE. It is clear that the advantage of the proposed scheme is more noticeable with a higher zenith angle. The error floor with low zenith angle (e.g. 2°) is lower than that of the high zenith angle (e.g. 6°). Moreover, the effect of the number of iterations for the proposed method is analyzed, where $N_{\lambda }=4$, $\alpha =5$, and $\phi =45$, $\theta =20$. The obtained results are shown in Figure 6. It can be seen that the relative error of the proposed method decreases with the increase in the number of iteration. The proposed method can significantly reduce the error after one iteration.

[figure(s) omitted; refer to PDF]

Finally, the advantage of the proposed IHA-ADI method is analyzed by using a rectangular waveguide with perfectly electric conducting (PEC) boundary. The cross-section of the tunnel is $60\lambda \ast 60\lambda$. A Gaussian source with an inclination angle of 0° and a half-power beamwidth of 15°is placed at the center of the y and z plane and the initial x plane with the frequency of 800 MHz ($N_{\lambda }=4$ and $\alpha =15$). The CPU computation time and memory consumption of all the methods are compared in Table 1; the proposed method approaches the efficiency of the MF-ADI method, and it significantly less time and memory consuming than the CN-PE method. Figure 7 indicates that the proposed method has better agreement with the analytical results. The proposed method is more accurate than the MF-ADI method.

Table 1

Time consumption and memory consumption.

Method	Memory (KB)	CPU time (s)
CN-PE	256.8	32.152
MF-ADI	2.0	0.988
IHA-ADI (1 iteration)	2.1	1.177
IHA-ADI (2 iteration)	2.1	1.283
IHA-ADI (3 iteration)	2.1	1.515

[figure(s) omitted; refer to PDF]

4. Conclusion

In this paper, we proposed an iterative high-accuracy ADI method for solving the 3D-PE in radio wave propagation. The proposed method is unconditionally stable and computationally efficient. The effectiveness of the proposed method is demonstrated by several numerical examples. Compared with the CN-PE and the MF-ADI methods, the proposed method achieves better results in terms of numerical accuracy.

Acknowledgments

This work was supported in part by the National Natural Science Foundation of China (Grant no. 61901532) and the Natural Science Foundation of Guangdong Province (Grant no. 2015A030312010).

References

[1] M. Levy, "Parabolic equation methods for electromagnetic wave propagation," IET, vol. 5, 2000.

[2] D. Lee, A. D. Pierce, E. C. Shang, "Parabolic equation development in the twentieth century," Journal of Computational Acoustics, vol. 8 no. 4, pp. 527-637, DOI: 10.1142/s0218396x00000388, 2000.

[3] G. Apaydin, L. Sevgi, Radio Wave Propagation and Parabolic Equation Modeling, 2017.

[4] A. E. Barrios, "A terrain parabolic equation model for propagation in the troposphere," IEEE Transactions on Antennas and Propagation, vol. 42 no. 1, pp. 90-98, DOI: 10.1109/8.272306, 1994.

[5] A. E. Barrios, "Parabolic equation modeling in horizontally inhomogeneous environments," IEEE Transactions on Antennas and Propagation, vol. 40 no. 7, pp. 791-797, DOI: 10.1109/8.155744, 1992.

[6] G. Apaydin, L. Sevgi, "A novel split-stepparabolic-equation package for surface-wave propagation prediction along multiple mixed irregular-terrain paths," IEEE Antennas and Propagation Magazine, vol. 52 no. 4, pp. 90-97, DOI: 10.1109/map.2010.5638238, 2010.

[7] M. A. Lentovich, V. A. Fock, "Solution of propagation of electromagnetic waves along the earth’s surface by the method of parabolic equations," acad.sci.ussr.j.phys, vol. 10, 1946.

[8] D. D. Wang, X. L. Xi, Y. R. Pu, J. F. Liu, L. L. Zhou, "Parabolic equation method for Loran-C ASF prediction over irregular terrain," IEEE Antennas and Wireless Propagation Letters, vol. 15, pp. 734-737, DOI: 10.1109/lawp.2015.2471079, 2016.

[9] L. Huang, X. Wu, Z. Li, Y. Lu, M. Wang, Y. Long, "A parallel FDTD/ADI-PE method for ultralarge-scale propagation modeling of ils signal analysis," IEEE Antennas and Wireless Propagation Letters, vol. 19 no. 12, pp. 2245-2249, DOI: 10.1109/lawp.2020.3029193, 2020.

[10] G. Apaydin, L. Sevgi, "Two-way propagation modeling in waveguides with three-dimensionalfinite-element and split-stepfourier-based pe approaches," IEEE Antennas and Wireless Propagation Letters, vol. 10, pp. 975-978, DOI: 10.1109/lawp.2011.2167492, 2011.

[11] Q. Guo, C. Zhou, Y. Long, "Greene approximation wide-angle parabolic equation for radio propagation," IEEE Transactions on Antennas and Propagation, vol. 65 no. 11, pp. 6048-6056, DOI: 10.1109/tap.2017.2748222, 2017.

[12] Q. Guo, Y. Long, "Pade second-order parabolic equation modeling for propagation over irregular terrain," IEEE Antennas and Wireless Propagation Letters, vol. 21,DOI: 10.1109/lawp.2017.2749603, 2017.

[13] L. M. Sergeyevich, "Non-local boundary conditions for split-step pad approximations of the helmholtz equation with modified refractive index," IEEE Antennas and Wireless Propagation Letters, vol. 1 no. 1–1, 2018.

[14] M. S. Lytaev, "Numerov-padé scheme for the one-way helmholtz equation in tropospheric radio-wave propagation," IEEE Antennas and Wireless Propagation Letters, vol. 19 no. 12, pp. 2167-2171, DOI: 10.1109/lawp.2020.3026626, 2020.

[15] R. Martelly, R. Janaswamy, "An ADI-PE approach for modeling radio transmission loss in tunnels," IEEE Transactions on Antennas and Propagation, vol. 57 no. 6, pp. 1759-1770, DOI: 10.1109/tap.2009.2019891, 2009.

[16] Z. He, H. Zeng, R. S. Chen, "Two-Way propagation modeling of expressway with vehicles by using the 3-D ADI-PE method," IEEE Transactions on Antennas and Propagation, vol. 66 no. 4, pp. 2156-2160, DOI: 10.1109/tap.2018.2805018, 2018.

[17] X. Zhang, C. D. Sarris, "A high-accuracy ADI scheme for the vector parabolic equation applied to the modeling of wave propagation in tunnels," IEEE Antennas and Wireless Propagation Letters, vol. 13, pp. 650-653, DOI: 10.1109/lawp.2014.2313737, 2014.

[18] X. Zhang, C. D. Sarris, "Error analysis and comparative study of numerical methods for the parabolic equation applied to tunnel propagation modeling," IEEE Transactions on Antennas and Propagation, vol. 63 no. 7, pp. 3025-3034, DOI: 10.1109/tap.2015.2421974, 2015.

[19] S. Wang, F. L. Teixeira, C. Ji, "An iterative ADI-FDTD with reduced splitting error," IEEE Microwave and Wireless Components Letters, vol. 15 no. 2, pp. 92-94, DOI: 10.1109/lmwc.2004.842835, 2005.

[20] B. Welfert, "Analysis of iterated ADI-FDTD schemes for Maxwell curl equations," Journal of Computational Physics, vol. 222 no. 1,DOI: 10.1016/j.jcp.2006.05.038, 2007.

[21] K. Y. Jung, F. L. Teixeira, "An iterative unconditionally stable LOD–FDTD method," IEEE Microwave and Wireless Components Letters, vol. 18 no. 2, pp. 76-78, DOI: 10.1109/lmwc.2007.915026, 2008.

[22] F. B. Jensen, W. A. Kuperman, M. B. Porter, H. Schmidt, S. McKay, "Computational Ocean acoustics," Computers in Physics, vol. 9 no. 1,DOI: 10.1063/1.4823373, 1995.

[23] J. C. Strikwerda, Finite Difference Schemes and Partial Differential Equations, 2004.