key: cord-027118-2xm8nkmi
authors: Sevastianov, Leonid A.; Lovetskiy, Konstantin P.; Kulyabov, Dmitry S.
title: An Effective Stable Numerical Method for Integrating Highly Oscillating Functions with a Linear Phase
date: 2020-06-15
journal: Computational Science - ICCS 2020
DOI: 10.1007/978-3-030-50417-5_3
sha: 
doc_id: 27118
cord_uid: 2xm8nkmi

A practical and simple stable method for calculating Fourier integrals is proposed, effective both at low and at high frequencies. An approach based on the fruitful idea of Levin, which allows the use of the collocation method to approximate the slowly oscillating part of the antiderivative of the desired integral, allows reducing the calculation of the integral of a rapidly oscillating function (with a linear phase) to solving a system of linear algebraic equations with a triangular or Hermitian matrix. The choice of Gauss-Lobatto grid nodes as collocation points let to increasing the efficiency of the numerical algorithm for solving the problem. To avoid possible numerical instability of the algorithm, we proceed to the solution of a normal system of linear algebraic equations.

The initial formulation of the method of numerical integration of highly oscillating functions by Levin and his followers suggests a possible ambiguity in finding the antiderivative: any solution to the differential equation without boundary (initial) conditions can be used to calculate the desired value of the integral.

Levin's approach [1] to the integration of highly oscillating functions consists in the transition to the calculation of the antiderivative function from the integrand using the collocation procedure in physical space. In this case, the elements of the degenerate [2] differentiation matrix of the collocation method [3] are a function of the coordinates of the grid points, the matrix elements are calculated using very simple formulas. In books [3, 4] various options for the implementation of this method are considered, many applied problems are solved.

The method proposed by Levin both in the one-dimensional and in the multidimensional case was published by him in articles [1, 9] , and then he was thoroughly studied in [10] . The method is presented in great detail in the famous monograph [4] , which describes the evolution of numerical methods for integrating rapidly oscillating functions over the past fifteen years.

There are a large number of works using various approaches in order to propose fast and stable methods for solving systems of linear algebraic equations (SLAE) that arise when implementing the collocation method. However, many of them [5] [6] [7] encounter difficulties in solving the corresponding systems of linear equations.

In particular, the use in specific implementations of the Levin collocation method in the physical space of degenerate Chebyshev differentiation matrices, which also have eigenvalues differing by orders of magnitude, makes it impossible to construct a stable numerical algorithm for solving the resulting SLAEs. The approach to solving the differential equation of the Levin method, described in [5, 6, 8] , is based on the approximation of the solution, as well as the integrand phase and amplitude functions in the form of expansion into finite series in Chebyshev polynomials. Moreover, to improve the properties of the algorithms, and hence the matrices of the corresponding SLAEs, three-term recurrence relations are used that connect the values of Chebyshev polynomials of close orders. However, these improvements are not enough to ensure stable calculation of integrals with large matrix dimensions.

In our work, we consider a method of constructing a primitive, based on the spectral representation of the desired function.

We propose increasing the efficiency of the algorithm by reducing the corresponding system of linear equations to a form that is always successfully solved using the LU-decomposition method with partial selection of the leading element.

Consider the integral that often occurs in Fourier analysis -in applications related to signal processing, digital images, cryptography and many other areas of science and technology.

In accordance with the Levin method, the calculation of this integral reduces to solving an ordinary differential equation

As argued in [1] , the system (2) has a particular solution which is not rapidly oscillatory, and we shall look for an approximation to this particular solution by collocation with 'nice' functions, e.g. polynomials. If the unknown function p(x) is a solution of Eq. (2), then the result of integration can be obtained according to the formula

Below we will consider the special case of integration of a highly oscillating function with a linear phase, reduced to the standard form

This can be justified, in particular, by the fact that in many well-known publications [7, 11, 12] stable transformations are discussed in detail, which make it possible to proceed from a general integral with a nonlinear phase to an integral in standard form (on the interval [−1, 1]) with a linear phase.

In the paper by Levin [1] , to automatically exclude the rapidly oscillating component ce Àixx of the general solution p x ð Þ ¼ p 0 x ð Þ þ ce Àixg x ð Þ , it is proposed to search for a numerical solution (2) based on the collocation method, using its expansion in a basis of slowly oscillating functions, rather than using difference schemes (or methods of the Runge-Kutta type).

In this case, the following statement is true [2] :

Statement. The solution of Eq. (2) obtained using the Levin collocation method is a slowly oscillating function O x À1 ð Þ for x ) 1.

