key: cord-0058885-gcoyitbg
authors: Lytaev, Mikhail S.
title: Automated Selection of the Computational Parameters for the Higher-Order Parabolic Equation Numerical Methods
date: 2020-08-24
journal: Computational Science and Its Applications - ICCSA 2020
DOI: 10.1007/978-3-030-58799-4_22
sha: fa59a72805899dd49a689b2686c5bc113e4729a0
doc_id: 58885
cord_uid: gcoyitbg

This study is devoted to the features of the numerical methods for the parabolic wave equation. While seeking a numerical solution, it is necessary to select a set of computational parameters of the numerical method. The choice of the computational parameters affects the speed and accuracy of the calculations. Automation of the choice of computational parameters is useful when applying mentioned numerical methods in complex software systems, where the user cannot select them manually. In this paper, we consider a finite-difference split-step Padé method for the one-way Helmholtz equation. A discrete dispersion relation based algorithm for finding the optimal computational parameters of the numerical method is presented.

Approximation of the elliptic Helmholtz equation by the parabolic equation (PE) has gained a considerable popularity in the wave propagation study. In 1946, Soviet scientists M. A. Leontovich and V. A. Fock [21] proposed a wave-based approach to solving the tropospheric radio wave propagation problem using the PE. The original analysis of the PE method has been based on some physical considerations. Subsequently, a series of papers appeared [3, 9, 10, 15] , where the PE method is derived and analyzed more strictly from a mathematical point of view, relying mainly on the theory of operators.

The PE method is well suitable for wave propagation modeling in inhomogeneous waveguides. PE method is used for tropospheric radio wave propagation [22] , propagation of the acoustic waves in the underwater sound channel [18, 29] , acoustic and electromagnetic scattering [17, 32] , optics [14] , quantum mechanics [11, 26] and geophysics [35] . The main advantage of this method is the possibility of effective numerical implementation. Several software systems have been implemented based on the PE method: AREPS [5] , RAM [13] , CARPET, PETOOL [28] . A significant number of PE codes for underwater acoustics are collected in [30] .

Any numerical method has many artificial parameters, such as the geometry of the computational grid, the size of the integration area, the approximation order, the convergence threshold and so on. These parameters directly affect the accuracy and speed of calculations. Despite a significant number of works devoted to numerical methods for PE, they practically do not pay attention to the problem of choosing the computational parameters depending on the propagation conditions, as well as describing the limits of applicability of these numerical methods. It is usually assumed that the researcher can pick them up manually based on his experience. However, this approach is completely inapplicable when using these methods in complex software systems. Operators of such systems cannot manually select computational parameters for each propagation case. There is a need for a method that is able to automatically determine the computational parameters depending on the input data and the required accuracy of calculations. The purpose of the present research is to develop a method for automating the finding of the computational parameters of PE numerical methods depending on the input data. The subject of the study is the finitedifference rational approximations of the one-way Helmholtz equation [7, 22] .

The paper is organized as follows. Section 2 describes the numerical scheme under consideration with emphasis on spectral representation of the propagation operator. In Sect. 3, a discrete dispersion analysis of the numerical scheme is performed. Section 4 is devoted to the synthesis of the automated computational parameter selection algorithm. Numerical examples are given in Sect. 5.

A number of wave propagation problems can be reduced to the Helmholtz equation in the following form [22] 

where ψ(x, z) is a complex-valued wave field function; -m 2 (x, z) is the squared refractive index; -k = 2π λ is the wavenumber; -λ is the wavelength.

Equation (1) may be defined in an unbounded domain. We are seeking the solution of Eq. (1) on a finite domain Ω = [0, x max ] × [0, z max ]. Nevertheless, the bounds of Ω may be transparent [12, 19, 25] .

The wave field is generated by the initial condition

where ψ 0 (z) is a known function, which usually corresponds to the far-field aperture pattern of the wave field radiator. We proceed from the assumption that the generated waves propagate mostly in a preferred direction along the longitudinal coordinate x.

Evolution solution of Eq. (1) for the outgoing field can be written as follows [22] u

where

Thus, the step-by-step numerical method is based on the approximation of the pseudo-differential propagation operator (2) . The following designations are used in this paper z j = jΔz, x n = nΔx, u n j = u(x n , z j ). Using the definition of pseudo-differential operator [34] and Fourier transform by variable z, expression (2) can be rewritten as follows

whereũ

To simplify expression (4), we use the Padé approximant [1, 7] , written in the following form

