The Laguerre finite difference oneway equation solver
Abstract
This paper presents a new finite difference algorithm for solving the 2D oneway wave equation with a preliminary approximation of a pseudodifferential operator by a system of partial differential equations. As opposed to the existing approaches, the integral Laguerre transform instead of Fourier transform is used. After carrying out the approximation of spatial variables it is possible to obtain systems of linear algebraic equations with better computing properties and to reduce computer costs for their solution. High accuracy of calculations is attained at the expense of employing finite difference approximations of higher accuracy order that are based on the dispersionrelationshippreserving method and the Richardson extrapolation in the downward continuation direction. The numerical experiments have verified that as compared to the spectral difference method based on Fourier transform, the new algorithm allows one to calculate wave fields with a higher degree of accuracy and a lower level of numerical noise and artifacts including those for nonsmooth velocity models. In the context of solving the geophysical problem the poststack migration for velocity models of the types Syncline and Sigsbee2A has been carried out. It is shown that the images obtained contain lesser noise and are considerably better focused as compared to those obtained by the known Fourier Finite Difference and PhaseShift Plus Interpolation methods. There is an opinion that purely finite difference approaches do not allow carrying out the seismic migration procedure with sufficient accuracy, however the results obtained disprove this statement. For the supercomputer implementation it is proposed to use the parallel dichotomy algorithm when solving systems of linear algebraic equations with blocktridiagonal matrices.
rmkRemark \newproofpfProof \newproofpotProof of Theorem LABEL:thm2
1 Introduction
The wave equation describes waves propagating in all directions. A oneway wave equation, also known as paraxial or parabolic, wave equation is an equation describing only downward propagating waves or only upward propagating waves for a positive or a negative vertical component of propagation velocity, respectively. Mathematical models based on the wave equation are often considered in problems of seismic prospecting [1, 2, 3], ocean acoustics [4, 5, 6] as well as for setting nonreflecting boundary conditions [7, 8, 9], whereas initially the ideas were formulated for solving the problem of electromagnetic wave propagation along the Earth’s surface[10].
For solving the oneway equation it is possible to use finite difference methods [11, 12]. To this end, the pseudodifferential squareroot operator was first approximated by an optimized series which has its origin in a continued fraction expansion [11, 13, 14, 15]. A finite difference method can handle strong lateral velocity variations, however it is diplimited and has significant numerical dispersions and artifacts. A conventional approach to reducing numerical dispersions is in applying finite difference schemes of a high order of accuracy, but when solving the oneway equation there arise essential computational problems. First, there is a numerical instability of calculations in using schemes of high orders of accuracy for the approximation of spatial derivatives towards the wave field extrapolation. Another difficulty is in that for decreasing the time of calculation, a practical approach, independently proposed in [16] and in [17], to solving extrapolation equation (4) is implemented by the MarchukStrang splitting [18]. However for an inhomogeneous velocity model such an approach provides only the first accuracy order as related to a step in depth because the operators to be split are not mutually commutative. In order to attain second order of accuracy, the splitting procedure should be specially symmetrized [18] thus twice increasing computer costs. The abandonment of using the MarchukStrang splitting procedure as regards this algorithm demands higher computer costs, especially for a 3D case when it is necessary to use iterative procedures for solving systems of linear algebraic equations (SLAEs). Another approach to reducing numerical dispersion and computing artifacts when applying the finite difference approximation is in the use of the filtration procedures [19] as well as of the nonlinear procedures [20] similar to the TVDschemes concepts for hyperbolic equations.
For reducing the numerical dispersion when solving geophysics problems the Fourierbased methods were developed [21, 22, 23]. These methods only approximately take into account the lateral velocity variations as being obtained on the basis of nonidentical mathematical transformations thus unpredictedly affecting the quality of solution. Nevertheless, for real practical tasks it has been just these methods that make possible to carry out calculations with sufficiently large mesh sizes. In this case numerical artifacts generated by the Fourierbased methods are essentially smaller than those generated by finite difference methods of second order of accuracy. On the other hand, the advantage of the finite difference approximation is in that it possesses convergence and is spatially localized, thus allowing the evaluation of the quality of solutions that are obtained with a sequence of imbedded meshes. In [24], an alternative approach to obtaining solutions similar to that of the oneway equation is proposed. This approach is a combination of various ideas related to the oneway equations, halfspace stiffness relation, special finite element discretization and complex coordinate stretching. However, the question of constructing a stable difference method of a high accuracy order for the classical oneway equation is still an open question and urgent both from the standpoint of applications and methods.
In this paper we propose a new finite difference algorithm of high order of accuracy for solving the oneway equation. As opposed to the generally accepted approach we use the integral Laguerre transform [25] instead of Fourier transform with respect to time. As for the approximation of spatial derivatives we employ only finite difference schemes. The approach in question was investigated in solving dynamic problems for the acoustics and elasticity equations [26, 27, 28, 29]. Applying the Laguerre transform as well as Fourier transform with respect to time enables us not to consider timecentering of the method of calculating the wave field that is required in analyzing finite difference approximations of the temporal derivative. Also, the Laguerre transform makes possible to obtain SLAEs with considerably better computational properties, than with using Fourier transform. Decreasing computer costs for solving SLAEs allows the abandonment of the MarchukStrang splitting and hence an increase both in stability and accuracy of the algorithm as a whole. In this case the approximation error relatively the mesh size towards the extrapolation of the wave field is of fourth accuracy order thus considerably reducing the numerical dispersion as compared to algorithms accurate to second order. As a result, a widespread opinion that a purely finite difference approach does not allow obtaining high quality solutions to the oneway equation will be invalidated.
2 Description of the method
2.1 Governing wave equations
The oneway equation can be obtained from the acoustics equation
(1) 
where is the field variable, is the wave velocity, the vertical direction is the extrapolation direction, i.e., the direction of oneway propagation and the positive axis is directed downward, i.e., toward increasing depth. Applying Fourier transforms to equation (1) over and yields
(2) 
where is a wave component at the radial frequency , is the horizontal wave number. Then, the 2D oneway acoustic wave equation for upcoming waves in the frequency domain can be expressed as:
(3) 
where is the imaginary unit. In practice, the squareroot operator is approximated by an optimized series which has its origin in a continued fraction expansion that can be represented by the ratios of polynomials [11].
(4) 
and the polynomial coefficients for the propagation angle should be optimized[13, 30].
Let us write down the oneway equation for the spatialtemporal domain . To this end it is convenient to introduce the auxiliary functions[14, 15] of the form
so that equation (4) can be rewritten as
Applying the inverse Fourier transform to the two latter equations, we come to a system of equations determining the downcontinuation process
(5a)  
(5b) 
2.2 Temporal approximation
Let us consider the direct and the inverse Laguerre transforms[31]
(6) 
Here are the orthonormal Laguerre functions, is Laguerre polynomial degree, is the order of the Laguerre functions, and is the transformation parameter. Assuming , we can show [31, 32] that
(7) 
where for we have
(8) 
Then applying the Laguerre transform (6) to equations (5a),(5b) we obtain the following system of equations for the calculation of the th coefficient of expansion:
(9a)  
(9b) 
where , and the index denotes number of a term of series (6).
The expansion parameters from (4),(5) for are chosen the following: , , , , , , for which, as shown in [30, 2], such approximation is valid up to the angles of degrees. It is also possible to select large values , thus asymptotically approaching to , however in this case computer costs will essentially increase.
The method in question can be considered as analog of the spectraldifference one based on Fourier transform, but the part of ”frequency” is played by the parameter , which determines the degree of the Laguerre polynomials. An important property is in that opposed to the Fourier method the operator of problem (9) is independent of the parameter separating harmonics . Hence it appears possible to solve in the domain the common system of equations for different righthand sides whereas this property is not valid in the domain . After carrying out the spatial approximation of equations (9), the multiple solution to the SLAE with a common matrix allows one to construct an effective computing procedure for solving a difference problem.
A disadvantage of the Laguerre transform is the absence of the fast transformation algorithm. However, taking into account the fact that input data are set only along the upper surface (), and the inverse transformation is done for a fixed time instant, the total cost of the direct and inverse transformations appears to be minor as compared to that needed for calculation of coefficients of decomposition (6) from the solution to problem (9).
2.3 Spatial approximation
To approximate equation (9a) we will use unconditionally stable the CrankNicolson scheme[33] of second order of accuracy:
(10a)  
(10b) 
where the difference operator is of the form
(11) 
Formally, the presence of a singular component in the solution [14, 15] does not allow attaining a high order of convergence for finite difference schemes, but for practical purposes the calculations have shown (see Section Numerical Experiments) that schemes of high orders of accuracy make possible to obtain solutions with a considerably lesser number of mesh nodes because such schemes more precisely reproduce the dispersion law. In deriving equations (5) the velocity model was assumed to be homogeneous, although it also yields satisfactory results for inhomogeneous media. In the latter case this model correctly keeps kinematics of waves, but not their amplitudes. For a wide range of problems such an approximate oneway model is admissible, as the correct account of amplitudes essentially increases computer costs [34]. Based on the abovesaid for approximating it is reasonable to use the dispersionrelationshippreserving method (DRP) by Tam and Webb [35], whose main idea is in the following. According to the Fourier derivative rule, values of optimized coefficients in (11) are defined as solution to the problem of minimizing the error functional in the space of wave numbers
where is some weight function, is a maximum accurate wave number. This approach and its various modifications [36, 37, 38] make possible to decrease the number of mesh nodes and to preserve high accuracy of calculations as compared to conventional difference schemes obtained with the Taylor expansion in series [39].
2.4 Solution of the SLAEs
Solving the difference equations as a stage of the algorithm proposed demands high computer costs, therefore its applicability depends on the speed of solving the corresponding SLAEs. As opposed to Fourier transform, the coefficients of the Laguerre expansion in series (6) are dependent in a recurrent manner (see (9)). Hence, for a fixed for different it is required to solve SLAEs many times with a common real matrix and different righthand sides. This allows the use of the methods based on decomposition, where factorization of a matrix for all is done only once.
Let us write down the difference problem (10) in the form of a SLAE as
(12) 
where is the unit matrix. Employing the Schur complement [40], the mesh functions can be calculated through the solution to the following reduced SLAE
(13) 
where the functions are defined as
(14) 
Making use of the matrix property [41]
(15) 
we can make a simplification
(16) 
where
Multiplying equation (16) by the matrix and taking into consideration the commutative property of , we obtain the governing equation for the calculation of mesh functions
(17) 
Matrix (17) is banded and can be explicitly represented without calculation of the matrices thus allowing us to apply efficient algorithms for solving SLAEs based on decomposition. Hence, we have to calculate different decompositions for (17), where is the number of mesh nodes in the direction . Note, also, that opposed to (16), for calculating the vector in the righthand side (17) the inversion of the matrix is not required.
After the calculation of the mesh functions before turning to calculating the functions , the functions can be expressed as
(18) 
Equation (18) is obtained applying the property of (15) to (14). The function should be calculated at the integer mesh nodes, because this simplifies the implementation of the Richardson extrapolation algorithm considered in the next Section.
2.5 Increasing accuracy of the spatial approximation via the Richardson extrapolation
According to the studies [27, 28, 29] it is needed to have a few tens of mesh nodes per a minimum wavelength in order to provide a reasonable accuracy of calculating schemes of second order for real spatialtemporal scales. In order to reduce computer costs, we increase the accuracy of the algorithm up to based on the Richardson extrapolation [42, 43].
Let the mesh functions , , defined on the meshes with the mesh steps and , be solutions to problem (10). Then the linear combination
(19) 
approximates the solution to problem (9) on the mesh accurate to . The Richardson extrapolation occurs more rarely than difference schemes of high accuracy orders. However, preliminary calculations have shown that multistep methods of the AdamsMoulton and the AdamsBashfort types do not provide stability even for small steps . The implicit RungeKutta methods of high accuracy orders are not efficient because they demand multiple calculation of the righthand side of the equation to be solved. The instability of methods of high orders is due to the presence of a singular component in the solution, but according to computational experiments the Richardson extrapolation makes possible to stabilize such instability.
Fourth order downwardcontinuation algorithm with the global correction To calculate the mesh functions accurate to , the following is necessary:

