Abstract
Ability to evaluate contour integrals is central to both the theory and the utilization of analytic functions. We present here a complex plane realization of the Euler–Maclaurin formula that includes weights also at some grid points adjacent to each end of a line segment (made up of equispaced grid points, along which we use the trapezoidal rule). For example, with a |$5\times 5$| ‘correction stencil’ (with weights about two orders of magnitude smaller than those of the trapezoidal rule), the accuracy is increased from |$2$|nd to |$26$|th order.
1. Introduction
With increasing use of computational methods in the study and the utilization of analytic functions, the need arises for numerically evaluating contour integrals both along open curve segments and around closed paths. If the analytic function is fast to evaluate at arbitrary points, Gaussian quadrature is a viable approach for integrating along straight line segments. This paper focuses on the case when the cost of function evaluations is nontrivial but values are available on an equispaced grid.1 Assuming the line segments (including their endpoints) to be free of singularities we show here how integrals along them can be evaluated very accurately and at very low cost.
One approach for integrating along a finite part of a grid line (i.e., with equispaced data) would be to use an enhanced Gregory-type formula as described in Fornberg (2020) (based on ideas in Fornberg & Reeger, 2019). These formulas are numerically stable (all weights positive) also for high orders of accuracy. The recommended |$10$|th-order-accurate schemes given in Fornberg (2020) (i.e., with error |$\mathcal{O}(h^{10})$| where |$h$| is the node separation) require, however, at least 11 nodes along each straight line segment.
The approach introduced here permits very high orders of accuracy while also removing the constraint on the segment length. In the present context (of having the analytic function available not only at grid points along each line segment but also at some grid points surrounding the end of each line segment), use of a |$3\times 3$| Cartesian grid ‘correction stencil’ brings the standard trapezoidal rule (TR) error from |$\mathcal{O}(h^{2})$| to |$\mathcal{O}(h^{10})$|, while a |$5\times 5$| stencil suffices for |$\mathcal{O}(h^{26})$| (which is a much higher order of accuracy than is reached in most grid-based contexts). Figures 1(a,b) illustrate end-of-segment stencils (with green circles), while Fig. 2(a,b) show hexagonal grid counterparts providing accuracies |$\mathcal{O}(h^{8})$| and |$\mathcal{O}(h^{20})$|, respectively.2
Two examples of Cartesian grid stencils for correcting the TR, which is here applied along a horizontal path (red dots).
Two examples of hexagonal grid stencils for correcting the TR.
On Cartesian grids one would typically follow integration paths that are made up of straight line segments with |$90^\circ$| angles from segment to segment, while for hexagonal grids, the angles would be either |$60^\circ$| or |$120^\circ$|.
2. The Euler–Maclaurin formula
3. Finite difference formulas in the complex plane
In view of Cauchy’s integral formula |$f^{(k)}(z)=\frac{k!}{2\pi \textrm{i}}\oint _{\varGamma }\frac{f(\xi )}{(\xi -z)^{k+1}}\,\textrm{d}\xi $| it is natural to approximate the derivatives in (2.1) with finite difference stencils that extend about an equal distance in all directions from the origin.3
3.1 Cartesian grid
3.1.1 |$3\times 3$| stencils
The description below for the |$3\times 3$| stencil case of Fig. 1(a) generalizes immediately to the remaining cases of Figs 1 and 2.
The |$3\times 3$| Cartesian correction case: weights for the leading odd-order derivatives and for their combination according to the EM formula (2.1)
| Constants in the stencil entries | Order of derivative | EM combination | |||
|---|---|---|---|---|---|
| 1 | 3 | 5 | 7 | ||
| |$c_{1}$| | |$8$| | |$16$| | |$2$| | |$4$| | |$6044$| |
| |$c_{2}$| | |$1$| | |$-1$| | |$-1$| | |$1$| | |$821$| |
| |$c_{3}$| | |$-1$| | |$1$| | |$1$| | |$1$| | |$-779$| |
| |$c_{4}$| | |$-8$| | |$16$| | |$-2$| | |$-4$| | |$-7556$| |
$$\begin{array}{c} \textrm{All entries to be multiplied by} \end{array}$$ | |$\dfrac{1}{40h}$| | |$\dfrac{3}{40h^{3}}$| | |$\dfrac{3}{h^{5}}$| | |$\dfrac{63}{h^{7}}$| | |$\dfrac{h}{403200}$| |
| Constants in the stencil entries | Order of derivative | EM combination | |||
|---|---|---|---|---|---|
| 1 | 3 | 5 | 7 | ||
| |$c_{1}$| | |$8$| | |$16$| | |$2$| | |$4$| | |$6044$| |
| |$c_{2}$| | |$1$| | |$-1$| | |$-1$| | |$1$| | |$821$| |
| |$c_{3}$| | |$-1$| | |$1$| | |$1$| | |$1$| | |$-779$| |
| |$c_{4}$| | |$-8$| | |$16$| | |$-2$| | |$-4$| | |$-7556$| |
$$\begin{array}{c} \textrm{All entries to be multiplied by} \end{array}$$ | |$\dfrac{1}{40h}$| | |$\dfrac{3}{40h^{3}}$| | |$\dfrac{3}{h^{5}}$| | |$\dfrac{63}{h^{7}}$| | |$\dfrac{h}{403200}$| |
The |$3\times 3$| Cartesian correction case: weights for the leading odd-order derivatives and for their combination according to the EM formula (2.1)
| Constants in the stencil entries | Order of derivative | EM combination | |||
|---|---|---|---|---|---|
| 1 | 3 | 5 | 7 | ||
| |$c_{1}$| | |$8$| | |$16$| | |$2$| | |$4$| | |$6044$| |
| |$c_{2}$| | |$1$| | |$-1$| | |$-1$| | |$1$| | |$821$| |
| |$c_{3}$| | |$-1$| | |$1$| | |$1$| | |$1$| | |$-779$| |
| |$c_{4}$| | |$-8$| | |$16$| | |$-2$| | |$-4$| | |$-7556$| |
$$\begin{array}{c} \textrm{All entries to be multiplied by} \end{array}$$ | |$\dfrac{1}{40h}$| | |$\dfrac{3}{40h^{3}}$| | |$\dfrac{3}{h^{5}}$| | |$\dfrac{63}{h^{7}}$| | |$\dfrac{h}{403200}$| |
| Constants in the stencil entries | Order of derivative | EM combination | |||
|---|---|---|---|---|---|
| 1 | 3 | 5 | 7 | ||
| |$c_{1}$| | |$8$| | |$16$| | |$2$| | |$4$| | |$6044$| |
| |$c_{2}$| | |$1$| | |$-1$| | |$-1$| | |$1$| | |$821$| |
| |$c_{3}$| | |$-1$| | |$1$| | |$1$| | |$1$| | |$-779$| |
| |$c_{4}$| | |$-8$| | |$16$| | |$-2$| | |$-4$| | |$-7556$| |
$$\begin{array}{c} \textrm{All entries to be multiplied by} \end{array}$$ | |$\dfrac{1}{40h}$| | |$\dfrac{3}{40h^{3}}$| | |$\dfrac{3}{h^{5}}$| | |$\dfrac{63}{h^{7}}$| | |$\dfrac{h}{403200}$| |
These |$3\times 3$| stencil corrections have increased the TR accuracy from |$\mathcal{O}(h^{2})$| to |$\mathcal{O}(h^{10})$|.
If the grid line that we integrate along is not directed straight to the right, i.e., the forward step |$h$| from each data point to the next is not positive, we need only to use the actual (negative or complex) |$h$| in (3.2). If we consider two segments joined end to end in the same direction, the corrections shown in (3.2) cancel each other, as ideally should be the case (since there then is no path corner present). These same observations apply in all further stencil cases we will consider (larger stencils, hexagonal nodes, etc.).
3.1.2 |$5\times 5$| stencils
With this correction stencil (having 24 nontrivial entries) the accuracy has improved from |$\mathcal{O}(h^{2})$| for TR to |$\mathcal{O}(h^{26})$|. Figure 3 illustrates the magnitude of the quadrature weights when this correction stencil is applied at the left end of a segment.
Illustration of the magnitudes of the weights in the |$5\times 5$| case schematically shown in Fig. 1(b). Entries belonging to the TR formula are shown in red, to the correction stencil in green and the overlapping entries in blue. Compared to the TR weights (1/2 at each end and 1 otherwise) the correction stencil weights are two (or more) orders of magnitude smaller.
3.2 Hexagonal grid
3.2.1 7-node stencil
3.2.2 19-node stencil
4. Numerical verification
Magnitude and phase angle display of the test function (4.1) with the two integration paths marked over the magnitude surface.
4.1 Comparison between present grid-based methods
We consider two different closed paths both illustrated in Fig. 4:
the rectangular path with corners |$z=\{1,1+\,\mathrm{i},-1+\,\mathrm{i},-1\}$| (using a Cartesian grid) and
the equilateral triangular path with corners |$z=\{1,\sqrt{3}\,\mathrm{i},-1\}$| (using a hexagonal grid),
and integrate around these in the positive direction. In either case the analytic result is |$\int _{C}f(z)\,\textrm{d}z=4\pi \,\mathrm{i}$| (since the only singularity inside the paths is the residue-2 pole at |$z=0.4\:(1+\,\mathrm{i})$|). Figure 6 shows how the numerical errors in these quadratures decrease with the step size |$h$| on the equispaced grids.5 For |$h$| relatively large the exponentially decaying error from the TR dominates.6 There is then a clear transition point after which the algebraic order errors from the ends of the segments become dominant. The slopes of these straight lines confirm the algebraic convergence orders quoted above. Figure 5 illustrates in a different way the quadrature accuracy by showing a grid density that is sufficient in this test case for full double precision (|$10^{-16}$|) quadrature accuracy.
Real and imaginary parts of the same function |$f(z)$| as illustrated previously in Fig. 4, here shown with a sufficient grid density for the Cartesian |$5\times 5$| node method to return the integral to full double precision |$10^{-16}$| accuracy. As before, TR points are marked in red and correction stencil points in green.
Log-log plot of the absolute error when using the four different quadrature implementations.
4.2 Comparisons against Gaussian and Clenshaw–Curtis quadrature
This paper is about efficient contour integration when function values are already available on an equispaced grid (making the cost of the present quadrature methods negligible). It can nevertheless be of interest to compare against standard very-high-order methods, such as Gaussian and Clenshaw–Curtis quadratures (below abbreviated GQ and CC, respectively) in the case when no grid data is available, and all required function values need to be computed. The two methods GQ and CC are surveyed in Trefethen (2008) (with brief MATLAB codes provided).
As a test problem we now consider |$\int _{-1}^{1}f(z)\,\textrm{d}z$| with the same test function (4.1) as above, giving the result shown in Fig. 7. We show here the error as a function of the total number of function evaluations needed, with those for the Cartesian 3|$\times $|3 and 5|$\times $|5 nodes cases being entirely grid based, whereas the evaluation points for both GQ and CC are the special ones needed for these particular methods.7 The difference in the initial exponential convergence rates between ‘Cartesian |$5\times 5$| nodes’ and ‘GQ’ seen in Fig. 7 can be attributed to GQ having its node density ‘depleted’ by a factor |$2/\pi $| near the interval center. The Cartesian grid method can also in some sense be said to have a depletion of nodes along the segment, as it uses additional nodes near the interval ends.
Comparisons between present Cartesian grid-based methods and the Gaussian and Clenshaw–Curtis methods in the case when no grid data is available, and all function evaluations need to be carried out just for the purposes of quadrature.
For accuracy levels around |$10^{-10}$| the Cartesian |$3\times 3$| node approach is seen to compete well. If one wishes to reach either double (|$10^{-16}$|) or quad (|$10^{-33}$|) precision (dotted lines in Fig. 7) we can note that CC is (in this special case) not competitive. This situation is in other applications not typical. For example the second line in the abstract of Trefethen (2008) notes that ‘experiments show that the supposed factor-of-2 advantage of Gauss quadrature is rarely realized.’ Analytic functions are in this context ‘exceptional’, with the factor-of-2 difference very apparent. Thorough discussions on CC and GQ can be found in several additional recent studies, e.g., Weideman & Trefethen (2007), Hale & Trefethen (2008).
The very-high-order GQ methods’ nodes and weights can be precomputed once and for all, so that does not necessitate any recurring cost. However, even in this case when function evaluations need to be done for quadrature purposes only, the present methods (based on TR with end corrections) are seen to be quite competitive with GQ.
4.3 Some comments on Birkhoff–Young formulas
The concept behind (4.2) (and its generalizations) is quite different from the present TR end correction approach. Integrating over a longer interval requires these formulas to be repeated every |$2h$| along the line of integration (i.e., they are not designed to be interval end corrections). Later enhancements described in the literature focus mostly on increasing the accuracy by also using non-equispaced nodes. For an up-to-date overview see Milovanović (2017).
Figure 8 compares the Cartesian node cases (shown before as blue curves in Fig. 5) against the BY 5-point and 9-point schemes. For |$h\rightarrow 0$| (to the right in the figure) the slopes will reflect the schemes’ algebraic order of accuracy (for both approaches largely a matter of the chosen stencil size). A more fundamental difference arises if errors along the rectangle sides dominate those near the corners (in this figure seen to the left for |$h$| relatively large). Along these sides the BY schemes lose several orders of magnitude in accuracy, since there they do not reduce to the regular TR method.8 Whenever the number of nodes along a side exceeds 11, the Cartesian |$5\times 5$| node scheme uses fewer function values than either of the two BY 5-point and 9-point schemes. The additional nodes the BY schemes involve along the straight sides of the integration path reduce the accuracy rather than improve it.
Accuracy comparison for the Cartesian grid test case of integrating around a rectangle. Comparison between repeated application of two BY formulas and the proposed end correction methods.
5. Conclusions
For numerical quadrature based on equispaced data many contemporary textbooks focus on the trapezoidal rule and use then the Newton–Cotes approach to reach higher orders. The even more classical Gregory methods9 were commonly described in numerical texts up to about 50 years ago (e.g., Fröberg, 1965), but are seldom described nowadays, in spite of being advantageous in several key regards. It was noted in Fornberg & Reeger (2019) and Fornberg (2020) that the Gregory approach furthermore can be made to ‘bypass’ the Runge phenomenon and gives positive quadrature weights also for quite high accuracy orders (to |$20$|th order and beyond). In the present study we have utilized the further assumptions that (i) the integrands are analytic functions, and (ii) function values are available also at some equispaced grid points surrounding each end of the straight line segments (rather than only at equispaced grid points within the line segment). With these assumptions we have here demonstrated that the EM formula offers excellent opportunities for very-high-accuracy calculations in a way that does not seem to have been previously utilized. The correction concept described here applies equally well to Cartesian as to hexagonal equispaced grids. It can readily be generalized to still higher orders of accuracy than we have given closed-form weights for here.
Footnotes
The function may have already been evaluated on a grid for graphical display purposes (Wegert, 2012; Fornberg & Piret, 2020), or the fastest evaluation method for the function might itself be grid based (e.g., the Painlevé functions; see Fornberg & Weideman, 2011). Yet another application case would be if integrals are needed along a large number of different paths.
As we will see later the correction will always be zero at the stencil center point. Hence, this point is not marked by a green ring in these figures.
Exactly equal distance (allowing a circular contour |$\varGamma $|, and thus an alternate way of approximating derivatives) is possible only for very small grid-based stencils—in this study only the 7-node stencil on a hexagonal grid.
Or equivalently (in view of the coefficient symmetries) the correction stencil at the segment end is the same as the one at the segment start, rotated by |$180^{\circ }$| around its center.
Computed using the Advanpix: Multiprecision Computing Toolbox for MATLAB; Advanpix LLC, Yokohama, Japan.
The somewhat slower exponential rate in the hexagonal path case is merely a reflection of this (triangular) path passing somewhat closer to the nearest pole, as seen in Fig. 4. Better than algebraic convergence of TR when not affected by end-of-segment effects follows from the EM formula. Exponential convergence of the TR for periodic problems was analyzed by Poisson in the 1820s (Poisson, 1827). For a modern survey of the TR in this case see Trefethen & Weideman (2014).
In the node count for the Cartesian 3|$\times $|3 and 5|$\times $|5 nodes cases we include in the number of function evaluations both the TR nodes along the line as well as one set of end correction nodes, since these latter sets typically will be used twice (for each line segment sharing a path corner).
Creating a BY-type scheme, for example, with a |$3\times 4$| weight stencil approximating the integral just between its two interior points, and repeating this every |$h$| (rather than every |$2h$|) sideways provides a better approximation along the line interior (however, still falling short of the simpler TR scheme).
First described in a letter written by James Gregory in 1670 (Rigaud, 1841). Although likely first developed (by Cotes) around 1707, the first mention in the literature of the Newton--Cotes formulas dates to 1722.
References