Methods for calculating and analyzing Padé approximation coefficients for various elementary functions are described in [1] . For the numerical calculations of coefficients a 1 , . . . , a p , b 1 , . . . , b p , it is recommended to use computer algebra and symbolic calculations.

It will be shown later that the approximation (5) affects the maximum propagation angle and the maximum longitudinal step size Δx.

Substituting (5) into (4) and introducing auxiliary functions v n 1 , v n 2 , . . . , v n p−1 , the propagation operator in spectral domain can be approximately written as a system of equations

By applying the inverse Fourier transform to each row of system (6), we obtain the following system of differential equations

This system can be formally written as a function of operator L

Thus, the use of the propagation operator (2) at each step is reduced to the solution of the system of one-dimensional Eqs. (7) .

The system of Eqs. (7) is solved sequentially from top to bottom. Each row is a similar differential equation of the form

where u is a known function, which is obtained on a previous step, v is unknown function.

Thus, the original differential equation is reduced to a semi-discrete one. Its solution can be represented by a sequence of one-dimensional differential equations of the form (8).

Further, we construct a numerical scheme for Eq. (8). We take into consideration a uniform partition in the z direction with step size Δz and denote u(jΔz) = u j . Consider the second order difference operator

Operator δ 2 can be represented as follows

Thus, owing to the definition of function of differential operator [34] , δ 2 can be written as a function

Now we will construct the inverse operator. We will express the differential operator through the second difference operator. Consider the inverse function

Second derivative operator can be represented in the following form

This expression allows one to construct an arbitrary accuracy finite-difference approximation of the differential operator. Indeed, by decomposing expression (9) according to Taylor's formula, we get

Leaving the first two terms in the decomposition (10) and using Padé approximation, we get the following expression for the second derivative [20] 

where β = 0, for 2-d order approximation, 1 12 , for 4-th order approximation. Substituting this approximation into (8), we obtain the following difference equation

where m 2 j = m 2 (nΔx, jΔz). Equation (12) on a finite grid with discrete boundary conditions can be represented as a tridiagonal system of linear algebraic equations that admits an effective numerical solution using the tridiagonal matrix algorithm [31] .

Using an approximation of order higher than four in (10) is possible, but results in less sparse matrices.

In the case of a homogeneous medium (m(x, z) ≡ 1), longitudinal and transversal discretization can be joined using Padé approximation. Operator L in this case coincides with the second derivative operator. Substituting the dependence of the differential operator on the second difference (9) into the propagation operator (2) and applying Padé approximation to the resulting operator, which now depends on the second difference operator δ 2 , we obtain the following expression

Next, we apply the algorithm described above for a second-order scheme (β = 0), but using the new coefficients a l and b l .

It should be noted that another widely used method for solving PE is the split-step Fourier method [28] . The split-step Fourier method does a better job of approximating the diffraction part of the propagation operator and is preferred in cases when the boundary conditions do not have a significant effect on propagation. At the same time, this method meets serious issues with modeling both artificial and natural boundary conditions [24, 37] . Also, the split-step Fourier method is not applicable for accounting for strong spatial variations of the refractive index m 2 .

Note that Padé-[1/1] approximation of the propagation operator is equivalent to the Crank-Nicolson numerical scheme [22] for PE, thus the considered numerical method is a generalization of it.

The above reasoning can also be applied to the linear Schrodinger equation. It is sufficient to replace the square root operator √ 1 + L with 1 + L/2 in the expression (2).

When constructing the numerical scheme in the previous section, several additional parameters appeared to be defined. Namely: the grid step along the longitudinal and transversal coordinate (Δx and Δz) and Padé approximation coefficients a 1 . . . a p , b 1 

It is well known that the Crank-Nicolson numerical scheme and finitedifference Padé approximations considered in this paper are unconditionally stable [22] . Thus, the solution of the finite-difference equation tends to the exact solution as the sampling of the computational grid increases. At the same time, reducing the cells of the computational grid increases the required time and memory. Therefore it is of great interest from a computational point of view to determine the maximum longitudinal and transversal increments, which at the same time provide acceptable accuracy. It is common to use adaptive step size algorithms [31] to solve this problem, which consists in gradually reducing the grid steps until the required accuracy is achieved. However, since the numerical scheme under consideration requires the selection of several parameters simultaneously, an alternative method is required. This study proposes an alternative solution based on the analysis of discrete dispersion relations.

We analyze the dispersion relations [4, 16] for the discrete numerical scheme, described in the previous section. Keeping in mind replacement (3), consider the discrete plane wave

