Generalizing the trapezoidal rule in the complex plane

Abstract

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.

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,

$$ \frac{1}{2\pi i}{\int}_{\varGamma}f(z) \pi\cot\pi z dz=\left\{ \begin{array}{c} \textrm{sum of resi-}\\ \textrm{dues inside }\varGamma \end{array}\right\} =\sum\limits_{k=0}^{N-1}f(k)=\left\{ \begin{array}{c} \textrm{TR result for}\\ {{\int}_{0}^{N}}f(x)dx \end{array}\right\} . $$

(18)

We next note:

• 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. With z = x + iy, it holds that¹⁶

  $$ \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)

Fig. 7

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

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.

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

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

$$ g(z;\alpha,\upbeta,\gamma)=\pi\left( \alpha\cot(\pi(z+i))+\upbeta\cot(\pi z)+\gamma\cot(\pi(z-i))\right), $$

(20)

we want the approximation

$$ \int f_{P}(z)dz\approx{\int}_{\varGamma}f(z) g(z;\alpha,\upbeta,\gamma)dz=\alpha T_{-}+\upbeta T_{0}+\gamma T_{+} $$

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

$$ \begin{array}{@{}rcl@{}} -i\pi(\alpha+\upbeta+\gamma)-2i\pi e^{+2i\pi x}\left( \alpha e^{-2\pi}+\upbeta e^{0}+\gamma e^{+2\pi}\right) & =-i\pi,\\ +i\pi(\alpha+\upbeta+\gamma)+2i\pi e^{-2i\pi x}\left( \alpha e^{+2\pi}+\upbeta e^{0}+\gamma e^{-2\pi}\right) & =+i\pi, \end{array} $$

implying

$$ \left[\begin{array}{ccc} 1 & 1 & 1\\ e^{-2\pi} & 1 & e^{+2\pi}\\ e^{+2\pi} & 1 & e^{-2\pi} \end{array}\right]\left[\begin{array}{c} \alpha\\ \upbeta\\ \gamma \end{array}\right]=\left[\begin{array}{c} 1\\ 0\\ 0 \end{array}\right], $$

(21)

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).

Fig. 9

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