CC BY
3
Discrete & Continuous Models & Applied Computational Science 2023, 31 (2) 128-138
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-2023-31-2-128-138 EDN: XAUSJA
Quadratures with super power convergence
Aleksandr A. Belov1'2, Maxim A. Tintul1, Valentin S. Khokhlachev1
1 M. V. Lomonosov Moscow State University, 1, bld. 2, Leninskie Gory, Moscow, 119991, Russian Federation
2 RUDN University, 6, Miklukho-Maklaya St., Moscow, 117198, Russian Federation
(received: April 25, 2023; revised: May 5, 2023; accepted: June 26, 2023)
Abstract. The calculation of quadratures arises in many physical and technical applications. The replacement of integration variables is proposed, which dramatically increases the accuracy of the formula of averages. For infinitely smooth integrand functions, the convergence law becomes super power. It is significantly faster than the power law and is close to exponential one. For integrals with bounded smoothness, power convergence is realized with the maximum achievable order of accuracy.
Key words and phrases: trapezoid rule, exponential convergence, error estimate, asymptotically sharp estimates
1. Introduction
Applied tasks. In many physical problems, it is required to approximate integrals that are not taken in elementary functions. Here are some examples:
1. Calculation of special functions of mathematical physics: Fermi-Dirac functions equal to the moments of the Fermi distribution, gamma function, cylindrical functions and a number of others.
2. Calculation of Fourier coefficients of a given function, Fourier and Laplace transforms.
3. Numerical solution of integral equations, both correctly posed and incorrect.
4. Solving boundary value problems for partial differential equations (including eigenvalue problems) written in integral form, etc.
Such integrals must be calculated with high accuracy up to computer round-off errors.
© Belov A. A., Tintul M.A., Khokhlachev V.S., 2023
This work is licensed under a Creative Commons Attribution 4.0 International License https://creativecommons.org/licenses/by-nc/4-0/legalcode
Calculation of quadratures. Commonly, trapezoid, mean and Simpson methods on a uniform grid are used for grid calculation of quadratures. The majorant error estimation is well known for these methods. For trapezoid and mean formulas it is 0(h2), for Simpson's formula it is 0(h4), where h is the grid step. There are ways to improve accuracy: calculation on a set of thickening grids and extrapolation refinement by the Richardson method, refinement by the Euler-Maclaurin formula, etc. [1, 2]. All these methods give a power dependence of the error on the grid step 0(hm).
If the integrand is periodic and the integral is calculated over the full period, then the dependence of the error on the step becomes exponential instead of power-law ~ exp(-1/h) [3-5]. This means that when the step is halved, the number of correct characters in the answer approximately doubles. This convergence rate is much faster than the power one. However, the corresponding class of integrand functions is rather narrow. Attempts have been made in the literature to expand this class [6-9], but they were considered unsuccessful [7].
In the present paper, an approach is proposed that dramatically accelerates the convergence of the mean rule. It is based on a special substitution of integration variables. The integrand function may be non-periodic. If it is infinitely smooth, then the proposed replacement provides super power convergence of the quadrature. This convergence rate is significantly faster than the power-law one and is close to the exponential one.
If the integrand has bounded smoothness, then the proposed method gives a power convergence with the maximum achievable order of accuracy.
The proposed approach does not require a priori information about the nature of the integrand function and is uniformly applicable to a wide range of tasks. The class of integrand functions, for which the super power convergence of quadratures is realized, is significantly expanded.
2. Change of integration variables
Consider the integral
1= i f(x)dx. (1)
Jo
Let us perform the variable change in two stages. First, using fractional polynomial transformation t(x), we map the segment x £ (0,1) to the straight line t £ +<x)). Then we map this line to the segment £ £ (0,1) using the transformation £(£), whose derivatives tend to zero near £ = 0 and £ = 1 faster than any degree .
Such substitutions can be made in various ways. In this paper, the following transformation was considered
№ = Pl--^ x(f) = 1 +1 th(^ (2)
where A, B, a are constants. The mapping (2) is shown in figure 1 as z(£) dependence. It is almost linear in the middle of the segment, but at its ends,
the derivatives of x^ quickly tend to zero. It is also possible to implement the replacement (2), in which the error function &(Bt) is taken instead of the hyperbolic tangent.
S
Figure 1. Variable transformation (2). Parameters A, B, a are equal to unity
After mapping (2), the integral takes the form
1= T f(№, № = f{x[W]}xtimk(0. (3)
Periodic continuation. Let us show that the new integrand f(£) admits an infinitely smooth periodic continuation beyond the boundaries of the segment
£e(0,1).
The expression t^(^) ~ (1 — £)-a-1 has poles at the ends of the
segment £ = 0 and £ = 1. However, for £ ^ 0 + 0 and £ ^ 1 — 0, the derivative xt ~ exp(—£-a (1 — ) tends to zero significantly faster. As a result, xtt^ ^ 0 when striving for points £ = 0 and £ = 1 from inside the
segment. Therefore, f(^) vanishes at the boundaries of the segment.
Similarly, it can be shown that all derivatives of this function tend to zero at £ ^ 0 + 0 and £ ^ 1 — 0. For example, the first derivative has the form
/_£ = fxX212 + fxtit\ + fxt t^. (4)
All derivatives of dtm/d^m ~ at the boundaries of the seg-
ment have poles that are multiplied by the expression ~ exp(—£-a(1 — ) in various degrees. Therefore, for £ ^ 0 + 0 and £ ^ 1 — 0 we have ft ^ 0.
The same is true for higher derivatives df™/dThus, the integrand function f can be periodically continued infinitely smoothly beyond the boundaries of the segment £ £ (0,1).
3. Mean rule convergence
On the segment £ £ (0,1), we introduce a uniform grid with a step h = 1/N. Half-integer nodes are denoted by £n+1/2 = (n — 1/2)h, n = 1, ...,N. We write the mean rule quadrature
N
In = Y.hf(tn+i/2). (5)
n=1
The following statement holds.
Theorem 1.
A) If f(x) is infinitely smooth on the segment x £ (0,1), then the quadrature (5) has super power convergence.
B) If f(x) has j continuous derivatives on x £ (0,1), the (j + 1)-th derivative has a discontinuity at the point x = a £ (0,1), and this point is a grid node, then the quadrature (5) has power convergence. The order of accuracy is j + 2 if j is even, and j + 3 if j is odd. This order of accuracy is maximal for a given smoothness of the integrand function.
Proof. Let us prove the statement A). The power part of the mean rule error is described by the Euler-Maclaurin formula [1]. It contains the differences of odd derivatives at the ends of the integration segment
to
5 = Y bkh2k ( f(2k-1) (1) — f(2k-1) (0)) , bk = const. (6)
k= 1
As noted above, due to the variable transform (2), the derivatives f(k) (£) — 0 for £ —y 0 + 0 and £ — 1 — 0. All summands in the sum of (6) vanish. Therefore, there are no power terms left in the error of the mean rule, and the convergence turns out to be super power one.
Let us prove the statement B). Under these assumptions, the power-law contribution to the error of the mean formula has the form
to
S = Yhh2k (f(2k-1)(1) — f(2k-1)(0)) +
k=1
K
+ Y bkh2k (f(2k-1)(a — 0) — f(2k-1) ( a + 0)) . (7)
k=1
The first sum in (7) is similar to (6). After the variable transform (2), it turns to zero.
The second sum is the error resulting from the singularity at the point a. If 2 k — 1 < j, then by virtue of continuity, the right and left limit values of
derivatives of the order of 2k —1 are the same f(2k 1 (a — 0) = f(2k 1 (a + 0).
What is the limit of summation of K? Since f is discontinuous at the point a, and only odd derivatives are included in (7), two cases are possible. If j is odd, then 2K — 1 = j + 2. Then 5 = 0(h^+3). If j is even, then 2K — 1 = j + 1, and 5 = 0(№+2). Obviously, this order of accuracy is the maximum, i.e. it cannot be improved. The theorem is proved. n
Note. The literature describes [6-9] variable substitutions similar to (2). In these works, trapezoid and Simpson formulas were used, in which one needs to calculate the integrand function at the boundary points. However, after
variable change (2), the integral function f(^) has essentially singular points within the boundaries of the segment £ = 0 and £ = 1. Therefore, calculating
f(0) and f(1) presents a problem; in particular, computer numbers overflow occurs.
To avoid this, in [6] it was proposed to cut the integration segment, i.e. instead of £ G (0,1), consider £ G (e, 1 — e), where £ is some small number. Such a cutting introduced a significant error, and it was not possible to realize superstellar convergence. The authors of [7] conducted numerical experiments and found that this approach is inferior in quantitative accuracy to Simpson's formula without replacing variables. Therefore, this approach was considered unpromising [7].
We use the mean rule that does not require calculating f(0) and f(1). Therefore, the described difficulty does not arise, and super power convergence is realized.
Infinitely smooth integrand. As an example, consider a test integral with a known exact value
The integrand is infinitely smooth.
The calculation was carried out on a set of grids with different N = 2,4,8,... On each grid, the mean rule quadrature and its error A = II — IN|, equal to the difference between the numerical and exact integrals, were calculated. Figure 2 shows a graph of the error A depending on the number of grid steps N. The scale of the graph is semi-logarithmic. At this scale, exponential convergence corresponds to a straight line, and a power-law curve corresponds to a logarithmic curve.
Dark circles correspond to the calculation with the replacement of variables (2), light circles correspond to the calculation without it. One can see that the proposed replacement of variables dramatically increases accuracy: already at N ~ 100, the error is A ~ 10-14, which is comparable to rounding errors. The gain in accuracy compared to the calculation without replacing variables reaches 10 orders of magnitude. The convergence rate is somewhat inferior to the exponential one, but cardinally exceeds the power one.
4. Method validation
'0
-
-8
-10-
-12-
-14 -
-16
0
40
80
120
N
Figure 2. The error of the mean rule quadrature in the test (8). Notation are explained in the text
Due to the presence of essentially singular points of f at the boundaries of the segment, the dependence of the error on the number of steps is nonmonotonic and alternating [10, 11]. In this graph, this can be seen by the non-monotonic behavior of the curve. Local minima correspond to the change of the error sign.
Therefore, the proposed replacement dramatically increases the accuracy of the mean rule quadrature. We recommend it for wide application in practical computing.
Integrand function with bounded smoothness. Often in applications, it is necessary to calculate integrals from piecewise given spline approximations and interpolants. They have limited smoothness. So, the simplest linear interpolation is continuous, but has discontinuities of the first derivative. The cubic spline is continuous along with the second derivative, and the third derivative experiences a discontinuity.
As an example, consider the integral of the function
for integers 1 < m < 5. The function (9) has a m — 1 continuous derivative, and the m-th derivative experiences a discontinuity. The exact values of the integral I are known, they are listed in table 1.
The calculation was carried out on several thickening grids. They were chosen so that the feature x = 0.5 was a node. For example, it is enough to take only even N for this. The resulting errors depending on the number of steps are shown in figure 3. The scale of the graph is double logarithmic.
(9)
Table 1
Test (9)
m I Q
1 1 + 2e0'5 - e 2
2 1 - 8e0'5 + 5e 4
3 1 + 48e0'5 - 29e 4
4 1 - 384e0'5 + 233e 6
5 1 + 3840e0'5 - 2329e 6
Therefore, the power convergence corresponds to a straight line whose slope is equal to the order of accuracy. The numbers near the lines are the values of m.
—j—i—i—i—i—i—i—i—i—i—i—i—i—i—i—i—i—i—i—i—[—
0,0 0,5 1,0 1,5 2,0
igN
Figure 3. The error of the mean rule quadrature in the test 9. Notation are explained in the text
It can be seen that on sufficiently detailed grids, the curves for each m tend to straight lines, i.e. power convergence is realized. The corresponding orders of accuracy of q are given in the table 1. They are completely consistent with Theorem 1.
On coarse meshes, the behavior of curves is irregular. The error depends on N nonmonotonically, changes sign and decreases significantly faster than the power law. Apparently, the law of convergence on coarse meshes corresponds to the super power one (similar to the figure 2). Justification of this consideration is beyond the scope of the present work.
For comparison, the figure 3 shows the calculation error according to the mean rule without variable change. For all the m considered, the errors were approximately the same, so we showed them with one line. It is indicated by an asterisk (*). This line corresponds to a power convergence with a second order of accuracy. It can be seen that for m > 2 (i.e., if there is at least one continuous derivative), the proposed replacement increases the order of accuracy and sharply reduces the quantitative error. In Fig. 3, the accuracy gain was from 3 to 8 orders of magnitude.
We performed a similar calculation using grids in which the singularity x = 0.5 did not get into the node. To do this, it is enough to take only odd N. This case does not fall under Theorem 1. Nevertheless, the theorem turned out to be true for it as well. The resulting errors were similar to the figure 3. In particular, the convergence rate for the considered m turned out to be the same. The quantitative accuracy was somewhat worse than figure 3. This was most noticeable for m = 1. For other m, the accuracy decreasing turned out to be insignificant.
5. Conclusion
In this paper, a special transformation of the integration variable is proposed, which dramatically increases the accuracy of the mean rule quadrature. For infinitely smooth integrand functions, convergence becomes super power one. For functions of bounded smoothness, the convergence law remains power-law, but the maximum achievable order of accuracy is realized.
Let us conduct a qualitative comparison of the proposed approach with other methods for improving the accuracy of quadratures listed in the introduction.
None of them provides super power convergence. Therefore, for infinitely smooth functions, the proposed approach provides obviously higher accuracy.
The use of Euler-Maclaurin corrections requires a large amount of a priori information about the integrand. It is necessary to accurately calculate high derivatives and a priori set the number of corrections to be taken into account. Therefore, the maximum order of accuracy is realized if the smoothness class of the integrand function is known.
On the contrary, the proposed approach is uniformly applicable to integrals both infinitely smooth and having bounded smoothness. One does not need to know the smoothness class in advance.
It is possible to increase the order of accuracy using Richardson extrapolation only on sufficiently detailed grids on which theoretical convergence is already being implemented, but rounding errors have not yet been achieved. On coarse grids, the use of extrapolation can even degrade accuracy.
In the proposed method, even on coarse grids, convergence is observed, and quite fast. A quantitative comparison of the Richardson extrapolation and the proposed method for bounded smoothness functions is beyond the scope of this paper.
Acknowledgments
This work did not receive specific funding.
References
[1] N. N. Kalitkin and E. A. Alshina, Numerical Methods. Vol. 1: Numerical Analysis [Chislennye Metody. T. 1: Chislennyi analiz]. Moscow: Akademiya, 2013, in Russian.
[2] N. N. Kalitkin, A. B. Alshin, E. A. Alshina, and V. B. Rogov, Computations with Quasi-Uniform Grids [Vychisleniya na kvaziravnomernykh setkakh]. Moscow: Fizmatlit, 2005, in Russian.
[3] 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.
[4] 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.
[5] 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.
[6] T. Sag and G. Szekeres, "Numerical evaluation of high-dimensional integrals," Math. Comp., vol. 18, pp. 245-253, 1964. DOI: 10.1090/ S0025-5718-1964-0165689-X.
[7] A. Sidi, "Numerical evaluation of high-dimensional integrals," International Series Numer. Math., vol. 112, pp. 359-373, 1993. DOI: 10.1007/ 978-3-0348-6338-4_27.
[8] M. Iri, S. Moriguti, and Y. Takasawa, "On a certain quadrature formula," International Series Numer. Math., vol. 17, pp. 3-20, 1987. DOI: 10. 1016/0377-0427(87)90034-3.
[9] M. Mori, "An IMT-Type Double Exponential Formula for Numerical Integration," Publ. Res. Inst . Math. Sci. Kyoto Univ., vol. 14, no. 3, pp. 713-729, 1978. DOI: 10.2977/prims/1195188835.
[10] 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.
[11] 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, M. A. Tintul, V. S. Khokhlachev, Quadratures with super power
convergence, Discrete and Continuous Models and Applied Computational
Science 31 (2) (2023) 128-138. DOI: 10.22363/2658-4670-2023-31-2-128-138.
Information about the authors:
Belov, Aleksandr A. — Candidate of Physical and Mathematical Sciences, Researcher of Faculty of Physics, M. V. Lomonosov Moscow State University; Assistant Professor of Department of Applied Probability and Informatics of Peoples' Friendship University of Russia named after Patrice Lumumba (RUDN University) (e-mail: [email protected], phone: +7(495)9393310, ORCID: https://orcid.org/0000-0002-0918-9263, ResearcherlD: Q-5064-2016, Scopus Author ID: 57191950560)
Tintul, Maxim A. — Master's Degree Student of Faculty of Physics, M. V. Lomonosov Moscow State University (e-mail: [email protected], phone: +7(495)9393310, ORCID: https://orcid.org/0000-0002-5466-1221)
Khoklachev, Valentin S. — Master's Degree Student of Faculty of Physics, M. V. Lomonosov Moscow State University (e-mail: [email protected], 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-2023-31-2-128-138 EDN: XAUSJA
Квадратуры со сверхстепенной сходимостью
А. А. Белов1'2, М. А. Тинтул1, В. С. Хохлачев1
1 Московский государственный университет им. М. В. Ломоносова, Ленинские горы, д. 1, стр. 2, Москва, 119991, Россия 2 Российский университет дружбы народов, ул. Миклухо-Маклая, д. 6, Москва, 117198, Россия
Аннотация. Вычисление квадратур возникает во многих физических и технических приложениях. В статье предложена замена переменных интегрирования, кардинально повышающая точность формулы средних. Для бесконечно гладких подынтегральных функций закон сходимости становится сверхстепенным. Он существенно быстрее степенного и близок к экспоненциальному. Для подынтегральных функций с ограниченной гладкостью реализуется степенная сходимость с максимально достижимым порядком точности.
Ключевые слова: формула трапеций, экспоненциальная сходимость, оценки точности, асимптотически точные оценки
