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

Authors: Xiaoping Wu; Zihao Li; Zhixi Liang; Yunliang Long
Journal: International Journal of Antennas and Propagation
 (2023)

---

## 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.


---

## Body

## 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(1)∂2ψ∂x2+∂2ψ∂y2+∂2ψ∂z2+k2n2ψ=0,where ψ 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(2)ux,y,z=ψx,y,ze−jk0x.Considering the same assumption as in [1],(3)∂2ψ∂x2≪k0∂ψ∂x.Substituting (2) into (1) results in(4)∂ψ∂x=12jk0∂2ψ∂y2+∂2ψ∂z2.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](5)1−ry4jk0δy−rz4jk0δzui,jk+1=1+ry4jk0δy+rz4jk0δzui,jk,where i, j, and k are the ordinal numbers of discrete points along y, z, and x directions.(6)ry=∆x∆y2,rz=∆x∆z2,δyui,jk=ui−1,jk−2ui,jk+ui+1,jk,δzui,jk=ui,j−1k−2ui,jk+ui,j+1k.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 termry/4jk0rz/4jk0δyδ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(7)1−ry4jk0−112δy1−rz4jk0−112δzui,jk+1=1+ry4jk0+112δy1+rz4jk0+112δzui,jk.Equation (7) can be solved by two splitting steps given by(8a)1−ry4jk0−112δyui,jk+1/2=1+ry4jk0+112δyui,jk,(8b)1−rz4jk0−112δzui,jk+1=1+rz4jk0+112δzui,jk+1/2.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

Addingry/4jk0rz/4jk0δyδz on the left side and the right side of (7), we obtain(9)1−ry4jk0−112δy1−rz4jk0−112δzuk+1=1+ry4jk0+112δy1+rz4jk0+112δzuk−ry4jk0rz4jk0δyδzuk+1−uk.We consider the following:(10)Q=1−ry4jk0−112δy1−rz4jk0−112δz,R=1+ry4jk0+112δy1+rz4jk0+112δz+ry4jk0rz4jk0δyδzuk,T=−ry4jk0rz4jk0δyδz.Then, (13) can be rewritten as(11)Quk+1=Tuk+1+R.(11) is a linear system that can be iteratively solved in the following form:(12)Qun+1k+1=Tunk+1+R,where n denotes the nth iteration.If the initial guess,u0k+1, is set to uk, (7) can be recovered from (12). Equation (7) is the first iterative solution of (12). The scheme in (12) is straightforward. For the n+1th iteration, we obtain(13)1−ry4jk0−112δy1−rz4jk0−112δzun+1k+1=1+ry4jk0+112δy1+rz4jk0+112δzuk−ry4jk0rz4jk0δyδzunk+1−uk.Afterwards, (13) is split into two steps:(14a)1−ry4jk0−112δyun+1tmp=1+ry4jk0+112δy1+rz4jk0+112δzuk+ry4jk0rz4jk0δyδzunk+1−uk,(14b)1−rz4jk0−112δzun+1k+1=un+1tmp.(17) is the proposed iterative high-accuracy ADI method for the 3D parabolic equation which is depicted in Figure 1.Figure 1
The proposed iterative high-accuracy ADI method for 3D-PE. The green circular dots represent the elements that are updated in (14a) at a time. The blue circular dots represent the elements that are updated in (14b) at a time. Step 1 is explicit in the z direction and implicit in the y direction. Step 2 is explicit in the y direction and implicit in the z direction.


