Generalizing the trapezoidal rule in the complex plane

In computational contexts, analytic functions are often best represented by grid-based function values in the complex plane. For integrating periodic functions, the spectrally accurate trapezoidal rule (TR) then becomes a natural choice, due to both accuracy and simplicity. The two key present observations are (i) the accuracy of TR in the periodic case can be greatly increased (doubling or tripling the number of correct digits) by using function values also along grid lines adjacent to the line of integration and (ii) a recently developed end correction strategy for finite interval integrations applies just as well when using these enhanced TR schemes.

1. If an analytic functions is somewhat costly to evaluate, it may be desirable to pay a one-time cost of evaluating over a grid and then re-use this data for multiple purposes, such as for graphical displays (e.g., [7, 12]) and, as focused on here, for highly accurate contour integrations (assuming some vicinity of the integration paths to be free of singularities). In some cases, the fastest evaluation method for the function is naturally grid based (e.g., the Painlevé functions [8]).

2. The present application is atypical in that Clenshaw-Curtis requires roughly twice as many evaluation points as Gaussian quadrature for matching accuracies; c.f., [10] and the discussion in Section 4.2 of [4].

3. First observed by Poisson in the 1820s [9]

4. Although many papers discuss and give error estimates for the periodic TR, there seem to be few—if any—previous references on utilizing this type of supplementary information for improving the accuracy. A somewhat related procedure was used in [3] for approximating derivatives/Taylor coefficients (c.f., its equations (5)–(9)).

5. Repeated application of Birkhoff-Young type formulas (first introduced in [1]) also produce multi-line TR-like formulas in the periodic case. However, as noted in [4], Section 4.3, their weights become such that they reduce rather than improve the accuracy compared with the regular TR.

6. c.f., [13], Example 5, and [7], page 353.

7. We recall the assumption that grid-based data is more economically available than “fresh” function evaluations. If not, it would be more cost effective to double and triple the number of nodes in the 1-line TR than to apply the 3-line and 5-line TR versions, respectively.

8. Appendix A in [5] sketches three derivations with references to still further ones.

9. As was done in [4]

10. The last parameter in the NSolve statement specifies the precision (number of decimal digits) to be used in this calculation of the coefficients.

11. Advanpix, Multiprecision Computing Toolbox for MATLAB, http://www.advanpix.com/, Advanpix LLC., Yokohama, Japan.

12. Figure 6 in [4] shows similarly the nodes used when integrating along the rectangle in the Cartesian node case.

13. One higher than the number of nodes in the correction stencil.

14. In Section 4.3 of [4], end corrected 1-line TR is compared against Gaussian and Clenshaw-Curtis quadratures for the test problem \({\int \limits }_{-1}^{1}f(z)dz\).

15. When expressed as derivatives of odd orders at the end point(s)

16. After the leading term, both expansions are geometric progressions that are readily summed in closed form, giving \(\pi \cot (\pi (x+iy))\) in their respective regions of convergence.

Discussions with André Weideman are gratefully acknowledged.

Authors and Affiliations

1. Department of Applied Mathematics, University of Colorado, Boulder, CO, 80309, USA

   Bengt Fornberg

Corresponding author

Correspondence to Bengt Fornberg.

Publisher’s note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Appendix. Derivation of multi-line periodic TR formulas by means of contour integration

For notational simplicity, we set h = 1 and let the problem be periodic over 0 ≤ x ≤ N. With the path Γ illustrated in Fig. 7a,

We next note:

1. 1.

   The value of \({\int \limits }_{P}f(z)dz\) along the period line P = [0, N] is unaffected if this path is translated sideways or up/down (assuming f(z) is singularity-free in the region of interest), and

2. 2.

   With z = x + iy, it holds that^{Footnote 16}

   $$ \pi \cot\pi z = \left\{ \begin{array}{ll} -i\pi \left( 1+2e^{+2i\pi x}e^{-2\pi y} + 2e^{+4i\pi x}e^{-4\pi y}+2e^{+6i\pi x}e^{-6\pi y} + \ldots\right) & \textrm{if }y > 0\\ +i\pi \left( 1+2e^{-2i\pi x}e^{+2\pi y} + 2e^{-4i\pi x}e^{+4\pi y}+2e^{-6i\pi x}e^{+6\pi y} + \ldots\right) & \textrm{if }y < 0 \end{array}\right. . $$

   (19)

Contours and pole locations in the contour integration-based derivations of the periodic 1-line and 3-line TR (trapezoidal rule) methods

Full size image

Equation (18) tells that TR would give the exact value for \({\int \limits } f_{P}(z)dz\) if it was not for \(\pi \cot \pi z\) approaching ∓ iπ only exponentially fast (19) rather than instantaneously when y deviates from zero, c.f., Fig. 8.

The imaginary part of \(\pi \cot \pi z\), displayed from a viewpoint in the 4^th quadrant. This figure is the same as Figure 5.25 (b) in [7]. Part (a) there illustrates the real part converging to zero in both half-planes, and part (c) illustrates the function’s magnitude and phase angle

Full size image

The different multi-line TR formulas follow now in the same manner as described below for the 3-line Cartesian case. Consider for this case the contour integration illustrated in Fig. 7b, with the three rows of poles given the weights α, β, and γ, respectively. Defining

we want the approximation

to be as accurate as possible. From the argument above, this happens if g(z; α,β, γ) approaches as fast as possible − iπ for y increasing past y = 1 and + iπ for y decreasing below y = − 1. Given (19), this occurs when

implying

with the solution \(\left \{ \alpha ,\upbeta ,\gamma \right \} =\frac {1}{(2\sinh \pi )^{2}}\left \{ -1,2\cosh 2\pi ,-1\right \} \), matching (4). The system (21) is equivalent to the one previously solved when obtaining (4) from (1), (2), and (3).

Figure 8 illustrates \(\text {Im}(\pi \cot \pi z)\), i.e., Im(g(z; α,β, γ)) in the original 1-line TR case of {α, β, γ} = {0, 1, 0}. With the 3-line TR coefficients {α,β, γ}≈{− 0.001874, 1.00375,− 0.001874}, the residues at the poles along the real axis are very slightly increased, and there will be very small residue poles also along y = ± 1. Figure 9a shows exactly the same data as Fig. 8, but limited in the y-direction to [0.7, 1.6]. The counterpart illustration for the 3-line TR is seen in Fig. 9b. The row of poles along y = 1 have negative residues, making them visually appear as if they were shifted in the x-direction. For increasing y, these poles perfectly cancel out the dominant Fourier mode in the x-direction (as seen along the left edge in the two plots).

Illustrations of Im(g(z; α, β, γ)) for the 1-line and 3-line TR formulations, showing how the latter more rapidly approaches − π for increasing y (since the dominant Fourier mode in the x-direction has been eliminated). Part a shows the same data as Fig. 8, restricted to 0.7 ≤ y ≤ 1.6

Full size image

Reprints and permissions