Let us consider in more detail the problem of finding the antiderivative integrand, or rather, the approximating polynomial p x ð Þ, satisfying condition (2) in a given number of points on the interval [−1, 1]. Consider the spectral method of finding an approximating function in the form of expansion in a finite series

in the basis of Chebyshev polynomials of the first kind T k x ð Þ f g 1 k¼0 , defined in the Hilbert space of functions on the interval [−1, 1].

The application of the collocation method to solve the problem p 0 x ð Þ þ ixp x ð Þ ¼ f x ð Þ leads to the need to fulfill the following equalities for the desired coefficients

at the collocation points x 0 ; x 1 ; . . .; x n f g . The last statement is equivalent to the fact that the coefficients c k ; k ¼ 0; . . .; n should be a solution to the system of linear algebraic equations of the collocation method:

We represent the values of the derivative of the desired function (polynomial) at the collocation points in the form of the product Dp ¼ p 0 of the matrix D by the vector of values of p. Recall that the Chebyshev differentiation matrix D has the standard representation in the physical space [3] (7) we reduce it to a system of linear algebraic equations

Here E is an identity matrix, f is a vector of values of the amplitude function on the grid. Denote by B the differentiation matrix in the frequency (spectral) space [13] , whose coefficients are explicitly expressed as

where 0 i; j n and r i ¼

& Denote by T the Chebyshev matrix of mapping a point (vector) from the space of coefficients to the space of values of the function [14] . Given that p ¼ Tc is the vector of values of the desired function (also in physical space), the components of the derivative vector can be written as Dp ¼ TBc [14] . As a result, we obtain the system of linear algebraic equations equivalent to system (9),

which is valid for an arbitrary grid on the interval [−1, 1]. We write Eq. (11) in detail 

where to reduce the formulas we used the notation T kj ¼ T k x j À Á ; k; j ¼ 0; . . .; n.

The product of a non-degenerate matrix T by a non-degenerate triangular matrix B þ ixE is a non-degenerate matrix. Therefore, the system of linear algebraic Eqs. (12) has a unique solution. Statement 1. The solution of this system of linear algebraic equations with respect to the coefficients c ¼ c 0 ; c 1 ; . . .; c n ð Þallows us to approximate the antiderivative function in the form of a series (5) and calculate the approximate value of the integral by formula (4).

System (12) is valid for an arbitrary grid on the interval [−1, 1]. However, consideration of the collocation problem on a Gauss-Lobatto grid allows significant simplification of this system of linear algebraic equations. First, we multiply the first and last equations from (12) by 1= ffiffi ffi 2 p to obtain an equivalent "modified" system with a new matrixT (instead of T), which is good because it has the property of discrete "orthogonality" and, therefore, is non-degenerate. Therefore, multiplying it on the left by its transposed one gives the diagonal matrix: We use this property and multiply the reduced (modified) system (12) on the left by the transposed matrixT T , thereby reducing it to the upper triangular form. Indeed, the matrix of the resulting system is calculated as the product of the diagonal matrix by the triangular matrix, which, in turn, is the sum of the Chebyshev differentiation matrix in the spectral space and the diagonal matrix.

Since the matrixT T is non-degenerate, the new system of linear algebraic equations is equivalent to system (12) and has a unique solution.

Taking into account the specific values of the Chebyshev polynomials on the Gauss-Lobatto grid [15] , simplifies the system, bringing it to the form ð13Þ

. . .; n and symbol P 00 denotes a sum in which the first and last terms are additionally multiplied by 1/2. By the Kronecker-Capelli theorem, the system of linear algebraic Eqs. (13) with a square matrix and a non-zero determinant is not only solvable for any vector of the right-hand side, but also has a unique solution. (13) has a stable solution.

Statement 3. To solve system (13) , no more than ( $ n 2 =4) operations of addition/subtraction and multiplication/division with a floating point are required.

Algorithms for solving systems of linear equations such as the Gauss method or the LU-decomposition work well when the matrix of the system has the property of diagonal dominance. Otherwise, standard solution methods lead to the accumulation of rounding errors. A stable solution to the system may be provided by the LU-decomposition method with a partial choice of a leading element. For the triangular matrix (13) , there is no forward pass (in which the leading element is selected) of the LU-decomposition, and the back pass of the method is always implemented without the selection of leading element. A solution to system (13) can still be unstable in the case when ix j j 2n. Passing to the solution of the normal system, that is, to the problem of minimizing the residual Ac Àf 2 , multiplying the system of Eqs. (13) on the left by the Hermitian conjugate matrix

we transform the matrix of the system (13) to the Hermitian form.