### 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μ of the iteration matrix(15)S=K−B−1J−A−1AB=K−B−1J−A−1ABK−B−1K−B=K−B−1CACBK−B,satisfy μ≤1, where J=1+1/12δy, A=ry/4jk0δy, K=1+1/12δz, B=rz/4jk0δz, CA=J−A−1A, and CB=K−B−1B. S is similar to CACB. 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 λJ≥0. The eigenvalues Cjω of CA then satisfy(16)Cjω=jωλJ−jω=ωλJ+ω2<1.The matrixCACA∗=J−A−1AA∗J∗−A∗−1=−J−A−1A2J+A−1=−J−A2−1A2=-CA2 has eigenvalues 0≤−Cjω2=ω2/λJ2+ω2<1 so that(17)CA2=maxλλCACA∗=maxωω21+ω2<1,μ≤CACB≤CACB<1.The CN-PE and the MF-ADI methods are unconditionally stable, which is estimated via von Neumann [23] by the amplification factor G≤1. The stability of the proposed IHA-ADI method can be estimated via von Neumann, amplification factor G satisfies.(18)G=1−Ry−1/12sin2σ1∆y/21−Rz−1/121+Ry+1/12sin2σ1∆y/21+Rz+1/12sin2σ2∆z/2−RyRzsin2σ1∆y/2sin2σ2∆z/2sin2σ2∆z/2sin2σ2∆z/2−RyRzsin2σ1∆y/2sin2σ2∆z/2≤1,where Ry=ry/4jk0 and Rz=rz/4jk0. 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(19)ψx,y,z=Ae−jkxx+kyy+kzz,where A is a constant coefficient, and ux,y,z is expressed as(20)ux,y,z=Ae−jkxx+kyy+kz−k0z,where k0=kx2+ky2+kz2.The normalized numerical wavenumberk0¯ can be extracted by substituting (20) into a different PE scheme. By substituting into different PE scheme, the normalized numerical wavenumber k0¯ 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(21)kx¯=k0¯cosφsinθ,ky¯=k0¯sinφsinθ,kz¯=k0¯cosθ.The numerical dispersion relations for the proposed iterative high-accuracy ADI method are given by(22a)ejk¯0,1nλ1−cosθαπ/Nλ2=1+jα−13sin2k¯0λsinφsinθ2Nλ1+jα−1/3sin2k¯0λcosφsinθ/2Nλ1−jα+1/3sin2k¯0λcosφsinθ/2Nλ+α2sin2k¯0λsinφsinθ/2Nλ∗sin2k¯0λcosφsinθ/2Nλ∗Mn−11−jα+1/3sin2k¯0λcosφsinθ/2Nλ,(22b)ejk¯0,2nλ1−cosθαπ/Nλ2=11−jα+1/3sin2k¯0λsinφsinθ/2Nλ,where α=∆x/k0∆z2, ∆z=∆y, and Nλ=λ/∆z=2π/k0∆z, M0=1,(23)Mn=1+jα−1/3sin2k¯0λsinφsinθ/2Nλ1+jα−1/3sin2k¯0λcosφsinθ/2Nλ1−jα+1/3sin2k¯0λcosφsinθ/2Nλ1−jα+1/3sin2k¯0λsinφsinθ/2Nλ,+α2sin2k¯0λsinφsinθ/2Nλ1−jα+1/3sin2k¯0λcosφsinθ/2Nλsin2k¯0λcosφsinθ/2NλMn−1−11−jα+1/3sin2k¯0λsinφsinθ/2Nλ.


## 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(1)∂2ψ∂x2+∂2ψ∂y2+∂2ψ∂z2+k2n2ψ=0,where ψ 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(2)ux,y,z=ψx,y,ze−jk0x.Considering the same assumption as in [1],(3)∂2ψ∂x2≪k0∂ψ∂x.Substituting (2) into (1) results in(4)∂ψ∂x=12jk0∂2ψ∂y2+∂2ψ∂z2.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](5)1−ry4jk0δy−rz4jk0δzui,jk+1=1+ry4jk0δy+rz4jk0δzui,jk,where i, j, and k are the ordinal numbers of discrete points along y, z, and x directions.(6)ry=∆x∆y2,rz=∆x∆z2,δyui,jk=ui−1,jk−2ui,jk+ui+1,jk,δzui,jk=ui,j−1k−2ui,jk+ui,j+1k.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 termry/4jk0rz/4jk0δyδ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(7)1−ry4jk0−112δy1−rz4jk0−112δzui,jk+1=1+ry4jk0+112δy1+rz4jk0+112δzui,jk.Equation (7) can be solved by two splitting steps given by(8a)1−ry4jk0−112δyui,jk+1/2=1+ry4jk0+112δyui,jk,(8b)1−rz4jk0−112δzui,jk+1=1+rz4jk0+112δzui,jk+1/2.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

Addingry/4jk0rz/4jk0δyδz on the left side and the right side of (7), we obtain(9)1−ry4jk0−112δy1−rz4jk0−112δzuk+1=1+ry4jk0+112δy1+rz4jk0+112δzuk−ry4jk0rz4jk0δyδzuk+1−uk.We consider the following:(10)Q=1−ry4jk0−112δy1−rz4jk0−112δz,R=1+ry4jk0+112δy1+rz4jk0+112δz+ry4jk0rz4jk0δyδzuk,T=−ry4jk0rz4jk0δyδz.Then, (13) can be rewritten as(11)Quk+1=Tuk+1+R.(11) is a linear system that can be iteratively solved in the following form:(12)Qun+1k+1=Tunk+1+R,where n denotes the nth iteration.If the initial guess,u0k+1, is set to uk, (7) can be recovered from (12). Equation (7) is the first iterative solution of (12). The scheme in (12) is straightforward. For the n+1th iteration, we obtain(13)1−ry4jk0−112δy1−rz4jk0−112δzun+1k+1=1+ry4jk0+112δy1+rz4jk0+112δzuk−ry4jk0rz4jk0δyδzunk+1−uk.Afterwards, (13) is split into two steps:(14a)1−ry4jk0−112δyun+1tmp=1+ry4jk0+112δy1+rz4jk0+112δzuk+ry4jk0rz4jk0δyδzunk+1−uk,(14b)1−rz4jk0−112δzun+1k+1=un+1tmp.(17) is the proposed iterative high-accuracy ADI method for the 3D parabolic equation which is depicted in Figure 1.Figure 1
The proposed iterative high-accuracy ADI method for 3D-PE. The green circular dots represent the elements that are updated in (14a) at a time. The blue circular dots represent the elements that are updated in (14b) at a time. Step 1 is explicit in the z direction and implicit in the y direction. Step 2 is explicit in the y direction and implicit in the z direction.