Based on the cubic splines interpolate values of the functions , preset on the mesh , into nodes .

Based on the cubic splines interpolate values of the functions ,, preset on the mesh , into semiinteger nodes .

On the mesh , applying equation (17), calculate the solution .

On the mesh , applying equation (17), calculate the solution .

Based on the Richardson extrapolation, correct the mesh function with (19).

After the calculation of the mesh functions for all depths , the functions can be defined from the solution to equations (18).

Turn to the calculation of , etc. coefficients of the expansion of the Laguerre series.
Also, we can formulate the second algorithm where the extrapolated value of (19) is calculated immediately after calculating the function , , , then the extrapolated value is used at the next step for calculating , , . This algorithm will be called ”the downwardcontinuation algorithm with the local correction”. If for calculating the solutions and one uses a numerically stable algorithm (the CrankNicolson scheme in our case), the extrapolation procedure based on the global correction will also be numerically stable, whereas in the general case, the stability is not ensured with the use of the local correction [43].
2.6 Parallel consideration
The possibility of the efficient use of modern supercomputer systems is one of obligatory demands imposed upon the development of new numerical algorithms. Therefore let us consider specific features of the parallel implementation of the approach proposed. As opposed to Fourier transform, applying the Laguerre transform does not allow the calculation in parallel of the function for different , as th coefficient of series expansion (6) recurrently depends on the th coefficient (9). Hence, it is needed to parallelize the algorithm at the stage of solving problem (10). In this case, the main difficulties are in solving SLAEs (17) and (18), whereas the parallel operation of multiplication by the matrices is not consuming from the communication standpoint. Let us consider the 2D data decomposition (Fig. 11a), where computing nodes contain processors having the shared memory access allowing implementation of innode interprocess communications considerably faster than the communications among the nodes. For solving SLAEs (17),(18) with banded matrices as a parallel algorithm it is reasonable to use the parallel dichotomy algorithm, which was developed for tridiagonal matrices [44] and blocktridiagonal matrices [28]. With respect to the number of arithmetical operations, the dichotomy algorithm is comparable with other available algorithms; however the time needed for interprocess communications is considerably less in the dichotomy algorithm as compared to other algorithms. This is because the implementation of the dichotomy process on a supercomputer reduces to calculating the sum of series for distributed data. The commutative and associative properties of addition enable a considerable reduction in the total computation time with the use of interprocessor interaction optimization algorithms. First, the parallel dichotomy algorithm was developed for solving SLAEs with the same tridiagonal matrix but different righthand sides. In [45], the dichotomy algorithm was applied to solving SLAEs with the Toeplitz tridiagonal matrices. It was shown that the Toeplitz tridiagonal matrices were able to effectively solve SLAEs both with one and several righthand sides. In [27, 29], the dichotomy algorithm was applied to implement a spectraldifference method to calculate acoustic and elastic wave fields. The authors could effectively use from up to processors for one calculation and to obtain a highly accurate numerical solution of the dynamic problem of elasticity theory. Thus, the dichotomy algorithm is a powerful instrument for solving SLAEs with tridiagonal and blocktridiagonal (banded) matrices.
The parallel dichotomy algorithm makes possible to use thousands of processors, however its efficiency depends on the number of simultaneously solved equations and the size of a SLAE matrix. However, such dependence is considerably weaker as compared to other parallel algorithms for solving tridiagonal and blocktridiagonal SLAEs. Due to dependence of the solution on it is impossible to solve problems (17),(18) for all depths in parallel, which would allow reducing communication costs at the expense of striping the interprocessor exchanges. Therefore the parallelization for direction should be done within one computational node profiting from the shared memory for the fast data exchange among processors thus allowing a considerable reduction of communication costs.
To increase the number of processors participating in one calculation it seems reasonable to carry out parallelization in direction. In this case calculations can be effectively implemented by the conveyor principle, i.e. as soon as at the node number the functions are calculated, their values and those of their second derivatives at the lower boundary of a subdomain are communicated to the computational node number . At the node with number , the calculation of functions with numbers starts, while at the node –functions with numbers . The value of the second derivative at the boundary of two subdomains is required for setting boundary conditions when constructing the cubic interpolation spline within the Richardson extrapolation procedure. Taking into account the fact that the number of terms of series (6) essentially exceeds the number of computational nodes, the load of such computational conveyor will be optimal. The conveyor communication costs will be minor because only one internode exchange for calculation of the functions for one value is required.
3 Numerical experiments
Let us discuss a few tests that would allow the evaluation of the quality of the solution to be obtained as compared with other known algorithms. Numerical procedures were implemented in Fortran90 using the MPI library. The calculation was performed on the ”Lomonosov” supercomputer of the Moscow State University. The supercomputer comprises Intel Xeon X5570 forecore processors operating at GHz in the Infiniband QDR communication environment. Each computational node contains two processors and 12 GB of RAM.
The coefficients of difference scheme (11) were chosen as follows[38]: , , , , , , , to provide the twelfth approximation order. The calculations have shown that as compared to classical difference schemes of high accuracy orders based on the Taylor expansion, the DRP approach makes possible to decrease the number of nodes in the difference scheme approximately by the factor of two with comparable accuracy. Taking into consideration the fact that computer costs for decomposition of matrices (17),(18) are proportional to the third degree of the matrix band width and forward and backward substitution of decomposition have costs proportional to the second degree, the effect of applying the DRP schemes is significant. Also, this allows one to decrease the demands for a minimum required random access memory for storing the matrices and .
3.1 Impulse Response
In the first test we illustrate analyzing the accuracy by impulse responses. For the calculation we used a homogeneous medium model with the speed and the size , the mesh steps being and . A point source was located at the center of the upper surface; the time dependence was as follows:
(20) 
where . As compared to the Fourier transform, where the basis functions are uniquely defined, the two parameters, and , should be set for using the Laguerre transform (6). These parameters were experimentally chosen based on the analysis of the convergence rate of the Fourier–Laguerre series for the shifted function with , where is the upper boundary of the time interval for which the wave field is calculated. The parameters and are chosen such that the function with in the meanquadratic norm is approximated accurate to . It should be regarded that as the value of grows, the values of the spatial derivatives of the functions will also increase. For this reason, the size of the mesh has to be simultaneously decreased in order to prevent growth of the error of the spatial approximation as the parameter grows. The number of addends in series (6) was for ; the expansion parameters were and .
Figures 1a,b present snapshots of the wave field from the point source for equation (4), which were obtained employing the finite difference algorithm (FD) of second accuracy order using the MarchukStrang splitting. This software procedure has been implemented in the available program package Seismic Unix. It is evident (Figs. 1a,b) that such an approach does not provide any accuracy of calculations due to multiple numerical noise and artifacts. In this case the transfer from the mesh with the step (Fig. 1a) to the mesh with the step (Fig. 1b) does not allow an increase in calculation accuracy up to the acceptable level.
A conceptually other situation arises (Figs. 1c,d,f,e) when the proposed Laguerre finite difference method (LFD) is used. It is clear that the new approach excludes the appearance of any numerical noise and the approximation error manifests as numerical dispersion which, as expected, is more distinct for the CrankNicolson scheme of second accuracy order (Figs. 1c,d). In Figs. 1e,f it is clear that applying the Richardson extrapolation makes possible to essentially improve the quality of solution. In addition, as for the second order scheme and for the fourth order of accuracy the transfer to a finer mesh allows diminishing the numerical dispersion. Thus, the key point in this case is the presence of convergence with decreasing the mesh step, which is not observed for calculations (Figs. 1a,b).
In spite of the absence of smoothness in the solution due to singularities [14, 15], the use of the method of fourth order of accuracy allows a decrease in numerical dispersion as compared to the scheme of second order of accuracy. However one should take into account the fact that with an equal number of mesh points the algorithm based on the Richardson extrapolation requires three times as much arithmetical operations as the CrankNicolson scheme. This is explained by the necessity of solving equations (10) on two meshes with the steps and . The question arises whether it is possible just to increase the number of mesh nodes by the factor of three and to make use of the second order scheme to obtain a result compatible in accuracy with the fourth order scheme. Fig. 2 shows the dependence of a wave field on the coordinate along the straight line ”Slice” (Fig. 1f). As is obvious, when the mesh step for the method of second accuracy order is three times smaller than that for the fourth order method (computer costs being approximately the same) the second order method yields less accurate results. Thus, it is reasonable to apply more complicated procedures based on the Richardson extrapolation, requiring large computer costs. It is also evident that the algorithm which is based on the local extrapolation inserts additional dispersion errors but better preserves amplitudes. Nevertheless, we are not going to use it because it is unstable for inhomogeneous media and inserts distinct artifacts associated with numerical dispersion into the image of a wave field.
When implementing the Richardson extrapolation algorithm we considered selecting the interpolation procedure and its influence on the solution to be solved. When using the linear interpolation the LFD algorithm demonstrated a very strong anisotropic dissipation in direction such that when passing a distance of several wavelengths the wave amplitude was tending to zero. The use of the parabolic or the cubic Lagrange interpolation somewhat improved the situation, but preservation of the amplitude was worse than for the cubic spline interpolation. Using the Lagrange interpolation of high orders as well as using splines of the orders exceeding the third caused instability of calculation. Thus, the choice of the cubic interpolation spline can be considered to be appropriate as other interpolation algorithms are either unstable or essentially dissipative.
Although the stability of the algorithm proposed has not been proved yet due to the difficulty in making the analytical analysis. In view of the following: asymmetry of matrices, the presence of singularities in a solution, applying the Richardson extrapolation procedure spline interpolation, etc. the algorithm in question in numerous experiments for different velocity models appeared to be more stable provided and a local contrast of the velocity medium model is not quite distinct. A higher level of stability as compared to other algorithms is explained by the fact that approximation errors in the direction of the wave field extrapolation are dominantly of dissipative character.
3.2 Migration procedure
To form the image of the Earth’s interior by the recorded seismic data one can use different algorithms of seismic migration [2, 11, 46, 47] including those of the wave field downward continuation algorithms. Let us consider the poststack migration based on the model of explosive boundaries [1]. This model allows the evaluation of the algorithm proposed as compared to other known methods. As many studies have dealt with the problem in question, we will focus on its computational features. As there is no analytical solution, we will make a comparison by the FD method, with the numericalanalytical Fourier Finite Difference (FFD) method [23] and the PhaseShift Plus Interpolation method (PSPI) [22] that are widespread in computational seismic prospecting. The two latter algorithms were developed instead of the FD method because it does not provide the required accuracy of calculation, and what is more important the computational stability when solving the oneway equation. However the FFD and PSPI methods approximately calculate the solution to the oneway equation for an inhomogeneous velocity model. The reason is that for the PSPI method the solution is obtained as a combination of several solutions for inhomogeneous models with different velocities. In the context of the FFD method, the input operator is represented as sum of two operators, one of them corresponding to a homogeneous medium with some averaged in horizontal velocity (usually, a minimum velocity is used as the reference velocity). The other operator represents a varying velocity correction. In the context of the MarchukStrang splitting procedure the first operator is inversed by analytical methods and the solution for correction is defined by finite difference methods of second order of accuracy. Although such a change is not similar to the original problem in terms of mathematics and is unstable for strongly contrast media, such an approach makes possible to considerably decrease computer costs. On the other hand, impossibility of evaluating the level of numerical errors on a sequence of imbedded meshes is an essential disadvantage of the FFD algorithm.
3.2.1 Syncline model
Theoretical seismograms (Fig. 3b) for the syncline model (Fig. 3a) were obtained with the help of the Gaussian beams algorithm [48, 49] implemented in the package Seismic Unix. For setting the boundary condition on the upper surface, the function for zerooffset section . Was expanded in series (6) with the parameters , and for . The calculations were carried out on meshes with the steps , and . According to the model of explosive boundaries the calculation velocities were set to be half the true velocity of the medium model.
As well as in the calculation of the wave field from the point source (Fig. 1) inadmissible level of numerical noise in the context of the migration problem for the FD method is observed (Fig. 4). In the FFD method, the noise level is much lower (Fig. 5), however, there are numerical artifacts and, in addition, the second horizon is insufficiently focused. When turning from the mesh with the step (Fig. 5a) to the mesh with the step (Fig. 5b) manifestation of artifacts remained as previously. Thus, the FFD algorithm can generate artifacts but does not allow their a posteriori evaluation on a sequence of imbedded meshes.
A conceptually different situation is observed with the LFD algorithm (Fig. 6), where the absence of highfrequency noise in the images obtained is distinct. The finite difference approximation errors for direction manifest as effect of numerical dispersion, whereas for direction the errors are of dissipative nature. The influence of these errors in turning to a finer mesh rapidly decreases, and the solution convergence of the mesh step makes possible to separate approximation errors from those of input data and inaccuracy of a model as it is. Hence, analyzing the results obtained when solving practical tasks can be essentially simplified. Additionally, calculations for a finer mesh with the step for all the three algorithms were carried out. The FD method for a small mesh step appeared to be unstable, whereas the quality of the image for the FFD method remained unchanged because the level of noise and artifacts was comparable with the calculation results with the mesh steps and . For the LFD method, the image in Fig. 6d does not visually differ from that obtained with the mesh step (Fig. 6c). Thus, with allowance for the convergence of the LFD algorithm, the result obtained with the mesh step can be considered to be sufficiently accurate.
3.2.2 Sigsbee Model
For Sigsbee2Amodel [50] (Fig. 7a), theoretical seismograms (Fig. 7b) were calculated using the algorithm of explosive boundaries implemented in the Madagascar package [51]. For setting the boundary condition, the function for zerooffset section was expanded in series (6) for , parameters and for . The calculations were done on the meshes with the steps , and .
According to the results of computational experiments, the FD method appeared to be unstable for the given velocity model and selected numerical parameters. Different degrees of smoothing the velocity medium model did not allow attaining the calculation stability for this method, therefore the FD method within this test is excluded from further consideration.
The PSPI method from the Seismic Unix package has shown unsatisfactory results because of bad energy focusing (Fig. 8a). This method was unable to correctly map boundaries of the complex salt shape (Fig. 8b). However note that the point type diffractors underlying the salt inclusion are correctly focused. The FFD method as compared to the PSPI has allowed obtaining a clearer image of the salt inclusion (Fig. 9b), but the image additionally contains numerous artifacts (Fig. 9a) due to approximate character of calculation formulas of the FFD method. The fact is, as well as in the previous test, an increase in resolving power of the mesh does not bring about an increase in calculation accuracy.
As opposed to the FFD and the PSPI methods, a conceptually different situation arises with the use of the LFD method: as is seen in Fig. 10, the LFD method does not contain additional artifacts but only those which are simultaneously present for the FFD and the PSPI algorithms, i.e. the artifacts being present in input data of the problem (zerooffset section). Using the DRP approach to determine coefficients of the finite difference approximation makes possible to attain high accuracy of calculation and to decrease numerical dispersion for direction. For direction, typical is numerical dissipation (Fig. 10a,b) that brings about smoothing an image. However, the numerical dissipation can be decreased at the cost of choosing a smaller step (Fig. 10c,d). Fourth approximation order due to the Richardson extrapolation allows the use of admissible steps for attaining the required calculation accuracy. Thus, errors caused by the LFD algorithm are controlled by choosing the step of a spatial mesh, whereas the choice of computational parameters for the PSPI and the FFD methods is a more ambiguous task, and decreasing a spatial mesh step for these methods does not really affect the imagery quality.
Thus, in the context of the poststack migration using the method in question allows obtaining a more accurate image of the Earth’s interior as compared to widespread mesh algorithms of solving the oneway equation. The method proposed can be also implemented for the prestack migration with a view to obtain high quality images. As opposed to numericalanalytical approaches requiring fulfillment of Fourier transform both with respect to time and with respect to space, the use of purely difference spatial approximation will enable us in the future to carry out calculations, first, for nonuniform meshes and, second, to implement a curvilinear boundary for the upper surface. The convergence of the new algorithm of the mesh step makes possible to evaluate the solution accuracy on a sequence of imbedded meshes thus separating approximation errors from errors of setting input data of the problem. For most of numerical migration methods it is necessary to preliminarily smooth a velocity medium model in order to provide the calculation stability. This brings about errors in the wave field kinematics and presents a severe problem in choosing a smoothing procedure. As the computational experiments show, the proposed LFD method for considered models and chosen computational parameters does not require preliminary smoothing of the medium model function, which counts in favor of a higher degree of stability as compared to existing mesh methods.
3.3 Performance consideration
As a matter of fact, a higher accuracy is not attained without supplementary computer costs. The singleprocessor implementation of the method proposed has shown that the LFD algorithm requires from to times more computational time than the PSPI, the FD, the FFD algorithms. The costs of the LFD algorithm for solving the SLAEs with banded matrices are essential, but the possibility of using optimized library functions for solving the SLAEs allows an increase in the method performance on the whole. Higher computer costs of the LFD algorithm can be justified because the FD method does not provide any accuracy and is weakly stable, and the FFD and the PSPI methods solve a certain approximate problem but not the oneway equation thus generating multiple noise and artifacts of the wave field.
For testing the performance of the method proposed the calculation time for a different number of mesh nodes and the number of processors was assessed. The use of multiprocessor systems makes absolute time costs insignificant because with a maximum number of processors the calculation time for all the test problems does not exceed several seconds (Table 1).
num proc.  
12  247  1008  44  119  944 
24  122  446  21  62  472 
48  76  294  14.8  41  313 
…  …  …  …  …  … 
216  11  46  2.4  6.2  49.6 
228  10  43  2.3  5.9  47.3 
240  9.7  40  2.2  5.6  44.9 
The dependence of the speedup value on the number of processors for all the tests is presented in Fig. 11b, indicating to the linear character of the dependence, for the first test attaining the superlinear speedup. A similar effect of superlinear speedup was observed when using the parallel dichotomy algorithm within the implementation of the alternating implicit direction method for solving the 2D Poisson equation[44]. Thus, the dichotomy algorithm implemented in the context of the approach proposed for solving the 2D oneway equation allows the effective use of a considerable number of processors. We can conclude that the chosen data decomposition of the problem (Fig. 11b) is fairly appropriate, as multiple exchanges for xdirection are instantly made at the expense of using the shared memory within one computational node, while slow internode communicative interactions for direction for solving problem (17) are done only once.
4 Conclusion
One of the major results of this research is in managing to construct a stable, efficient, purely finite difference method of high accuracy order for solving the oneway equation. By now this has presented considerable difficulties in numerical analysis due to the necessity of providing high accuracy and at the same time the calculation stability. Replacing Fourier transform by the Laguerre transform after employing the finite difference approximation of spatial derivatives, using algebraic transformations made possible to obtain the banded SLAEs appropriate for solutions by direct methods, in particular, by the parallel dichotomy algorithm. The rejection from using the splitting of the MarchukStrang type for the problem operator has excluded appearance of numerous computing artifacts and has increased both the order of approximation and the stability of the method as a whole. The oneway equation does not allow a correct description of the wave field amplitude for inhomogeneous media, hence the optimization of values of finite difference scheme coefficients aimed at decreasing the phase errors is a reasonable approach. The numerical experiments have shown that such DRP schemes make possible to decrease the mesh step in horizontal approximately by the factor of two as compared to classical difference schemes based on the Taylor series. Implementation of the Richardson extrapolation procedure has not only increased the accuracy of the algorithm up to fourth order relative to the mesh step in depth, but also has provided the stability of the method as compared to using schemes of the AdamsMoulton type. Applying implicit RungeKutta methods of high orders of accuracy is unreasonable, as in this case it is required to solve SLAEs of much higher dimensions. In addition, as well as for the Richardson extrapolation, the interpolation of the wave field values for auxiliary mesh nodes will be needed.
The algorithm proposed can be generated to solving 3D problems as well. To this end, the efficient iterative procedure for solving SLAEs should be developed. This is because the parallel dichotomy algorithm, being a direct method, is not intended for banded matrices of a large width. Also, it is necessary to reject not only the splitting of the MarchukStrang type for the problem operator, but also the splitting along horizontal directions that is used for decreasing computer costs, but generates numerical anisotropy. The efficient algorithm for solving the 3D oneway equation will be presented in the next publication in the near future.
5 Acknowledgments.
The author is grateful to Gennady Zhernyak for fruitful discussions of the paper.
References
 [1] J.F. Claerbout. Toward a unified theory of reflector mapping. Geophysics, 36(3):467–481, 1971.
 [2] O. Yilmaz and S.M. Doherty. Seismic Data Analysis: Processing, Inversion, and Interpretation of Seismic Data. Society of Exploration Geophysicist, 2001.
 [3] D.A. Angus. The oneway wave equation: A fullwaveform tool for modeling seismic body wave phenomena. Surveys in Geophysics, 35(2):359–393, 2013.
 [4] F.D. Tappert. Wave Propagation and Underwater Acoustics, chapter The parabolic approximation method, pages 224–287. Springer Berlin Heidelberg, Berlin, Heidelberg, 1977.
 [5] F.B. Jensen, W.A. Kuperman, M.B. Porter, and H. Schmidt. Computational Ocean Acoustics. SpringerVerlag New York, 2011.
 [6] D. Lee, A.D. Pierce, and ErC. Shang. Parabolic equation development in the twentieth century. Journal of Computational Acoustics, 08(04):527–637, 2000.
 [7] E.L. Lindman. ”Freespace” boundary vonditions for the time dependent wave equation. Journal of Computational Physics, 18(1):66 – 78, 1975.
 [8] B. Engquist and A. Majda. Absorbing boundary conditions for the numerical simulation of waves. Math. Comp., 31:629–651, 1977.
 [9] L.N. Trefethen and L. Halpern. Wellposedness of oneway wave equations and absorbing boundary conditions. Mathematics of Computation, 47(176):421–435, 1986.
 [10] M. Leontovich and V. Fock. Solution of the problem of propagation of electromagnetic waves along the earthâs surface by the method of parabolic equation. Acad. Sci. USSR J. Phys., 10:13–24, 1946. (Engl. transl).
 [11] J.F. Claerbout. Imaging the Earth’s Interior. Blackwell Scientific Publications, Inc., Cambridge, MA, USA, 1985.
 [12] D. Lee and G. Botseas. Ifd: An implicit finitedifference computer model for solving the parabolic equation. Technical Report 6659, Nav. Underwater Syst. Ctr, 1982.
 [13] L. Halpern and L.N. Trefethen. Wideangle oneway wave equations. J. Acoust. Soc. Am., 84(4):1397–404, 1988.
 [14] A. Bamberger, B. Engquist, L. Halpern, and P. Joly. Higher order paraxial wave equation approximations in heterogeneous media. SIAM Journal on Applied Mathematics, 48(1):129–154, 1988.
 [15] A. Bamberger, B. Engquist, L. Halpern, and P. Joly. Parabolic wave equation approximations in heterogenous media. SIAM Journal on Applied Mathematics, 48(1):99–128, 1988.
 [16] G.I. Marchuk. Some application of splittingup methods to the solution of mathematical physics problems. Applik Mat, 13(2):103–132, 1968.
 [17] G. Strang. On the construction and comparison of difference schemes. SIAM Journal on Numerical Analysis, 5(3):506–517, 1968.
 [18] G.I. Marchuk. Handbook of numerical analysis. In P. G. Ciarlet and J. L. Lions, editors, Handbook of Numerical Analysis, volume 1, pages 197–460. Elsevior Science Publishers BV, 1990.
 [19] C. Bunks. Effective filtering of artifacts for implicit finitedifference paraxial wave equation migration1. Geophysical Prospecting, 43(2):203–220, 1995.
 [20] T. Fei and K. Larner. Elimination of numerical dispersion in finitedifference modeling and migration by fluxcorrected transport. Geophysics, 60(6):1830–1842, 1995.
 [21] J. Gazdag. Wave equation migration with the phaseshift method. Geophysics, 43(7):1342–1351, 1978.
 [22] J. Gazdag and P. Sguazzero. Migration of seismic data by phase shift plus interpolation. Geophysics, 49(2):124–131, 1984.
 [23] D. Ristow and T. Rühl. Fourier finitedifference migration. Geophysics, 59(12):1882–1893, 1994.
 [24] M.N. Guddati. Arbitrarily wideangle wave equations for complex media. Computer Methods in Applied Mechanics and Engineering, 195:65 – 93, 2006.
 [25] M. Abramowitz and I.A. Stegun. Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables. Dover, New York, ninth dover printing, tenth gpo printing edition, 1964.
 [26] B.G. Mikhailenko. Simulation of seismic wave propagation in heterogeneous media. Siberian J. of Numer. Mathematics, 6:415–429, 2003.
 [27] A.G. Fatyanov and A.V. Terekhov. Highperformance modeling acoustic and elastic waves using the parallel dichotomy algorithm. J. Comp. Phys., 230(5):1992–2003, 2011.
 [28] A.V. Terekhov. A fast parallel algorithm for solving blocktridiagonal systems of linear equations including the domain decomposition method. Parallel Comput., 39(67):245–258, 2013.
 [29] A.V. Terekhov. Spectraldifference parallel algorithm for the seismic forward modeling in the presence of complex topography. Journal of Applied Geophysics, 115(0):206–219, 2015.
 [30] Myung W. Lee and Sang Y. Suh. Optimization of oneway wave equations. Geophysics, 50(10):1634–1637, 1985.
 [31] L. Debnath and D. Bhatta. Integral Transforms and Their Applications, Second Edition. Taylor & Francis, 2006.
 [32] B.G. Mikhailenko. Spectral laguerre method for the approximate solution of time dependent problems. Applied Mathematics Letters, 12:105–110, 1999.
 [33] J. Crank and P. Nicolson. A practical method for numerical evaluation of solutions of partial differential equations of the heatconduction type. Proc. Camb. Phil. Soc., 43(1):50–67, 1947.
 [34] Y. Zhang, G. Zhang, and N. Bleistein. Theory of trueamplitude oneway wave equations and trueamplitude commonshot migration. Geophysics, 70(4):E1–E10, 2005.
 [35] C.K.W. Tam and J.C. Webb. Dispersionrelationpreserving finite difference schemes for computational acoustics. Journal of Computational Physics, 107(2):262 – 281, 1993.
 [36] W. Liang, Y. Wang, and C. Yang. Determining finite difference weights for the acoustic wave equation by a new dispersionrelationshippreserving method. Geophysical Prospecting, 63(1):11–22, 2015.
 [37] C. Chu and P.L. Stoffa. Determination of finitedifference weights using scaled binomial windows. Geophysics, 77(3):W17–W26, 2012.
 [38] J.H. Zhang and Z.X. Yao. Optimized explicit finitedifference schemes for spatial derivatives using maximum norm. Journal of Computational Physics, 250:511 – 526, 2013.
 [39] A.A. Samarskii. The Theory of Difference Schemes. Marcel Dekker, 2001.
 [40] R.W. Cottle. Manifestations of the schur complement. Linear Algebra and its Applications, 8(3):189 – 211, 1974.
 [41] H.V. Henderson and S.R. Searle. On deriving the inverse of a sum of matrices. SIAM Review, 23(1):53–60, 1981.
 [42] L.F. Richardson. On the approximate arithmetical solution by finite differences of physical problems involving differential equations, with an application to the stresses in a masonry dam. Proceedings of the Royal Society of London A: Mathematical, Physical and Engineering Sciences, 83(563):335–336, 1910.
 [43] G. I. Marchuk and V. V. Shaidurov. Difference methods and their extrapolations. Appl. Math. Springer, New York, NY, 1983.
 [44] A.V. Terekhov. Parallel dichotomy algorithm for solving tridiagonal system of linear equations with multiple righthand sides. Parallel Comput., 36(8):423–438, 2010.
 [45] A.V. Terekhov. A highly scalable parallel algorithm for solving toeplitz tridiagonal systems of linear equations. Journal of Parallel and Distributed Computing, 87:102–108, 2016.
 [46] S.H. Gray, J. Etgen, J. Dellinger, and D. Whitmore. Seismic migration problems and solutions. Geophysics, 66(5):1622–1640, 2001.
 [47] B. Biondi. 3D Seismic Imaging. Society of Exploration Geophysicists, 2006.
 [48] M.M. Popov. A new method of computation of wave fields using gaussian beams. Wave Motion, 4, 1982.
 [49] V. Cerveny. Gaussian beam synthetic seismograms. J. Geophys., 58:44–72, 1985.
 [50] J. Paffenholz, B. McLain, J. Zaske, and P.J. Keliher. Subsalt multiple attenuation and imaging: Observations from the Sigsbee2B synthetic dataset, chapter 538, pages 2122–2125.
 [51] S. Fomel, P. Sava, I. Vlad, Y. Liu, and V. Bashkardin. Madagascar: opensource software project for multidimensional data analysis and reproducible computational experiments. Journal of Open Research Software 1:e8, 2013.