Although the system of linear equations became more filled, since instead of upper triangle three-diagonal matrix a system of linear equations with all matrix elements filled appeared, its computational properties are cardinally improved. The resulting matrix of a system of linear algebraic equations is Hermitian, its eigenvalues are real, and the eigenvectors form an orthonormal system. The method of LU-decomposition with a partial choice of the leading element, due to the properties of the resulting matrix, provides [17] the stability of the numerical algorithm for finding the only solution to the system.

Let us describe the sequence of operations of the presented algorithm for calculating the integral of a rapidly oscillating function of the form (1) with a linear phase.

Input data preprocessing.

1. If the integral is given on the interval a; b ½ , we pass to the standard domain of integration À1; 1 ½ by changing the variables x ¼ bÀa

Fill by columns the Chebyshev transformation matrix T from (12) using only one pass of the recursive method for calculating the values of Chebyshev polynomials of the first kind of the n-th order.

3. Calculate the vector of the right-hand side of system (13). 4. Fill in the elements of the sparse matrix (13) , which depend only on the dimension n and the phase value x. 5. If x j j[ 2n, then go to step 6. Otherwise go to step 7. 6. The matrix of system (13) is a matrix with a diagonal dominance and can be stably solved. The solution values at the boundary points are used to determine the desired antiderivative values. Go to step 8. 7. Multiply relation (13) on the left by the conjugate matrix to obtain a Hermitian matrix (14) with diagonal dominance. In this case, to determine the values of the antiderivative at the boundary points the normal solution is stably determined using the LU-decomposition with a partial choice of the leading element. 8. We calculate the values of the antiderivative at the ends of the interval using the formulas p 1 ð Þ ¼ P n j¼0 c j ; and p À1 ð Þ ¼ P n j¼0;jÀeven c j À P n j¼0;jÀodd c j . The desired value of the integral is obtained using the formula I f ; x ð Þ ¼ p 1 ð Þe ix À p À1 ð Þe Àix .

6 Numerical Examples Example 1. We give an example of calculating the integral when, for a good polynomial approximation of a slowly oscillating factor of the integrand, it is necessary to use polynomials of high degrees.

This integral is given by Olver ( [16] , p. 6) as an example of the fact that the GMRES method allows one to calculate the integral much more accurately than the Levin collocation method. However, in his article, solving the resulting system of linear algebraic equations requires O n 3 ð Þ operations, as in the Levin collocation method using the Gaussian elimination algorithm (Table 1) .

A comparison of our results at 40 interpolation points with the results of [16] shows a significant gain in accuracy: the deviation from the exact solution is of the order of 10 À17 compared with the deviation of the order of 10 −7 in Olver's article. The proposed algorithm to achieve an accuracy of 10 −13 in the calculation of the integral uses no more than 30 points (n 30Þ for x = 1,…, 100.

Example 2. As a second example, we consider the integral

from [16] , where the results of calculating the integrals depending on the number of approximation points are illustrated on Fig. 1 [16] . À Á and the integral can be written as:

Let us consider the calculation of this integral for various values of the parameter x using an algorithm that takes into account the linearity of the phase function ( Table 2 ).

The proposed algorithm to achieve an accuracy of 10 À16 when calculating the integral uses no more than 90 points (n 90) with x = 0.1, …, 100. A significant gain in the number of addition/subtraction and multiplication/division operations is achieved when the frequency value is greater than the number n, which ensures the diagonal dominance of the system of linear algebraic equations in the matrix (13) (Fig. 2) .

It is useful to compare the algorithm we developed for finding the integrals of rapidly oscillating functions with the results of [5] , which presents various and carefully selected numerical examples for various classes of amplitude functions.

Example 3. Consider the calculation of the integral with an exponential function as the amplitude

ð Þ e ixx dx; a ¼ 16; 64; x ¼ 20; 1000:

The exact value of the integral can be calculated by the formula I a; x ð Þ ¼ 2Ãe Àa sinh a þ ix ð Þ a þ ix ð Þ [5] . The plot of the deviation of the integral calculated by us from the exact one depending on the number of collocation points (absolute error) is shown in Fig. 3 . Comparison with the results of [5] shows that the accuracy of calculating the integrals practically coincides with that of [5] . Our advantage is the much simpler form of the matrix of a system of linear equations. In the best case, when |x| > n matrix of the system is a triangular matrix with a dominant main diagonal. If n > |x|, then the transition to the search for a normal solution to a system with a positively defined Hermitian matrix allows us to create a numerically stable solution scheme. Example 4. In this example [5] , the rapidly oscillating function e i2pax is considered as the amplitude one. It is clear that in this case, to achieve the same accuracy in calculating the integral as in the previous example, a larger number of collocation points will be required (Fig. 4) .

e i2pax e ixx dx; a ¼ 5; 10; x ¼ 20; 1000: ð19Þ Example 5. The example demonstrates the calculation of the integral in the case when the amplitude function is the generating function of the Chebyshev polynomials of the first kind.

1 À a 2 1 À 2ax þ a 2 e ixx dx; a ¼ 0:8; 0:9; x ¼ 20; 1000:

The behaviour of the amplitude function should lead to an almost linear dependence of the approximation accuracy on the number of points for various values of the parameters a and x. Numerical experiments carried out confirm this assertion (Fig. 5 ). Moreover, the accuracy of calculating the integrals is not inferior to the accuracy of the methods of [5] .

The example is rather complicated for interpolation by Chebyshev polynomials. To achieve acceptable accuracy (10 À18 ), the deviation of the calculated value of the integral from the exact one requires about 300 approximation points both for small values of x ¼ 20 and for large x ¼ 1000 (Fig. 6 ). Example 7. We give an example of integration when the amplitude function has second-order singularities at both ends of the integration interval

The value of this integral can be calculated in an analytical form: f x ð Þ ¼ 3pJ 2 x ð Þ=x 2 . We present the numerical values of the integral for various values of the frequency: I 20 ð Þ ¼ À0:00377795409950960, I 1000 ð Þ¼À2:33519886790130  10 À7 (Fig. 7) . Similar to the previous example, to achieve good accuracy in calculating the integral, it is necessary to consider a large number of collocation points. However, for this type of amplitude functions, the method presented in the article works reliably both in the case of low and high frequencies.

The given examples demonstrate that the dependence of the solution on the number of approximation points is similar to the dependence demonstrated in Olver [17] and Hasegawa [5] . The advantage of our approach is the simplicity of the algorithm and the high speed of solving the resulting very simple system of linear algebraic equations. If it is necessary to repeatedly integrate various amplitude functions at a constant frequency, multiple gains are possible due to the use of the same LU-decomposition backtracking procedure. 

A simple, effective, and stable method for calculating the integrals of highly oscillating functions with a linear phase is proposed. It is based on Levin's brilliant idea, which allows the use of the collocation method to approximate the antiderivative of the desired integral. Using the expansion in slowly oscillating polynomials provides a slowly changing solution of the differential equation.

The transition from a solution in physical space to a solution in spectral space makes it possible to effectively use the discrete orthogonality property of the Chebyshev mapping matrix on a Gauss-Lobatto grid. With this transformation, the uniqueness of the solution of the studied system is preserved, and its structure from a computational point of view becomes easier.

There are a large number of works using various approaches aimed to offer fast and effective methods for solving SLAEs that arise when implementing the collocation method. However, many methods [5, 6] encounter instability when solving the corresponding systems of linear equations. When using Chebyshev differentiation matrices in physical space, instability is explained primarily by the degeneracy of these matrices and the huge spread of eigenvalues of the matrix of the collocation method system. The approach to solving the differential equation based on the representation of the solution, as well as the phase and amplitude functions, in the form of expansion in finite series by Chebyshev polynomials and the use of three-term recurrence relations [5, 6, 8] also does not provide a stable calculation for n > |x|. To overcome instability, various methods of regularizing the systems under study are proposed.

In our work, we propose a new method for improving computational properties by preconditioning of the system in the spectral representation and by searching for its pseudo-normal solution. The proposed method has been reduced to solving a SLAE with a Hermitian matrix. A number of numerical examples demonstrates the advantages of the proposed effective stable numerical method for integrating rapidly oscillating functions with a linear phase.

Procedures for computing one-and two-dimensional integrals of functions with rapid irregular oscillations

Filon and Levin methods In: Computing Highly Oscillatory Integrals

Chebyshev Polynomials

Computing Highly Oscillatory Integrals

A user-friendly method for computing indefinite integrals of oscillatory functions

Stability and error estimates for Filon-Clenshaw-Curtis rules for highly oscillatory integrals

An improved Levin quadrature method for highly oscillatory integrals

A well-conditioned Levin method for calculation of highly oscillatory integrals and its application

Fast integration of rapidly oscillatory functions

Moment-free numerical integration of highly oscillatory functions

A comparison of some methods for the evaluation of highly oscillatory integrals

An alternative method for irregular oscillatory integrals over a finite range

A Practical Guide to Pseudospectral Methods

Regularized computation of oscillatory integrals with stationary points

Integration of highly oscillatory functions

Fast, numerically stable computation of oscillatory integrals with stationary points

GMRES shifted for oscillatory integrals

Acknowledgement. The publication has been prepared with the support of the "RUDN University Program 5-100".