## 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μ of the iteration matrix(15)S=K−B−1J−A−1AB=K−B−1J−A−1ABK−B−1K−B=K−B−1CACBK−B,satisfy μ≤1, where J=1+1/12δy, A=ry/4jk0δy, K=1+1/12δz, B=rz/4jk0δz, CA=J−A−1A, and CB=K−B−1B. S is similar to CACB. 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 λJ≥0. The eigenvalues Cjω of CA then satisfy(16)Cjω=jωλJ−jω=ωλJ+ω2<1.The matrixCACA∗=J−A−1AA∗J∗−A∗−1=−J−A−1A2J+A−1=−J−A2−1A2=-CA2 has eigenvalues 0≤−Cjω2=ω2/λJ2+ω2<1 so that(17)CA2=maxλλCACA∗=maxωω21+ω2<1,μ≤CACB≤CACB<1.The CN-PE and the MF-ADI methods are unconditionally stable, which is estimated via von Neumann [23] by the amplification factor G≤1. The stability of the proposed IHA-ADI method can be estimated via von Neumann, amplification factor G satisfies.(18)G=1−Ry−1/12sin2σ1∆y/21−Rz−1/121+Ry+1/12sin2σ1∆y/21+Rz+1/12sin2σ2∆z/2−RyRzsin2σ1∆y/2sin2σ2∆z/2sin2σ2∆z/2sin2σ2∆z/2−RyRzsin2σ1∆y/2sin2σ2∆z/2≤1,where Ry=ry/4jk0 and Rz=rz/4jk0. 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(19)ψx,y,z=Ae−jkxx+kyy+kzz,where A is a constant coefficient, and ux,y,z is expressed as(20)ux,y,z=Ae−jkxx+kyy+kz−k0z,where k0=kx2+ky2+kz2.The normalized numerical wavenumberk0¯ can be extracted by substituting (20) into a different PE scheme. By substituting into different PE scheme, the normalized numerical wavenumber k0¯ 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(21)kx¯=k0¯cosφsinθ,ky¯=k0¯sinφsinθ,kz¯=k0¯cosθ.The numerical dispersion relations for the proposed iterative high-accuracy ADI method are given by(22a)ejk¯0,1nλ1−cosθαπ/Nλ2=1+jα−13sin2k¯0λsinφsinθ2Nλ1+jα−1/3sin2k¯0λcosφsinθ/2Nλ1−jα+1/3sin2k¯0λcosφsinθ/2Nλ+α2sin2k¯0λsinφsinθ/2Nλ∗sin2k¯0λcosφsinθ/2Nλ∗Mn−11−jα+1/3sin2k¯0λcosφsinθ/2Nλ,(22b)ejk¯0,2nλ1−cosθαπ/Nλ2=11−jα+1/3sin2k¯0λsinφsinθ/2Nλ,where α=∆x/k0∆z2, ∆z=∆y, and Nλ=λ/∆z=2π/k0∆z, M0=1,(23)Mn=1+jα−1/3sin2k¯0λsinφsinθ/2Nλ1+jα−1/3sin2k¯0λcosφsinθ/2Nλ1−jα+1/3sin2k¯0λcosφsinθ/2Nλ1−jα+1/3sin2k¯0λsinφsinθ/2Nλ,+α2sin2k¯0λsinφsinθ/2Nλ1−jα+1/3sin2k¯0λcosφsinθ/2Nλsin2k¯0λcosφsinθ/2NλMn−1−11−jα+1/3sin2k¯0λsinφsinθ/2Nλ.


## 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:Figure2 illustrates the relative error of the three methods with respect to θ for one iteration, where α=5, Nλ=4, and φ=45. It can be seen that the relative error increases with the increase in θ. The relative error differences for the three methods are higher when θ 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 2
Relative error versus the zenith angleθ.Figure3 illustrates the relative error with respect to α when Nλ=4, θ=6, and φ=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 α. It can be deduced from section II that α 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 α.Figure 3
Relative error versusα.Figure4 illustrates the relative error function of the azimuth angle, where θ=6, Nλ=4, and α=5. It can be observed that the wave-propagation along the diagonal (φ≈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 4
Relative error versus the azimuth angleφ.To better illustrate the effect of the number of grid points per wavelengthNλ and the zenith angle θ on the obtained results, Figure 5 shows that the relative errors of the three methods for θ=2° or 6°, α=5, and φ=45. It can be seen that the relative error of the three methods decreases with the increase in Nλ, 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λ=4, α=5, and φ=45, θ=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 5
Relative error versusNλ.Figure 6
The effect of the number of iterations.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 is60λ∗60λ. 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λ=4 and α=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.
MethodMemory (KB)CPU time (s)CN-PE256.832.152MF-ADI2.00.988IHA-ADI (1 iteration)2.11.177IHA-ADI (2 iteration)2.11.283IHA-ADI (3 iteration)2.11.515Figure 7
Fields in a rectangular tunnel (Z = 50 m): (a) the analytical result, (b) the MF-ADI method, (c) the proposed IHA-ADI method (1 iteration), and (d) the proposed IHA-ADI method (2 iterations).
(a)(b)(c)(d)


## 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.



---

*Source: 9955888-2023-03-22.xml*