where k x = k cos θ is the horizontal wavenumber; -k z = k sin θ is the vertical wavenumber; -θ is the angle between the direction of the plane wave and the positive direction of the x-axis.

j , we substitute the expression (14) to the discrete scheme. Then, after some simple calculations, we derive the following dispersion relation

Dispersion relation for the original Helmholtz Eq. (1) is written as follows k 2

x + k 2 z = k 2 . Discrete horizontal wavenumber as a function of the computational parameters and the propagation angle θ is expressed as follows

Horizontal wavenumber for the original Helmholtz equation (1) is expressed by the following formula

Using expressions (15) and (16), one can estimate the phase error of the numerical scheme at each propagation step. The discrete dispersion relation (15) includes all the above mentioned computational parameters of the numerical scheme, which makes it possible to analyze it. Relation (15) gives an opportunity to analyze the accuracy of a particular scheme depending on the propagation angle. Figure 1 represents the dependence of the value k d x on the propagation angle for various approximations of the propagation operator. Corresponding relative error is depicted in Fig. 2 . One can easily see that the increased accuracy of the approximation leads to an increased maximum propagation angle.

A numerical scheme based on Padé approximation of the form [p/q] is stable under the condition 0 ≤ q −p ≤ 2 [6] . The scheme is A-stable under the condition p = q [6, 22] . Its disadvantage is that the corresponding values of k x are always real. As one can see from expression (16) , in the evanescent part of the spectrum (k z > k), the horizontal wavenumber of the original Helmholtz equation is purely imaginary. Thus, evanescent waves, which are formed when the waves are scattered on vertical inhomogeneities, do not attenuate when using the scheme [p/p]. This leads to artificial Gibbs oscillations of the solution [22, 25] . Schemes of the form [p − 1/p] and [p − 2/p] are L-stable, thus they suppress the occurring high-frequency oscillations. It can be seen that the scheme [3/4] really provides the fulfillment of the condition Im(k x ) > 0 for k z > k, i.e. evanescent waves will be damped. However, the behavior of k d x for scheme [3/4] is not the same as for original Helmholtz equation. In many propagation problems, the evanescent spectrum is not significant, and such regularization method is quite suitable. To correctly account for the evanescent spectrum, special Padé approximations of the square root should be applied [27, 36] . Their essence is to turn the branch cut of the square root to a certain angle α. Then, based on a real-valued approximation of the form [2]

where a l = 2 2p + 1 sin 2 lπ 2p + 1 , b l = cos 2 lπ 2p + 1 , 

The seeking of the optimal computational parameters consists of the following two stages: determining the maximum propagation angle θ max and solving the minimization problem. The analysis of the discrete dispersion relation demonstrates that increasing the maximum wave propagation angle increases the required density of the calculation grid and, accordingly, the calculation time. It should be noted that the required maximum propagation angle strongly depends on the propagation conditions, especially on boundary conditions. The angle between the propagating wave and the surface does not exceed a few degrees when propagating over a smooth surface, such as the surface of the sea. Inhomogeneities of the landscape require consideration of a wider angular spectrum. However, even in this case, the maximum angle rarely exceeds 10 • . Waves with propagation angle up to 90 • and evanescent waves must be considered in the presence of vertical obstacles that usually occur in the urban environment. Thus, based on the geometry of the propagation medium, one can determine the maximum propagation angle θ max .

