CC BY
5
Discrete & Continuous Models
#& Applied Computational Science 2021, 29 (3) 251-259
ISSN 2658-7149 (online), 2658-4670 (print) http://journals-rudn-ru/miph
Research article
UDC 519.872:519.217
PACS 07.05.Tp, 02.60.Pn, 02.70.Bf
DOI: 10.22363/2658-4670-2021-29-3-251-259
Asymptotically accurate error estimates of exponential convergence for the trapezoidal rule
Aleksandr A. Belov1,2, Valentin S. Khokhlachev1
1 M. V. Lomonosov Moscow State University 1, bld. 2, Leninskie Gory, Moscow, 119991, Russian Federation 2 Peoples' Friendship University of Russia (RUDN University) 6, Miklukho-Maklaya St., Moscow, 117198, Russian Federation
(received: August 19, 2021; accepted: September 9, 2021)
In many applied problems, efficient calculation of quadratures with high accuracy is required. The examples are: calculation of special functions of mathematical physics, calculation of Fourier coefficients of a given function, Fourier and Laplace transformations, numerical solution of integral equations, solution of boundary value problems for partial differential equations in integral form, etc. For grid calculation of quadratures, the trapezoidal, the mean and the Simpson methods are usually used. Commonly, the error of these methods depends quadratically on the grid step, and a large number of steps are required to obtain good accuracy. However, there are some cases when the error of the trapezoidal method depends on the step value not quadratically, but exponentially. Such cases are integral of a periodic function over the full period and the integral over the entire real axis of a function that decreases rapidly enough at infinity. If the integrand has poles of the first order on the complex plane, then the Trefethen-Weidemann majorant accuracy estimates are valid for such quadratures.
In the present paper, new error estimates of exponentially converging quadratures from periodic functions over the full period are constructed. The integrand function can have an arbitrary number of poles of an integer order on the complex plane. If the grid is sufficiently detailed, i.e., it resolves the profile of the integrand function, then the proposed estimates are not majorant, but asymptotically sharp. Extrapolating, i.e., excluding this error from the numerical quadrature, it is possible to calculate the integrals of these classes with the accuracy of rounding errors already on extremely coarse grids containing only ~ 10 steps.
Key words and phrases: trapezoidal rule, exponential convergence, error estimate, asymptotically sharp estimates
1. Introduction
Applied tasks. In many physical problems it is needed to calculate integrals, that cannot be obtained in terms of elementary functions. Here are some examples:
© Belov A. A., Khokhlachev V.S., 2021
This work is licensed under a Creative Commons Attribution 4.0 International License http://creativecommons.org/licenses/by/4-0/
1) Calculation of special functions of mathematical physics: the Fermi-Dirac functions, which are equal to the moments of the Fermi distribution, the Gamma function, cylindrical functions and a number of others.
2) Calculation of the Fourier coefficients of a given function, Fourier and Laplace transform.
3) Numerical solution of integral equations, both well-posed and ill-posed.
4) Solving boundary value problems for partial differential equations (including eigenvalue problems) written in integral form, etc.
Calculation of quadratures. Trapezoidal rule, rectangle rule and Simpson's rule are commonly used for grid computation of quadratures. Usually the error of these methods quadratically depends on the grid step, and a large number of steps is needed to obtain good accuracy.
However, there are a number of cases when the error of the trapezoidal rule depends on the grid step exponentially, i.e. when the step is reduced by half, the number of correct signs of the numerical result is approximately doubled. This rate of convergence is similar to that of Newton's method. Two such cases are known. These are: the integral of the periodic function over the full period and the improper integral of a function that decreases rapidly enough at infinity.
If the integrand has first order poles on the complex plane, then for such quadratures there are majorant error estimates of Trefethen and Weidemann [1], see also [2]-[10]. In [11], [12] the generalization of Trefethen and Weidemann estimates is built for the case when the nearest pole of an integrand function is multiple.
In this paper, new error estimates of exponentially convergent quadratures of periodic functions over the full period are described. Integrand function can have an arbitrary number of poles of an integer order on the complex plane. If the mesh is detailed enough and the profile of the integrand resolved well, then the proposed estimates are not majorant, but asymptotically accurate.
It is possible to calculate the integrals of the indicated classes with the accuracy of round-off errors even on extremely coarse grids containing only ~ 10 steps by extrapolation (i.e., subtraction) of this error from the numerical value of the quadrature.
2. Exponentially convergent quadratures
One of the classes of exponentially convergent quadratures are integrals of periodic functions over the full period. By replacing z = exp (2nix/X) we move from the integral over the period [0, X] to the integral over the unit circle |z| = 1 on the complex plane. We choose the bypass direction of this circle counter clockwise. In [1], the following statement is formulated and proved:
Theorem 1. Let u(z) be analytic in the ring Rr1 < |z| < R, where R > 1, and lu(z)l < M0. We introduce a uniform grid on the unit circle ^^ — exp (2nin/N) , n = 0, N. Consider the integral and the trapezoidal rule quadrature
1= I u(z)^~, TN = Jj^u(zn).
i n=1 \z\ = 1
Then the estimate for the quadrature error holds
4nM0 RN -1
A \T T | < 4KMo (
It is obvious that, by replacing z = exp (ix), theorem 1 holds for integral over full period of the function u (exp (ix)) on the real axis.
In the works [11], [12], it was shown that the dependence of the estimate (1) from N can be not majorant, but asymptotically accurate. This holds, if u (z) has only first order pole type singularities, and R is taken such that the closest singularity to the unit circle lies on the boundary of the ring R-i < \z\ < R. In this case, the integrand function increases significantly if one approaches the singularity from inside the ring. Thereby, the constant M0 loses its usual meaning from theorem 1. We carefully studied proof of the theorem 1 given in [1] and we found the possibility of significant strengthening the results of this theorem, under some additional conditions on the integrand function.
3. Calculating the error
Let us consider in detail the contour integral over the unit circle of a function that has one simple pole inside it and another simple one pole outside it. This case corresponds to the integral considered in [1]. Suppose the point a1 is inside, and the point a2 is outside \z\ = 1 and, the function u (z) is analytic in the ring R-1 < \z\ < R, where R = min{1/ \a1 \, \a2\}. Then the integral has the form
G= 4 g(z)dz= i 7-^-,dz = 2^ U(ai1
(z-ai)(z-a,2) (a-i -a2)'
\z\ = 1 \z\ = 1
We make one assumption for the sake of simplifying the calculations. Its effect on the result is weak. Let u(z) = 1, then we rewrite the integrand function in this form
g ^ = 1 =_1_+_1_
(z-ai )(z-a,2) (a-i - a,2) (z - ai) (a,2 - aj (z - av)'
Now we decompose each fraction in the Laurent series as the sum of the geometric progression
1 ^ at1 1 ^ zk* 9(z) = ^--r > . -r—i -~T--r > ,
K -o>2) f=0 zkl+i (a>2 -ai) f=0 ak22+i'
We use the grid zn, n = 0,N, which is introduced in theorem 1. Our goal is to obtain an explicit expression for the grid step Azn = zn+1 — zn
Azn = exp
2ni (n + 1) N
exp
f2mn\ _ (2m i 1 [~N~) N=o Zn [~N +U\N2
). (2)
Discarding the 0 (N 2) term in the expression for the grid step, we write the trapezoidal rule quadrature in the following form
N-l
Gn = ^g(zn )Az%
n=0
2ni
• N-l
E9(zn).
n=0
We substitute the representation of g(zn) by the sum of the series in the quadrature and then swap the series and the finite sum. Last step is allowed due to absolute convergence of the resulting double number series (each member of the double series of modules can be estimated by the corresponding member of an infinitely decreasing geometric progression, which has finite sum). The following expression for the quadrature formula is obtained
Gn =
2ni
N(ai — 0,2)
oo N-l Sj — — a 1
o N-l s2 + 1 Zn
EE £ + EE
S2 + 1
s, =0 n=0
s2=0 n=0 a2
(3)
To perform these transformations, we need the following well known result
N-l
E exp (±2^i
n=0
nk ~N
N, k is a multiple of N, otherwise.
We convert the second sum in square brackets in the formula (3)
o N-l s2 + l Zn
2Lc 2Lc s2+l
s2=0 n=0 a2
(s2 + 1) is a multiple of N,'
= ^ (*2 + 1)=NP2,
^2 = 1,<X)
"E
- HNP2 P2=l U2
1 = N- 1/a
N 2
1 — 1/a2 .
We convert the first sum in (3) ^ alJ (sl is a multiple of N,
s,=0 n=0
We get
Si = Npi, Pi =
N E ^ = N-
Pj=0
1 — a?.
Gn =
2ni
(a1 — 0-2 )
1 — a
11
+
N
a2 —1
n
Finally, we calculate the quadrature error
AN = G-GAj = 2m
N
1
(a1 - 0-2 )
1 -
1
N
1 - aI
+
1
(a2
a2 -1
The obtained result can be easily generalized to the case when the function u (z) does not equal identically to one. The derivation is similar but far too cumbersome. Let us formulate the final result.
Theorem 2. Let the point z = a1 be inside the unit circle, and let the point z = a2 be outside of it. Let the function u (z) be analytic on the entire complex plane, with the possible exception of an infinitely distant point, and u ^ 0 at the points z = a1 2. Then the trapezoidal rule for the integral G has the following error estimation
An = G - Gn = = 2ni
u (a1 )
K - a2)
1-
1-<j
+
u (a2) 2 - a,i )
La2 -1
(4)
Estimate (4) is not majorant, but asymptotically accurate. The only one approximation that was made is contained in the approximate expression for the grid step (2).
4. Validation
Calculations were carried out with the test integral having a known value
J =
sin (z)
(z-ai ) (z-a.2 )
dz = 2ni
. sin (a1 )
(5)
where a1 = 0.6 + 0.6i and a2 = 2 —i. In this case, 1/ la1l & 1.2 and la21 & 2.2, so the value R from theorem 1 equals 1/ la11. During the calculations, the following information was obtained: actual error, the Trefethen-Weidemann estimate (1), our estimate (4) and the error after extrapolation.
The figure 1 shows quadrature error versus number of grid steps in the semi-logarithmic scale. Here, the black dots represent the actual error, the white circles represent our estimate, and the black squares represent the Trefethen-Weidemann estimate with the constant M0 = 1. Recall that this constant loses its meaning from theorem 1, if the singularity lies on the boundary of the ring.
The plot shows that our estimate coincides with the actual error already at N > 4. The Trefethen-Weidemann estimate does not represent the initial part of the curve. It describes the curve starting from N = 15. This estimate is asymptotically accurate in N, but the true value of the constant M0 is unknown. Therefore, the Trefethen-Weidemann estimate cannot be used for extrapolation. Thereby, the error estimate constructed in this paper is much stronger than the Trefethen-Weidemann estimate.
These conclusions are also confirmed by the figure 2. Here, we plot the ratio of error estimates to actual accuracy versus number of grid steps. The number 1 corresponds to the Trefethen-Weidemann estimate and the number 2 is for our estimate. It can be seen that when N > 4 our estimate is almost indistinguishable from the actual error. Therefore, it can be excluded from the quadrature (i.e. extrapolated). This dramatically increases the accuracy of the calculation. One can also see that the Trefethen-Weidemann estimate significantly less accurate in assessing the dependence of the error on the number of nodes: the corresponding relation goes out to a constant on the much more detailed grids than the estimate (4).
Figure 1. Graph of convergence of the Figure 2. Ratio of theoretical estimates (1) trapezoidal formula. Symbols are described and (4) for the integral (5) to actual in the text accuracy versus number of grid steps.
Symbols see in the text
5. Extrapolation of the error
Let us exclude the error (4) from the calculated quadrature by the formula
Gn = Gn + An . (6)
This is equivalent to introducing some new quadrature formula. The error is shown in the figure 1 by white triangles. One can see that the speed of convergence of the quadrature (6) radically exceeds even the exponential one. The accuracy of round-off errors is achieved already at N ~ 15, which is ~ 10 times less than for the trapezoidal rule. Labor intensity of such computation is comparable to the complexity of explicit formulas. This approach is essentially new and exceeds the world level.
6. Conclusion
The described method is a powerful tool for solving physical problems. If one can find transformation of variables that reduce integral to one of the considered types, then the calculations are accelerated thousands of times. In this paper 1) a fundamentally new estimate of the error of quadrature
is constructed, it is asymptotically accurate. 2) Extrapolation procedure is proposed, which provides calculation of the quadrature with the accuracy of unit errors rounding, and it is already performed on very rough grids with the number of steps from 5-15.
Acknowledgments
This work was supported by grant MK-3630.2021.1.1.
References
[1] L. N. Trefethen and J. A. C. Weideman, "The exponentially convergent trapezoidal rule," SIAM Review, vol. 56, no. 3, pp. 385-458, 2014. DOI: 10.1137/130932132.
[2] J. Mohsin and L. N. Trefethen, "A trapezoidal rule error bound unifying the Euler-Maclaurin formula and geometric convergence for periodic functions," in Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences, vol. 470, 2014, p. 20130 571. DOI: 10.1098/ rspa.2013.0571.
[3] J. A. C. Weideman, "Numerical integration of periodic functions: A few examples," The American Mathematical Monthly, vol. 109, no. 1, pp. 2136, 2002. DOI: 10.2307/2695765.
[4] N. Eggert and J. Lund, "The trapezoidal rule for analytic functions of rapid decrease," Journal of Computational and Applied Mathematics, vol. 27, no. 3, pp. 389-406, 1989. DOI: 10.1016/0377-0427(89)90024-1.
[5] H. Al Kafri, D. J. Jeffrey, and R. M. Corless, "Rapidly convergent integrals and function evaluation," Lecture Notes in Computer Science, vol. 10693, pp. 270-274, 2017. DOI: 10.1007/978-3-319-72453-9_20.
[6] J. Waldvogel, "Towards a general error theory of the trapezoidal rule," in Springer Optimization and Its Applications. 2010, vol. 42, pp. 267282. DOI: 10.1007/978-1-4419-6594-3_17.
[7] E. T. Goodwin, "The evaluation of integrals of the form f °° f(x)e-x dx," Mathematical Proceedings of the Cambridge Philosophical Society, vol. 45, no. 2, pp. 241-245, 1949. DOI: 10. 1017/ S0305004100024786.
[8] N. N. Kalitkin and S. A. Kolganov, "Quadrature formulas with exponential convergence and calculation of the Fermi-Dirac integrals," Doklady Mathematics, vol. 95, no. 2, pp. 157-160, 2017. DOI: 10.1134/ S1064562417020156.
[9] N. N. Kalitkin and S. A. Kolganov, "Refinements of precision approximations of Fermi-Dirak functions of integer indices," Mathematical Models and Computer Simulations, vol. 9, no. 5, pp. 554-560, 2017. DOI: 10.1134/S2070048217050052.
[10] N. N. Kalitkin and S. A. Kolganov, "Computing the Fermi-Dirac functions by exponentially convergent quadratures," Mathematical Models and Computer Simulations, vol. 10, no. 4, pp. 472-482, 2018. DOI: 10.1134/S2070048218040063.
[11] A. A. Belov, N. N. Kalitkin, and V. S. Khokhlachev, "Improved error estimates for an exponentially convergent quadratures [Uluchshen-nyye otsenki pogreshnosti dlya eksponentsial'no skhodyashchikhsya kvadratur]," Preprints of IPM im. M. V. Keldysh, no. 75, 2020, in Russian. DOI: 10.20948/prepr-2020-75.
[12] V. S. Khokhlachev, A. A. Belov, and N. N. Kalitkin, "Improvement of error estimates for exponentially convergent quadratures [Uluchsheniye otsenok pogreshnosti dlya eksponentsial'no skhodyashchikhsya kvadratur]," Izv. RAN. Ser. fiz, vol. 85, no. 2, pp. 282-288, 2021, in Russian. DOI: 10.31857/S0367676521010166.
For citation:
A. A. Belov, V. S. Khokhlachev, Asymptotically accurate error estimates of exponential convergence for the trapezoidal rule, Discrete and Continuous Models and Applied Computational Science 29 (3) (2021) 251-259. DOI: 10.22363/2658-4670-2021-29-3-251-259.
Information about the authors:
Belov, Aleksandr A. — Candidate of Physical and Mathematical Sciences, Assistant professor of Department of Applied Probability and Informatics of Peoples' Friendship University of Russia (RUDN University); Researcher of Faculty of Physics, M. V. Lomonosov Moscow State University (e-mail: aa.belov@physics.msu.ru, phone: +7(495)9393310, ORCID: https://orcid.org/0000-0002-0918-9263, ResearcherID: Q-5064-2016, Scopus Author ID: 57191950560) Khoklachev, Valentin S. — Master's degree student of Faculty of Physics, M. V. Lomonosov Moscow State University (e-mail: valentin.mycroft@yandex.ru, phone: +7(495)9393310, ORCID: https://orcid.org/0000-0002-6590-5914)
УДК 519.872:519.217
PACS 07.05.Tp, 02.60.Pn, 02.70.Bf
DOI: 10.22363/2658-4670-2021-29-3-251-259
Асимптотически точные оценки экспоненциальной сходимости для формулы трапеций
А. А. Белов1 2, В. С. Хохлачев1
1 Московский государственный университет им. М. В. Ломоносова Ленинские горы, д. 1, стр. 2, Москва, 119991, Россия 2 Российский университет дружбы народов ул. Миклухо-Маклая, д. 6, Москва, 117198, Россия
Во многих прикладных задачах требуется экономичное вычисление квадратур с высокой точностью. Примерами являются: вычисление специальных функций математической физики, расчёт коэффициентов Фурье заданной функции, преобразования Фурье и Лапласа, численное решение интегральных уравнений, решение краевых задач для уравнений в частных производных в интегральной форме и т.д. Для сеточного вычисления квадратур обычно используют методы трапеций, средних и Симпсона. Обычно погрешность этих методов зависит от шага степенным образом, и для получения хорошей точности требуется большое число шагов. Однако существует ряд случаев, когда погрешность метода трапеций зависит от величины шага не квадратично, а экспоненциально. Такими случаями являются интеграл от периодической функции по полному периоду и интеграл по всей числовой прямой от функции, достаточно быстро убывающей на бесконечности. Если подынтегральная функция имеет полюса первого порядка в комплексной плоскости, то для таких квадратур справедливы мажорантные оценки точности Трефетена и Вайдемана.
В работе построены новые оценки погрешности экспоненциально сходящихся квадратур от периодических функций по полному периоду. Подынтегральная функция может иметь произвольное число полюсов целого порядка на комплексной плоскости. Если сетка достаточно подробная (разрешает профиль подынтегральной функции), то предлагаемые оценки являются не мажорантными, а асимптотически точными. Экстраполируя, то есть исключая эту погрешность из численной квадратуры, можно вычислять интегралы указанных классов с точностью ошибок округления уже на чрезвычайно грубых сетках, содержащих всего ~ 10 шагов.
Ключевые слова: формула трапеций, экспоненциальная сходимость, оценки точности, асимптотически точные оценки
