The stabilization of high-order multistep schemes for the Laguerre one-way wave equation solver

TL;DR

This paper develops stable, high-order spectral-difference methods for solving the one-way wave equation using Laguerre transforms, employing spline interpolation for stability and predictor-corrector algorithms to reduce computational costs.

Contribution

It introduces stabilizing procedures based on spline interpolation for high-order multistep schemes solving the Laguerre one-way wave equation, improving stability and efficiency.

Findings

Stable high-order schemes achieved for 1D and 2D wave equations.

Efficient predictor-corrector method reduces computational costs.

Validated accuracy and stability through seismic migration applications.

Abstract

This paper considers spectral-difference methods of a high-order of accuracy for solving the one-way wave equation using the Laguerre integral transform with respect to time as the base. In order to provide a high spatial accuracy and calculation stability, the Richardson method can be employed. However such an approach requires high computer costs, therefore we consider alternative algorithms based on the Adams multistep schemes. To reach the stability first for the 1D one-way equation and then for the 2D case, the stabilizing procedures using the spline interpolation were developed. This made possible to efficiently implement a predictor-corrector type method in terms of which a boundary value problem for high-order elliptic equations is substituted for a sequence of inversions of second order elliptic operators thus decreasing computer costs. To assess the accuracy and stability of difference approximations for the 1D one-way wave equation, the analytical solution based on the double Laguerre transform was obtained. This solution can be efficiently calculated if for summation one makes use of fast algorithms of computing a discrete linear convolution. For the 2D one-way wave equation the stability and accuracy of the procedures proposed have been studied on implementing the migration algorithm within a problem of seismic prospecting.

Tables

Table 1. Table 1: Dependence of the error value ‖ u exact − u h ‖ 2 / ‖ u exact ‖ 2 subscript norm superscript 𝑢 exact superscript 𝑢 ℎ 2 subscript norm superscript 𝑢 exact 2 \|u^{\mathrm{exact}}-u^{h}\|_{2}/\|u^{\mathrm{exact}}\|_{2} on the number of mesh nodes for different methods.

	AM5-I5	AM6-I7	CN	RK4	Richardson	AM5-D4
$N_x$
1000	0.32	0.18	1.47	0.99	6.04e-02	0.56
1500	6.67e-2	2.4e-2	1.51	0.92	1.13e-2	0.18
2000	1.72e-2	4.6e-3	1.38	0.6	3.5e-3	6.5e-2
3000	2.3e-3	4.18e-4	0.87	1.68	6.82e-4	1.33e-2
4000	5.6e-4	7.52e-5	0.53	5.46e-2	2.14e-4	4.2e-3
4500	3.1e-4	3.72e-5	0.43	3.39e-2	1.33e-4	2.6e-3

Equations

$$\frac{\partial \tilde{u}}{\partial z}=-\mathrm{i}\frac{\omega }{c}\sqrt{1-{\left(\frac{c{k}_x}{\omega }\right)}^2}\tilde{u},$$

$$\sqrt{1-{\left(\frac{c{k}_x}{\omega }\right)}^2}\approx \left[1-\sum _{s=1}^n\frac{{\beta }_s{c}^2{k}_x^2}{{\omega }^2-{\gamma }_s{c}^2{k}_x^2}\right],$$

$${\partial }_{t}v+c{\partial }_{x}v=0,t>0,x\in \mathbb{R}$$

$$\text{mathds}L\{g(t)\}={\bar{g}}_m={\int }_0^{\infty }g(t)l_m(\eta t)dt,g(t)=\text{mathds}{L}^{-1}\{{\bar{g}}_m\}=\sum _{m=0}^{\infty }{\bar{g}}_m{l}_m(\eta t),$$

$$\text{mathds}L\left\{\frac{d}{dt}g(t)\right\}=\frac{\eta }{2}{\bar{g}}_m+{\Phi }_1({\bar{g}}_m),$$

$${\Phi }_1({\bar{g}}_m)=\sqrt{\eta }g(0)+\eta \sum _{j=0}^{m-1}{\bar{g}}_j.$$

$$\begin{array}[]{ll}\displaystyle\left(\frac{\eta}{2}+c\partial_{x}\right)\bar{v}^{m}+\Phi_{1}(\bar{v}^{m})=0,\ m=0,1,...,\end{array}$$

$${\Phi }_1({\bar{v}}^m)=\eta {\bar{v}}^{m-1}+{\Phi }_1({\bar{v}}^{m-1}),$$

$$\left(\frac{\eta }{2}+c{\partial }_x\right){\bar{v}}^0+\sqrt{\eta }\phi (x)=0,$$

$$\left(\frac{\eta }{2}+c{\partial }_x\right){\bar{v}}^m=\left(-\frac{\eta }{2}+c{\partial }_x\right){\bar{v}}^{m-1},m=1,2,...$$

$$\frac{c\left({\bar{v}}_{j+1}^m-{\bar{v}}_j^m\right)}{h_x}+\frac{\eta }{2}{\bar{v}}_j^m=\frac{c\left({\bar{v}}_{j+1}^{m-1}-{\bar{v}}_j^{m-1}\right)}{h_x}-\frac{\eta }{2}{\bar{v}}_j^{m-1}.$$

$${\widetilde{\stackrel{ˉ}{v}}}^m=\frac{\exp (\mathrm{i}{k}_x{h}_x)-\beta -1}{\exp (\mathrm{i}{k}_x{h}_x)+\beta -1}{\widetilde{\stackrel{ˉ}{v}}}^{m-1}=G(k_x){\widetilde{\stackrel{ˉ}{v}}}^{m-1}.$$

$${\left|G(k_x)\right|}^2=\frac{{\left(\beta +1-\cos \left(k_x{h}_x\right)\right)}^2+{\sin }^2\left(k_x{h}_x\right)}{{\left(\beta -1+\cos \left(k_x{h}_x\right)\right)}^2+{\sin }^2\left(k_x{h}_x\right)}=\frac{A}{B}.$$

$$\frac{c\left({\bar{v}}_{j+1}^m-{\bar{v}}_j^m\right)}{h_x}+\frac{\eta }{2}{\bar{v}}_{j+1}^m=\frac{c\left({\bar{v}}_{j+1}^{m-1}-{\bar{v}}_j^{m-1}\right)}{h_x}-\frac{\eta }{2}{\bar{v}}_{j+1}^{m-1}.$$