Next, we need to solve the problem of computational parameter optimization. It is necessary to minimize the number of calculations of the Eq. (8), as well as the number of nodes of its discrete analog (12) . It corresponds to the minimization of the following expression g (Δx, Δz, a 1 . . . a p , b 1 

when

where ε is acceptable phase error at the distance x max from the source of radiation. Obviously, the speed of calculations is directly depend on the value of function g. The minimization algorithm is implemented by simply iterating over the various values of the computational parameters. This approach gives satisfactory results in terms of the time of the minimization procedure, however more advanced optimization methods may be applied in the future.

It should be noted that the dispersion analysis is carried out without taking into account the inhomogeneities of the refractive index. This assumption is not significant when the refractive index does not suffer strong spatial variations. However, the variations of the refractive index are often used to account for the spatial location of the obstacles [8, 24] . For path-through integration of the obstacles, it is noteworthy that the wavelength λ inside the medium with the refractive index n changes following the ratio

Therefore, it is necessary to proportionally reduce the size of the cells of the computational grid inside the obstacle. A non-uniform computational grid may be used for this purpose [33] .

Thus, simplified description of the algorithm is the following:

1. Estimate the maximum propagation angle θ max based on the geometry of the problem; 2. Minimize function (17) under condition (18) .

The first part of the algorithm significantly depends on the particular problem and its geometry. The second part does not depend on the specific problem.

In this section, we analyze the results of the presented algorithm for various values of the maximum propagation angle θ max and approximation orders. We consider the following Padé approximations: [1/1], [1/2] , [3/4] , [7, 8] and approximation of the form (13). 4-th order transversal approximation is used. The length of the computational domain x max = 2e4λ, width z max = 1e2λ. Acceptable phase error ε = 10 −3 . The results of the algorithm for the maximum propagation angle equal to 3 • , 10 • , 20 • and 45 • are given in Tables 1, 2, 3 and 4 respectively. These tables present the optimal cell sizes of the computational grid for various approximations, as well as the value of the function g. "-" means that reasonable computational parameters could not be selected.

Increasing the approximation order of the finite-difference numerical scheme has a significant impact on the performance. It is clearly seen that increasing the approximation order affects both the maximum propagation angle and the calculation speed. Padé-[7/8] approximation yields three orders faster numerical scheme than the standard Crank-Nicolson ([1/1]) approximation even in a narrow-angle case. Standard Crank-Nicolson numerical method can only be used for very small propagation angles. Propagation angles greater than 45 • can only be handled by the scheme (13). The numerical propagation method considered in this paper and the presented method for optimal computational parameter selection are implemented by the author as Python 3 open-source software library [23] . 

Discrete dispersion analysis of the finite-difference numerical scheme based on Padé approximation of the one-way Helmholtz equation was carried out. Based on this analysis, an algorithm for automatic computational parameter selection of the specified numerical scheme has been developed. The method is useful for the numerical solution of a class of wave propagation problems. The method is based on the theoretical analysis of the numerical scheme and uses a priori information about the propagation medium. The presented method is easier to implement and more effective than the adaptive selection of parameters. The presented method can be used in complex software systems that use wave propagation modeling.

Higher order paraxial wave equation approximations in heterogeneous media

Parabolic wave equation approximations in heterogenous media

Wide-angle alternatingdirection implicit finite-difference beam propagation method

Areps and temper-getting familiar with these powerful propagation software tools

Numerical Methods for Ordinary Differential Equations

A split-step pade solution for the parabolic equation method

A higher-order energy-conserving parabolic equqation for range-dependent ocean depth, sound speed, and density

Bremmer series that correct parabolic approximations

Higher-order parabolic approximations to timeindependent wave equations

Numerical solutions of the timedependent schrödinger equation in two dimensions

Discrete non-local boundary conditions for split-step padé approximations of the one-way helmholtz equation

Underwater Acoustic Modeling and Simulation

Exact transparent boundary condition for the parabolic equation in a rectangular computational domain

Uniform high-frequency approximations of the square root helmholtz operator symbol

Wide-angle beam propagation using padé approximant operators

Frequency-domain and time-domain solvers of parabolic equation for rotationally symmetric geometries

Computational Ocean Acoustics

Accurate artificial boundary conditions for the semi-discretized linear schrödinger and heat equations on rectangular domains

Numerical ocean acoustic propagation in three dimensions

Solution of the problem of propagation of electromagnetic waves along the earth's surface by the method of parabolic equation

Parabolic Equation Methods for Electromagnetic Wave Propagation. The Institution of Electrical Engineers

Python wave proragation library

Split-step padé approximations of the Helmholtz equation for radio coverage prediction over irregular terrain

Nonlocal boundary conditions for split-step padé approximations of the helmholtz equation with modified refractive index

Fourth order real space solver for the time-dependent schrödinger equation with singular Coulomb potential

Rational square-root approximations for parabolic equation algorithms

PETOOL: MATLAB-based oneway and two-way split-step parabolic equation tool for radiowave propagation over variable terrain

Transparent boundary conditions for iterative highorder parabolic equations

Ocean acoustics library

Multisector parabolic-equation approach to compute acoustic scattering by noncanonically shaped impenetrable objects

Nonuniform depth grids in parabolic equation solutions

Pseudodifferential Operators and Nonlinear PDE

The laguerre finite difference one-way equation solver

Complex padé approximants for wide-angle acoustic propagators

Applying the parabolic equation to tropospheric groundwave propagation: a review of recent achievements and significant milestones