$$\frac{c\left({\bar{v}}_{j+1}^m-{\bar{v}}_j^m\right)}{h_x}+\frac{\eta }{4}\left({\bar{v}}_{j+1}^m+{\bar{v}}_j^m\right)=\frac{c\left({\bar{v}}_{j+1}^{m-1}-{\bar{v}}_j^{m-1}\right)}{h_x}-\frac{\eta }{4}\left({\bar{v}}_{j+1}^{m-1}+{\bar{v}}_j^{m-1}\right).$$

$$\begin{array}[]{ll}v(0,t)=f(t),&t\geq 0,\\ v(x,0)=0,&x\geq 0,\\ f(0)=0.\end{array}$$

$${\bar{v}}^m=\frac{1}{3}\left(4{\bar{v}}^m({\Omega }_2)-{\bar{v}}^m({\Omega }_1)\right).$$

$$\begin{array}[]{ll}\displaystyle c\frac{\bar{v}^{m}_{i+1}-\bar{v}^{m}_{i}}{h_{x}}=-\frac{1}{720}\sum_{j=-3}^{1}\alpha_{j}\left(\frac{\eta}{2}\bar{v}^{m}_{i+j}+\Phi_{1}\left(\bar{v}^{m}_{i+j}\right)\right),\end{array}$$

$${\Phi }_1\left({\bar{v}}_i^m\right)=-\frac{\eta }{2}{\bar{v}}_i^{m-1}+c\frac{-{\bar{v}}_{i+2}^{m-1}+8{\bar{v}}_{i+1}^{m-1}-8{\bar{v}}_{i-1}^{m-1}+{\bar{v}}_{i-2}^{m-1}}{12h_x}+O\left(h_x^4\right).$$

$${\bar{v}}^m(x)=\sum _{j=0}^{\infty }{V}_j^m{l}_j(\kappa x),m=0,1,\mathrm{2...},$$

$$\left(\eta +c\kappa \right)V_0^m=\left(-\eta +c\kappa \right)V_0^{m-1}+2c\sqrt{\kappa }\left({\bar{f}}^m-{\bar{f}}^{m-1}\right),$$

$$\left(\eta +c\kappa \right)V_j^m+2c{\rm Y}(V_j^m)=\left(-\eta +c\kappa \right)V_j^{m-1}+2c{\rm Y}(V_j^{m-1}),$$

$${\rm Y}\left(V_j^m\right)=\kappa \sum _{i=0}^{j-1}{V}_i^m=\kappa V_{j-1}^m+{\rm Y}\left(V_{j-1}^m\right),$$

$$V_j^m\equiv 0,{\bar{f}}^m\equiv 0,\forall m<0.$$

$$\left(\eta +c\kappa \right)V_j^m+\left(\eta -c\kappa \right)V_j^{m-1}=\left(\eta -c\kappa \right)V_{j-1}^m+\left(\eta +c\kappa \right)V_{j-1}^{m-1},m=0,1,...;j=1,2,...$$

$$\left\{\begin{array}[]{ll}V^{m}_{0}=\kappa^{-1/2}\left(\bar{f}^{m}-\bar{f}^{m-1}\right),&m=0,1,...,\\ V^{m}_{j}=V^{m-1}_{j-1},&m=0,1,...;\;j=1,2,...\end{array}\right.$$

$${\bar{v}}^m(x)=\sum _{j=0}^{\infty }{V}_j^m{l}_j(\kappa x)=\sum _{j=0}^m{V}_0^{m-j}{l}_j(\kappa x),m=0,1,2,...$$

$$f(t)=\exp \left[-\frac{(2\pi f_0(t-t_0){)}^2}{{\delta }^2}\right]\sin (2\pi f_0(t-t_0)),$$

$$K(x)={\int }_0^{\infty }{v}^2(x,t)dt=\sum _{k=0}^{\infty }{\left[{\bar{v}}^k(x)\right]}^2,$$

$$\frac{\partial u}{\partial t}+c\frac{\partial u}{\partial z}-\sum _{s=1}^n\frac{\partial {\psi }_s}{\partial t}=0,$$

$$\frac{1}{c^2}\frac{{\partial }^2{\psi }_s}{\partial t^2}-{\gamma }_s\frac{{\partial }^2{\psi }_s}{\partial x^2}-{\beta }_s\frac{{\partial }^{2}u}{\partial x^2}=0,s=1,2,...n,$$

$$L\left\{\frac{d^2}{d{t}^2}g(t)\right\}={\left(\frac{\eta }{2}\right)}^2{\bar{g}}_m+{\Phi }_2({\bar{g}}_m),{\Phi }_2({\bar{g}}_m)\equiv {\eta }^2\sum _{j=0}^{m-1}(m-j){\bar{g}}_j.$$

Full text

1 Introduction

Mathematical models based on the one-way wave equation (OWWE) are often considered in problems of ocean acoustics [1, 2, 3], seismic prospecting [4, 5, 6, 7] , as well as for setting non-reflecting boundary conditions [8, 9, 10]

[TABLE]

where $\mathrm{i}=\sqrt{-1}$, $\tilde{u}\equiv\tilde{u}(k_{x},z,\omega)$ is a wave component at the angular frequency $\omega$, $k_{x}$ is the horizontal wave number, $c$ is the wave velocity, the vertical direction $z$ is the extrapolation direction, i.e., the direction of one-way propagation, and the positive axis $z$ is directed downward, i.e., toward increasing depth. The square-root operator can be formally represented by the Padé expansion [11, 12, 13, 14]

[TABLE]

where the coefficients $\gamma_{s},\beta_{s}$ for the propagation angle should be optimized[12, 15]. 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 one-way model is admissible as the correct account of amplitudes increases computer costs [16].

The fundamental problem of the downward continuation algorithms of wave fields is the instability. If the coefficients $\gamma_{s},\beta_{s}$ are real, then for the angles around $\pi/2$ the argument of the square root becomes lesser than zero, the left-hand side of equation (2) being complex, while the right-hand side is still real, hence causing inconsistency in the approximation. This results in an improper propagation of the evanescent mode which should exponentially decay. For stabilizing the real Padé approximation there are a few approaches [17, 18, 19] that allow suppressing unstable components of a wave field. On the other hand, setting the coefficients $\gamma_{s},\beta_{s}$ to be complex [20, 21, 22, 23], a better consistency of the right-hand and the left-hand sides of equation (2) can be attained. From the physical viewpoint this means the introduction of artificial dissipation that restricts an increase in instability for evanescent waves. However the presence of strong gradients of the function of velocity of a medium, the use of the Marchuk-Strang type splitting for decreasing computer costs [24, 25, 26] as well as the simulation of high-frequency wave fields, etc. can bring about both the instability of calculation and excess energy dissipation of waves to be simulated. This is explained by the fact that optimal values of the coefficients $\gamma_{s},\beta_{s}$ are selected based on the frozen coefficients principle for a homogeneous medium, while calculations are carried out for inhomogeneous velocity models employing difference approximations and different decompositions for differential equations. Also, it is well to bear in mind that with the same number of terms in series (2) the accuracy of the complex Padé approximation is somewhat lower than that real.

In [27], another approach was proposed that provides the stability and high accuracy of calculations for the real coefficients $\gamma_{s},\beta_{s}$. In this case the solution is sought for as a series in the Laguerre functions, while for increasing the spatial approximation accuracy the Richardson extrapolation[28] is used, which additionally restricts the instability of approximation (2). Earlier the stabilizing properties of the Richardson method for other problems were noted [29, 30]. However the Richardson method is not efficient because with its use it is needed to additionally solve the original problem on a mesh with a doubled number of nodes. We will consider a stabilizing procedure based on the spline-filtration providing the stability of the Adams type multistep schemes of a high accuracy order [31]. This will allow carrying out calculations with the fifth order of accuracy instead of the fourth order thus to some extent decreasing computer costs.

In addition to the instability of the OWWE another fundamental problem is solving poorly conditioned systems of linear algebraic equations (SLAEs) in a frequency domain [32, 33]. As opposed to the Fourier transform, after employing the Laguerre transform the matrix of SLAEs is real and well-conditioned. Coefficients of the Laguerre series expansion have a recurrent dependence, therefore for their calculation it is required to solve the SLAEs with the same matrix and different right-hand sides, for example, with a parallel dichotomy algorithm [34, 35, 36]. Even a higher efficiency of calculations can be attained if one turns from solving difference problems for elliptic high-order operators to a sequence of problems for a second order operator. To this end, having stabilized the Adams schemes instability we will consider a method of the predictor-corrector type of the fifth accuracy order [31].

2 The stability Analysis for a model 1D one-way wave equation

The aspects of stability in constructing a numerical method for solving the 2D OWWE occupy a highly important place. To investigate the stability let us first consider a model problem for the 1D OWWE:

[TABLE]

with the initial condition $v(x,0)=\varphi(x),\;(\varphi(0)=\varphi(1))$ and the periodic boundary condition $v(0,t)=v(1,t)$.

To solve problem (3), let us consider the direct and inverse Laguerre transforms [37] of a function $g(t)\in L_{2}(0,\infty)$

[TABLE]

where $l_{m}(\eta t)\equiv\sqrt{\eta}\exp({-\eta t/2})L_{m}(\eta t)$ are the orthogonal Laguerre functions, $L_{m}(t)$ is the Laguerre polynomial of $m$ degree and $\eta>0$ is the transformation parameter.

Setting $\lim_{t\rightarrow\infty}g(t)=0$, the following relations are valid [37, 38]

[TABLE]

where

[TABLE]

Making use of transformation (4) for equation (3), we obtain

[TABLE]

where the index $m$ denotes number of a term in series (4). Taking into consideration

[TABLE]

for studying the stability of difference schemes, let us turn to another form of equation (7):

[TABLE]

For solving (8b) let us consider the difference scheme of the first order of accuracy

[TABLE]

Substituting the solution in the form $\bar{v}^{m}_{j}=\tilde{\bar{v}}^{m}\exp(\mathrm{i}k_{x}jh_{x})$ into difference equation (9), obtain

[TABLE]

Here $\beta=\eta h_{x}/(2c)$, $G=G(k_{x})$ is called the amplification factor which is a complex function of the wavenumber $k_{x}$. A difference equation will be stable in the Von Neumann sense [39] if $|G(k_{x})|\leq 1$ $\forall$ $k_{x}$. For equation (10) let us estimate the value

[TABLE]

For $c>0$ obtain $A>0,B>0,A-B=4\beta\left(1-\cos(k_{x}h_{x})\right)\geq 0$, hence, $\left|G(k_{x})\right|^{2}\geq 1$ and scheme (9) will be unstable. For $c<0$ it can be shown that $\left|G(k_{x})\right|^{2}\leq 1$, which suggests the stability of the scheme.

Now let us consider another method of the first order of accuracy

[TABLE]

In a similar manner reducing scheme (11) to the form $\tilde{\bar{v}}^{m}=G(k_{x})\tilde{\bar{v}}^{m-1}$, we obtain that for $c>0$ the value $A>0,B>0,A-B=4\,{\beta}\left(\cos\left({k_{x}h_{x}}\right)-1\right)\leq 0$, and hence, $\left|G(k_{x})\right|^{2}\leq 1$. Thus, scheme (11) will be stable for $c>0$ and unstable for $c<0$.

Let us consider the Crank-Nicolson scheme (CN-Scheme) [40] of the second order of accuracy

[TABLE]

It is not difficult to obtain that for the scheme in question $|G(k_{x})|=1$ holds, therefore scheme (12) is unconditionally stable.

For multistep schemes of a high-order of accuracy of the Adams-Moulton type (AM scheme) the stability has been analytically studied. Having omitted cumbersome operations we will restrict ourselves to the main conclusions. The AM-schemes of the third and fourth orders of accuracy will be unstable for periodic boundary conditions if $c>0$ and stable if $c<0$. The schemes of the fifth order of accuracy and higher are absolutely unstable. The Von Neumann spectral feature is a necessary but not sufficient condition of stability [41], and in the case of non-periodic boundary conditions the stability is not preserved for the AM-schemes of higher than the second order of approximation. To overcome these difficulties we will consider the ways of stabilizing the multistep AM-schemes of a high-order of accuracy first for the 1D and then for the 2D OWWE.

3 The stabilization of high-order schemes for the 1D one-way wave equation

For equation (3) at $c>0$ instead of the periodic boundary conditions we consider initial and boundary conditions of the form:

[TABLE]

A stable spectral-difference algorithm for solving the 2D OWWE was proposed in [27]. It includes the Richardson extrapolation procedure that for problem (7), (13) and scheme (12) can be written down in the following form.

Algorithm 1-1. The Richardson extrapolation.

Let the auxiliary functions $\bar{v}^{m}(\Omega_{1})$, $\bar{v}^{m}(\Omega_{2})$ be defined on the meshes $\Omega_{1},\Omega_{2}$ with the mesh steps $h_{x}$ and $h_{x}/2$. To calculate the functions $\bar{v}^{m}$ accurate to $O(h_{x}^{4})$, the following is necessary:

Based on the cubic splines interpolate values of the function $\Phi_{1}(\bar{v}^{m})$, preset on the mesh $\Omega_{1}$, into nodes $\Omega_{2}$. 2. 2.

On the mesh $\Omega_{1}$, applying equation (12), calculate the solution $\bar{v}^{m}(\Omega_{1})$. 3. 3.

On the mesh $\Omega_{2}$, applying equation (12), calculate the solution $\bar{v}^{m}(\Omega_{2})$. 4. 4.

Based on the Richardson extrapolation, correct the mesh function with

[TABLE] 5. 5.

Turn to the calculation of the $(m+1)\text{th},\;\text{the }(m+2)\text{th}$, etc. coefficients of the expansion of the Laguerre series.

This technique of calculating the Laguerre series coefficients is stable and provides the fourth order of accuracy. However the necessity of calculating $\bar{v}^{m}(\Omega_{2})$ triples the common computer costs, therefore there arises a problem of constructing a more efficient method of no less than fourth order of accuracy.

Let us consider a difference approximation for equation (7) based on the Adams-Moulton multistep method of the fifth order of accuracy

[TABLE]

where the coefficients of the difference scheme are equal to $\alpha_{-3}=-19,\;\alpha_{-2}=106,\;\alpha_{-1}=-264,\;\alpha_{0}=646,\;\alpha_{1}=251$. This scheme is unstable according to the Von Neumann spectral property, but it can be stabilized when carrying out calculations in the following manner.

Algorithm 1-2. The stabilization of the Adams-Moulton scheme via the quintic spline filtration.

To calculate the functions $\bar{v}^{m}$ accurate to $O(h_{x}^{5})$ on the mesh $\Omega$ with the mesh step $h_{x}$, the following is necessary:

Let the number of nodes of the mesh $\Omega$ be odd. Construct the quintic splines for the function $\Phi_{1}({\bar{v}}^{m})$ using only odd nodes of the mesh. 2. 2.

Replace values of the function $\Phi_{1}(\bar{v}^{m}_{k})$ for even $k$ by their interpolated values (the quintic spline filtration). 3. 3.

Applying equation (14), calculate the solution $\bar{v}^{m}$. 4. 4.

Turn to the calculation of the $(m+1)\text{th},\;\text{the}(m+2)\text{th}$, etc. coefficients of the expansion of the Laguerre series.

Such an algorithm of calculations makes possible to stabilize the numerical instability of scheme (14) and to attain a higher approximation order as compared to the Richardson extrapolation. To stabilize the solution, instead of the quintic spline one can use other interpolation algorithms [42, 43]: the cubic spline interpolation, barycentric, Lagrangian, etc. However numerous computer-aided experiments have not revealed any advantages over splines because the procedures of constructing splines are efficient enough as compared to solving elliptic equations in the 2D case. In addition, the barycentric interpolation demands high computer costs and both the Lagrangian interpolation and the cubic spline-interpolation are more dissipative than the quintic splines. If splines are not being used, the nodes of an interpolating polynomial should be symmetrically placed regarding the node for which the interpolated value is calculated. Otherwise due to the asymmetry of interpolation nodes the profile of a wave is distorted or the instability of calculation arises.

We can offer another way of stabilizing the numerical instability of the AM-schemes, which does not demand the calculation of splines.

Algorithm 1-3. The stabilization of the Adams-Moulton scheme via inconsistent approximation.

Taking into account the equivalence of problems (7) and (8), for computing the values of the grid functions $\Phi_{1}\left(\bar{v}^{m}_{i}\right)$ instead of (6) use the following approximation

[TABLE] 2. 2.

Solve equation (7) through scheme (14). 3. 3.

Turn to the calculation of the $(m+1)\text{th},\;\text{the }(m+2)\text{th}$, etc. coefficients of the expansion of the Laguerre series.

The non-consistent approximation of the operator $\partial/\partial x$ for the right-and the left-hand sides of equations (7),(8), stipulates supplementary non-physical dissipation preventing the development of instability. However, if instead of the central approximation of the fourth order of accuracy in (15) one uses a higher order approximation or a non-central scheme, the stability is lost.

In addition to the above we have considered a stable algorithm based on the central differences of a high-order of accuracy for the operator $\partial/\partial x$. For the 1D OWWE such calculations were carried out up to the schemes accurate to the twelfth order. As a result, the stability and high-order of accuracy were attained. However this approach will not work for the 2D OWWE because of high computer costs that are required for solving SLAEs.

Thus, in addition to the Richardson method, other stable algorithms of a high-order of accuracy can be proposed. However for calculating of $\bar{v}^{m}$ the function $\bar{v}^{m-1}$ should be known throughout the whole calculation domain as its values are needed for implementing the stabilizing procedures.

4 The analytic solution via the Laguerre transforms for the 1D one-way wave equation

In order to assess the accuracy of the algorithms proposed, let us consider a fully analytical method for solving the 1D OWWE. To satisfy boundary conditions (13), we seek the solution to equation (8) in the form

[TABLE]

where the transformation parameter $\kappa>0$. Then, after applying the Laguerre spatial transform to equation (8) with allowance for initial boundary conditions (13) we have

[TABLE]

where

[TABLE]

[TABLE]

Taking (18) into account, equation (17b) takes the following form

[TABLE]

Since $c>0$, then selecting $\kappa=\eta/c$, we finally obtain

[TABLE]

Based on (16),(20), we can write down the solution to equation (8),(13) in the form

[TABLE]

The latter sum in (21) is a discrete linear convolution, therefore for a given $x$ the functions $\bar{v}^{m}(x),\;m=0,...,M$ can be computed in $O(M\log M)$ arithmetical operations based on the algorithm of the FFT [44]. Note that if we select $\kappa\neq\eta/c$, the solution for $\bar{v}^{m}$ will not be representative in the form of convolution thus increasing computer costs. The final solution to equation (3) in the time domain is calculated via the inverse Laguerre transform (4).

5 Numerical Experiments for the 1D one-way wave equation

For testing the methods proposed for the 1D OWWE, we used a homogeneous medium model with the speed $3000\;m/s$ and the size $7.5\;km$. For the test calculations we set boundary condition (13) depending on time as

[TABLE]

where $t_{0}=0.2s,\,\delta=4,\,f_{0}=30\text{Hz}$. As compared to the Fourier transform, where the basis functions are uniquely defined, the parameter $\eta$ for using the Laguerre transform (4) should be set. This parameter was experimentally chosen based on the analysis of the convergence rate of the Fourier-Laguerre series for the shifted function $f(t)$ with $t_{0}=T$, where T is the upper boundary of the time interval for which the wave field is calculated. The parameter $\eta$ is chosen such that the function f(t) with $t_{0}=T$ in the mean-quadratic norm is approximated accurate to $\varepsilon<10^{-10}$. The number of addends in series (4) was $n=2500$ for $T=2\;s$; the expansion parameters were $\eta=600$.

From Table 1 it is evident that with decreasing the mesh size by a factor of two, the error of the Adams and the Richardson methods is decreasing according to theoretical approximation order. For example, for meshes with the number of nodes $N_{x}=2000$ and $N_{x}=4000$ the values of the error of the AM-scheme of fifth order with the quintic spline interpolation (AM5-I5) has $31$ times difference, which almost corresponds to the fifth order of approximation. The AM-scheme of the sixth order with the seventh order spline interpolation (AM6-I7) demonstrates the sixth order of approximation, while the AM-scheme of the fifth order with formula (15) (AM5-D4) is only of the fourth order of approximation.

From Table 1 and Fig. 1 it also follows that with an equal mesh size, the Richardson extrapolation of the fourth order of approximation is more accurate as compared to the Adams methods of the fifth and sixth orders. There is no contradiction because the estimation of the accuracy of difference schemes includes a constant independent of the mesh size. This constant is smaller for the Richardson method as compared to the AM5-15 method as all splines are constructed on a mesh with a doubled step. In addition, the Richardson method requires the solution to a supplementary problem on the mesh $\Omega_{2}$, therefore it is more correct to compare the accuracy of calculations when the general number of nodes of the meshes $\Omega_{1},\Omega_{2}$ is equal to the number of nodes of the mesh $\Omega$ for other methods. Indeed, comparing solutions for the Richardson method with $N_{x}=1500$ for the mesh $\Omega_{1}$ and the Adams methods with $N_{x}=4500$, it is evident (Table  1) that the latter are significantly more accurate. Applying the spline-filtration procedure does not bring about a considerable loss in accuracy, otherwise the accuracy of the methods AM5-I5, AM6-I7 would be lower or compatible with the fourth order method AM5-D4, for which the filtration is not employed.

For comparison similar calculations were carried out for the CN-Scheme and the explicit Runge-Kutta method (the RK4-method) [31], which are of the second and fourth orders of accuracy, respectively. For the method RK4, values of the grid function $\Phi_{1}\left(\bar{v}^{m}\right)$ in semi-integer nodes of the mesh were calculated based on the quintic splines. The implementation of the implicit Runge-Kutta method of a high-order is not reasonable because of essentially higher computer costs as compared to the approaches proposed in the given study. From Table 1 and Figure 1 it is clear that the method RK4 and the CN-scheme are stable and converge to the analytical solution with decreasing a mesh size. However for large mesh steps, the method RK4 possesses a pronounced numerical dissipation, while on the contrary, the Crank-Nicolson scheme demonstrates a dispersive error. A low accuracy for large mesh steps makes the application of these algorithms disadvantageous in comparison with the Adams methods.

To evaluate dissipative properties of the methods proposed, let us consider the integral of the form

[TABLE]

where the latter equality is the Parseval relation.

For problem (3),(13) at $c=const>0$ with a sufficient number of terms in Laguerre series (4), $K(x)=const$ should hold with a good accuracy. Figure 2a,b shows that for the analytical method (20),(21) and for the CN-scheme the value $K(x)$ is preserved with the precision of a machine. However due to the numerical dispersion the solution obtained with the CN-scheme does not satisfy the 1D OWWE equation with some kind of accuracy for large mesh steps. The Richardson method at $N_{x}=1000$ (Fig. 2a) is less dissipative than the algorithms AM5-I5 and AM6-I7, while with an increase of the number of mesh nodes (Fig. 2b) the situation is contrary. This means that when solving the 2D OWWE by the Richardson method the stability will be stronger as compared to the Adams methods. The method RK4 is most dissipative among all under consideration, hence the initial impulse from a source has smoothed into the straight line (Fig. 1a). Thus, the explicit method RK4 and the CN-scheme, as was already mentioned, cannot be offered for the use within the Laguerre method.

6 The 2D one-way wave equation solver

6.1 Temporal approximation

Let us write down 2D OWWE (1),(2) for the spatial-temporal domain [13, 14]

[TABLE]

where $u\equiv u(x,z,t)$ is the field variable, $\psi_{s}\equiv\psi_{s}(x,z,t)$ are auxiliary functions.

Assuming $g(0)=\left.\frac{dg}{dt}(t)\right|_{t=0}=0$ and $\lim_{t\rightarrow\infty}g(t)=\lim_{t\rightarrow\infty}\frac{dg}{dt}(t)=0$, we can show [37, 38] that

[TABLE]

Then applying the Laguerre transform (4) to equations (24) we obtain the following system of equations for the calculation of the $m$-th coefficient of expansion:

[TABLE]

where $\tilde{\eta}=\eta/2$ and the index $m$ denotes number of a term in series (4). The polynomial coefficients $\gamma_{s},\beta_{s}$ for $n=3$ are chosen as follows:  $\gamma_{1}=0.972926132$, $\gamma_{2}=0.744418059$, $\gamma_{3}=0.150843924$, $\beta_{1}=0.004210420$, $\beta_{2}=0.081312882$, $\beta_{3}=0.414236605$, for which, as shown in [5, 15], such approximation is valid up to the angles of $89$ degrees.

6.2 The spatial approximation

Multistep schemes of the Adams type of a high-order for solving the 2D OWWE are of practical importance as opposed to the Richardson method which requires solving an auxiliary problem on the mesh $\Omega_{2}$ with a doubled number of nodes. To approximate equation (25) we will use the Adams-Moulton scheme of the fifth order of accuracy:

[TABLE]

where the coefficients of a difference scheme $\alpha_{-3}=-19,\;\alpha_{-2}=106,\;\alpha_{-1}=-264,\;\alpha_{0}=646,\;\alpha_{1}=251$ and the difference operator $\mathcal{L}_{x}$ is of the form

[TABLE]

For approximating $\partial^{2}/\partial x^{2}$ it is reasonable to use the dispersion-relationship-preserving method (DRP) by Tam and Webb [45], for which according to the Fourier derivative rule, $k_{j}\Longleftrightarrow-\mathrm{i}\partial_{j},$ values of the optimized coefficients $a_{n}$ in (27) are defined as solution to the problem of minimizing the error functional in the space of wave numbers. This approach and its various modifications [46, 47, 48] 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 [42]. To provide the twelfth approximation order the coefficients of difference scheme (27) were chosen as follows [48]: $a_{0}=-3.12513824$, $a_{1}=1.84108651$, $a_{2}=-0.35706478$, $a_{3}=0.10185626$, $a_{4}=-0.02924772$, $a_{5}=0.00696837$, $a_{6}=-0.00102952$.

6.3 The solution of the SLAEs

Let us write down the difference problem (26) in the form of a SLAE as

[TABLE]

where $\bar{\mathbf{\Theta}}_{k}^{m}=-\bar{\mathbf{U}}_{k}^{m}+\sum_{s=1}^{3}\bar{\mathbf{\Psi}}_{k}^{m,s}$ and $I$ is the unit matrix. Employing the Schur complement [49], the mesh functions $\bar{\mathbf{U}}_{k+1}^{m}$ can be calculated through the solution to the following reduced SLAE

[TABLE]

where

[TABLE]

Making use of the matrix property [50] for (29)

[TABLE]

multiplying the equation by the matrix $M_{1}M_{2}M_{3}$ and taking into consideration the commutative property of $M_{i}M_{j}=M_{j}M_{i}$, we obtain the governing equation for the calculation of the mesh functions $\mathbf{\bar{U}}^{m}_{k+1}$

[TABLE]

As opposed to the Fourier transform, the coefficients of the Laguerre expansion in series (4) are dependent in a recurrent manner (25). Hence, for a fixed $k$ for different $m$ it is required to solve SLAEs many times with a common real matrix and different right-hand sides. Matrix (31) is banded and can be explicitly represented without calculation of the matrices $M_{s}^{-1}$ thus allowing us to apply efficient algorithms for solving SLAEs based on $LU$-decomposition. Solving the SLAEs with banded matrices with a parallel algorithm, it is reasonable to use the parallel dichotomy algorithm [34, 36, 35], which was developed for tridiagonal matrices and block-tridiagonal matrices. With respect to the number of arithmetical operations, the dichotomy algorithm is comparable with other available algorithms; however the time needed for inter-process 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 inter-processor interaction optimization algorithms.

After the calculation of the mesh functions $\bar{\mathbf{\mathbf{U}}}^{m}_{k+1}$ before turning to calculating the functions $\bar{\mathbf{\mathbf{U}}}^{m}_{k+2}$, the functions $\bar{\Psi}^{m,s}_{k+1}$ should be calculated as

[TABLE]

Making use of the property (30), we arrive at

[TABLE]

For solving the 2D OWWE one needs not only the stabilization of the numerical instability of difference approximation for the operator $\partial/\partial z$, but also the instability of the real Padé approximation (2). The method AM5-I5 allows solving both these problems.

Algorithm 2-1. The Adams-Moulton downward-continuation procedure for the 2D OWWE.

To calculate the mesh functions $\bar{\mathbf{U}}^{m},\bar{\mathbf{\Psi}}^{m}$ accurate to $O(h_{x}^{\xi}+h_{z}^{5})$, the following is necessary:

Let the number of nodes of the mesh $\Omega$ in the direction $z$ be odd. For all $i$ for the functions $\Phi_{1}(\bar{u}^{m}_{ik})$, $\Phi_{2}(\bar{\psi}^{m,s}_{ik})$ construct, independently, the 1D quintic splines in the direction $z$, using only odd values of $k$. 2. 2.

Replace the values of the functions $\Phi_{1}(\bar{u}^{m}_{ik})$, $\Phi_{2}(\bar{\psi}^{m,s}_{ik})$ for even $k$ by their interpolated values (the quintic spline filtration). 3. 3.

For $k=4,...,K-1$

1. 3.1.

Applying equation (31), calculate the solution $\bar{\mathbf{U}}^{m}_{k+1}$. 2. 3.2.

Applying equation (32), calculate the solution $\bar{\mathbf{\Psi}}^{m,s}_{k+1},\;s=1,2,3$ 4. 4.

Turn to the calculation of the $(m+1)\text{th},\;\text{the}(m+2)\text{th}$, etc. coefficients of the expansion of the Laguerre series.

The above-considered way of stabilizing the solution allows the stability not only of the Adams-Moulton implicit schemes, but also the Adams-Bashfort explicit schemes that are stable for essentially lesser steps $h_{z}$. As consequence, the number of SLAEs to be solved of the form (31),(32) is multiply increased thus making the Adams-Bashfort method inefficient. Let us now consider the predictor-corrector method combining the computational efficiency of explicit and high stability of the implicit Adams schemes.

For equation (25a) as a predictor procedure we choose the Adams-Bashfort explicit method of the fifth order of accuracy

[TABLE]

where $\varrho_{0}=1901,\;\varrho_{-1}=-2774,\;\varrho_{-2}=2616,\;\varrho_{-3}=-1274,\;\varrho_{-4}=251$. In terms of correction we will use the Adams-Moulton scheme of the fifth order of accuracy (26a), where unknown values $\bar{\psi}^{m,s}_{ik+1}$ should be replaced by predicted values. In this case, the functions $\bar{u}^{m}_{ik+1}$ are explicitly expressed. If for correction we use scheme (26a) and substitute the predicted values both for $\bar{\psi}^{m,s}_{ik+1}$ , and for $\bar{u}^{m}_{ik+1}$ into the right-hand side, then for the sake of stability, smaller steps $h_{z}$ and simultaneously a larger number of corrections will be required.

Algorithm 2-2. The Predictor-Corrector downward-continuation procedure for the 2D OWWE.

To calculate the mesh functions $\bar{\mathbf{U}}^{m},\bar{\mathbf{\Psi}}^{m}$ accurate to $O(h_{x}^{\xi}+h_{z}^{5})$, the following is necessary:

Let the number of nodes of the mesh $\Omega$ in the direction $z$ be odd. For all $i$ for the functions $\Phi_{1}(\bar{u}^{m}_{ik})$, $\Phi_{2}(\bar{\psi}^{m,s}_{ik})$ independently construct the 1D quintic splines in the direction $z$ using only odd values of $k$. 2. 2.

Replace values of the functions $\Phi_{1}(\bar{u}^{m}_{ik})$, $\Phi_{2}(\bar{\psi}^{m,s}_{ik})$ for even $k$ by their interpolated values (the quintic spline filtration). 3. 3.

For $k=4,...,K-1$

1. 3.1.

Applying equation (34), calculate the predicted solution ${\breve{\mathbf{U}}}^{m}_{k+1}$. 2. 3.2.

Applying equation (33) with ${\breve{\mathbf{U}}}^{m}_{k+1}$, calculate the predicted solution $\breve{\Psi}^{m,s}_{k+1},\;s=1,2,3$ . 3. 3.3.

Applying equation (26a), substituting $\breve{\Psi}^{m,s}_{k+1}$ instead of $\bar{\Psi}^{m,s}_{k+1}$, calculate the corrected solution $\breve{\breve{\mathbf{U}}}^{m}_{k+1}$. 4. 3.4.

Applying equation (33) with $\breve{\breve{\mathbf{U}}}^{m}_{k+1}$, calculate the final solution for $\bar{\mathbf{\Psi}}^{m,s}_{k+1},\;s=1,2,3$ . 5. 3.5.

Applying equation (26a) with $\bar{\mathbf{\Psi}}^{m,s}_{k+1}$, calculate the final solution for $\bar{{\mathbf{U}}}^{m}_{k+1}$. 4. 4.

Turn to the calculation of the $(m+1)\text{th},\;\text{the }(m+2)\text{th}$, etc. coefficients of the expansion of the Laguerre series.

Thus, within the predictor-corrector method instead of indefinite non-symmetrical SLAE (31) it is necessary to solve SLAEs of the form (33) with sign-defined symmetric matrices of lesser dimensions that allows the use of efficient algorithms of the computational linear algebra and a decrease in the calculation time. As compared to the Marchuk-Strang method that is accurate to the second order, the predictor-corrector method is of the fifth order of accuracy. Further increase of approximation order is not reasonable because for providing the stability an essentially smaller step $h_{z}$ should be set.

Both for the Adams method and for the predictor-corrector method one needs initial values to start the calculation. It is required to use other methods such as the Richardson extrapolation or Crank-Nicolson scheme with a smaller step for obtaining such initial values.

7 Numerical experiments for the 2D one-way wave equation

Analytically it is really difficult to provide a strict mathematical substantiation of algorithms of a high-order of accuracy for the 2D OWWE. Therefore to confirm the efficiency of the approaches proposed, thorough testing is needed. Let us discuss a few tests that would allow the evaluation of the quality of the solution to be obtained as compared with the Richardson algorithm. Numerical procedures were implemented in Fortran-90 using the MPI library.

7.1 The impulse response

In the first test we illustrate analyzing the accuracy by the impulse responses. For the calculation, we used the homogeneous medium model with the speed $250\;m/s$ and the size $3.5\;km\times 1.5\;km$. The point source (22) with the parameters $t_{0}=0.2\,s,\,\delta=4,\,f_{0}=30\,\text{Hz}$ was located at the center of the upper surface. The number of addends in series (4) was $n=4000$ for $T=6\;s$; the expansion parameter was $\eta=600$.

A disadvantage of the Laguerre transform is the absence of the fast transformation algorithm. The implementation of the forward Laguerre transform (4) with the help of the method of least squares, the number of arithmetical operations will be of order $O(NP)$, where $P$ is the number of discrete points of the approximated function and $N\leq P$ is the number of terms in the Laguerre series needed for attaining the required accuracy in the norm $L_{2}$, whereas for the fast Fourier transform computer costs are essentially less and make up $O(N\log(N))$. However, taking into account the fact that input data are set only along the upper surface ($z=0$), 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 series (4) from the solution to problem (25). The numerical experiments have confirmed that the time needed for carrying out the Laguerre transform for the initial data is less than one percent of the total calculation time.

In Section 5, for the 1D OWWE it was shown that with an equal mesh size the Richardson extrapolation is more accurate than the AM5-I5-scheme, which is also valid for the 2D case. If one selects the step $h_{z}=1$ m for the mesh $\Omega_{1}$ for the Richardson method and the step $h_{z}=0.3$ m for the AM5-15 and PC5-15 methods, then the whole volume of calculations and calculation accuracy for all the three methods will be compatible (Fig. 3). However with a thorough consideration of values of the amplitudes along the straight line ”Slice”, it appears to be clear (Fig. 4) that the AM5-I5 and the PC5-I5 methods are more accurate when the number of nodes of the mesh $\Omega$ is equal to the number of nodes of the meshes $\Omega_{1},\Omega_{2}$.

A considerable accuracy and computational efficiency of the AM5-I5, PC5-I5 methods is attained at the expense of lesser stability as compared to the Richardson method. It was experimentally revealed that for the Richardson method the condition of stability is of the form $h_{z}/h_{x}\leq 1$, while for the AM5-I5 method it should be $h_{z}/h_{x}<0.4$, and for the PC5-I5 algorithm the stability is attained at $h_{z}/h_{x}<0.3$. The fact that minimum steps required for providing a good accuracy and stability almost coincide, allows us to propose the statement about the balance of the PC5-I5 method which is by $30-40\%$ more efficient as compared to the AM5-I5 method.

Additionally we have considered the AM6-I7 and the AM5-D4 methods, which are stable for the 1D OWWE but unstable for the 2D OWWE. This is because in addition to the numerical instability due to the choice of approximation for the operator $\partial/\partial z$, there is instability caused by the presence of a singular component in the solution to the OWWE, i.e. when denominators (2) are close to zero or vanish. In the AM5-D4 method, to approximate $\partial/\partial z$ different difference schemes were used for the right-hand and the left-hand sides of equations (7),(8) thus additing numerical dissipation, but this approach does not allow restricting the growth of the number of singular components for the 2D OWWE. Also, the AM6-I7 method demonstrates (Fig. 2) lesser dissipation as compared to the algorithm AM5-I5 and, hence, insufficient level of fictitious absorption does not allow stabilizing the instability for the 2D OWWE.

7.2 Migration procedures

Testing the algorithms based on the combination of the Laguerre transform and the Richardson extrapolation was done in [27], where on solving a problem of seismic prospecting, as an example, the post-stack migration procedure within the explosive boundaries model has been implemented [4]. As compared to the Finite Difference (FD) [11], Fourier Finite Difference (FFD) [17] and Phase Shift Plus Interpolation (PSPI) [51] methods the algorithm based on the Laguerre transform made possible to obtain a more qualitative solution. Let us now consider similar tests for the algorithms AM5-I5, PC5-I5 for solving the 2D OWWE.

7.2.1 The syncline model

Theoretical seismograms (Fig. 5b) for the syncline model (Fig. 5a) were obtained with the help of the Gaussian beams algorithm [52, 53] implemented in the package Seismic Unix. For setting the boundary condition on the upper surface, the function for the zero-offset section $\left.u(x,z,t)\right|_{z=0}=g(x,t)$ was expanded in series (4) with the parameters $n=2500$ and $\eta=800$ for $t\in[0,4]\;s$. The calculations were carried out on meshes with the steps $h_{x}=6\;m$, $h_{z}=3$ or $6$ m for the Richardson method $h_{x}=6\;m$ and $h_{z}=3$ or $1$ m for the methods AM5-I5, PC5-I5. According to the model of explosive boundaries, the calculation velocities were set to be half the true velocity of the medium model.

The possibility to steadily calculate the solution to the 2D OWWE for non-smooth velocity models is the advantage of the Richardson method while preliminary smoothing for providing stability is needed for the algorithms AM5-I5, PC5-I5, because application of methods of a high-order is based on the assumption that the solution and the velocity function possess the required smoothness. In order to restrict an increase in numerical instability with a minimum influence on the wave kinematics we use averaging of the form

[TABLE]

However with decreasing the step $h_{x}$ the growth of singular components of the field in the area of discontinuities of the velocity function increases, therefore a multiple application of formula (35) may be needed.

In Fig. 6a,b it is clear that the Richardson extrapolation is less accurate if the general number of nodes of the meshes $\Omega_{1},\Omega_{2}$ is equal to the number of nodes of the mesh $\Omega$ for the methods AM5-I5, PC5-I5. If the mesh steps are equal, the quality of the image obtained is approximately the same, but the Richardson method requires three times as many calculations as compared to other techniques. The algorithms AM5-I5 (Fig.6c,d) and PC5-I5(Fig. 6e,f) allow obtaining images that are of the same accuracy, which witnesses to the correctness of calculations because different ideas of using the Adams schemes underlie these methods. The spline-filtration procedure makes possible to provide the stability with lesser costs than the Richardson extrapolation and the method PC5-I5 is three times more efficient than the method AM5-I5.

7.2.2 The Sigsbee model

For Sigsbee2A model [54] (Fig. 7a), theoretical seismograms (Fig. 7b) were calculated using the algorithm of explosive boundaries implemented in the Madagascar package [55]. For setting the boundary condition, the function for zero-offset section$\left.u(x,z,t)\right|_{z=0}=g(x,t)$ was expanded in series (4) for $n=3500$, the parameter $\eta=300$ for $t\in[0,12]\;s$. The calculations were done on the meshes with the steps $h_{x,z}=12.5\;m$ for the Richardson method and $h_{x}=12.5\;m$, $h_{z}=4.16\;m$ for the methods AM5-I5, PC5-I5.

The velocity model Sigsbee is not smooth, therefore preliminary smoothing (35) for the stability of the algorithms AM5-I5 and PC5-I5 is required. If for the model Syncline a single smoothing to provide the stability is needed, for the model Sigsbee the smoothing procedure was applied three times. For the Richardson method the smoothing procedure of the velocity model because of a higher inner dissipation and stability, was not carried out.

The simulation results for the Richardson method are presented in Fig. 8; the results obtained for the methods AM5-I5, PC5-I5 are given in Figure 9, because for these methods there is no difference in solution. The Richardson method with the total number of nodes of the meshes $\Omega_{1},\Omega_{2}$ equal to the number of nodes of the mesh $\Omega$ for the methods AM5-I5, PC5-I5 possesses less accuracy for amplitudes, which is due to a stronger numerical dissipation. All the three algorithms have demonstrated the stability, however the multistep procedures of the predictor-corrector type are more efficient.

8 Conclusion

As a result of the conducted study it appeared possible to determine the reasons of the numerical instability for the schemes of a high-order of accuracy in solving the OWWE and to propose the stabilizing procedure based on the spline-filtration. This allowed us to implement the Adams multistep schemes and the predictor-corrector method accurate to the fifth order. A combination of the Laguerre transform and the predictor-corrector method reduces the original problem to solving the SLAEs with symmetrical well-conditioned real matrices, which is one of the main advantages of the methods proposed as compared to the classical Fourier approach. In the 2D case, direct algorithms for solving SLAEs are applied, which are less sensitive to the above-mentioned properties of matrices. But for the 3D problems, conditioning of matrices is defined by the rate of convergence of iterative procedures and hence the total calculation time.

In spite of the fact that the Richardson extrapolation procedure demands is more computationally expensive costs than the Adams multistep methods, one should not completely reject its application. First, the Richardson method possesses a greater numerical dissipation, which in many cases makes possible to calculate inhomogeneous velocity models without preliminary smoothing. Also, supplementary stability will not be redundant when considering the OWWE for an elastic model. Second, the Richardson method can be used to calculate initial values for multistep methods that are not self-starting. As a rule, for solving this problem the Runge-Kutta type schemes are employed, but within the Laguerre method such schemes do not provide the required accuracy due to a strong numerical dissipation.

A combination of the spline-filtration, the Adams multistep methods and the Laguerre transform is mutually complementary. Experimentally it was verified that the change of the Adams methods for the backward difference schemes does not provide the stability of calculation with the help of proposed stabilizing procedures, whereas the change of the Laguerre transform for the Fourier transform with respect to time makes the spline-filtration unreasonable. This is because in this case the solution for each harmonic is independently determined by initial conditions on the daily surface. On the contrary, the matter of coefficients of the Laguerre series and their recurrent dependence make possible to delicately remove unstable components of the wave field not additing new numerical artifacts. Thus, the considered ways of decreasing computer costs make the proposed methods of solving the OWWE to be promising for the calculation of applied problems.

9 Acknowledgments

The calculation was performed on the supercomputers of the Moscow and the Novosibirsk State Universities and of the Siberian supercomputer center. The work was supported by a grant from the Russian Ministry of Education no. MK-152.2017.5.