A fast simple algorithm for computing the potential of charges on a line

TL;DR

This paper introduces a rapid and straightforward algorithm to efficiently compute electrostatic potentials generated by charges positioned on a line, significantly simplifying calculations in relevant physics applications.

Contribution

The paper presents a novel, simple, and fast algorithm specifically designed for calculating potentials of line charges, filling a gap in computational physics methods.

Findings

Algorithm significantly reduces computation time.

Numerical tests demonstrate high accuracy.

Applicable to various charge distributions on a line.

Abstract

We present a fast method for evaluating expressions of the form $$ u_j = \sum_{i = 1,i \not = j}^n \frac{\alpha_i}{x_i - x_j}, \quad \text{for} \quad j = 1,\ldots,n, $$ where $\alpha_i$ are real numbers, and $x_i$ are points in a compact interval of $\mathbb{R}$. This expression can be viewed as representing the electrostatic potential generated by charges on a line in $\mathbb{R}^3$. While fast algorithms for computing the electrostatic potential of general distributions of charges in $\mathbb{R}^3$ exist, in a number of situations in computational physics it is useful to have a simple and extremely fast method for evaluating the potential of charges on a line; we present such a method in this paper, and report numerical results for several examples.

Tables

Table 1. Table 1. Key for column labels of Tables 2 , 3 , and 4 .

Label	Definition
$n$	number of points
$t_w$	time of algorithm of §3 without precomputation in seconds
$t_p$	time of precomputing exponentials for algorithm of §3 in seconds
$t_u$	time of algorithm of §3 using precomputed exponentials in seconds
$t_d$	time of direct evaluation in seconds
${ϵ}_r$	maximum absolute relative difference defined in (32)
$t_f$	time of FFT using precomputed exponentials (for time comparison only)

Table 2. Table 2. Numerical results for uniformly random points in [ 1 , 10 ] 1 10 [1,10] .

Table 3. Table 3. Numerical results for Chebyshev nodes on [ − 1 , 1 ] 1 1 [-1,1] .

Table 4. Table 4. Time comparison with FFT.

Equations

$$u_j=\sum _{i=1,i\ne j}^n\frac{{\alpha }_i}{x_i-x_j},\text{for}j=1,\dots ,n,$$

$$\frac{1}{r}={\int }_0^{\infty }{e}^{-rt}dt.$$

$$\frac{1}{r}\approx \sum _{j=1}^m{w}_j{e}^{-r{t}_j},$$

$$\frac{1}{r}\mapsto \frac{1}{r+r^{\prime }}$$

$$\frac{1}{r+r^{\prime }}\approx \sum _{j=1}^m{w}_j{e}^{-(r+r^{\prime })t_j}=\sum _{j=1}^m{w}_j{e}^{-r^{\prime }{t}_j}{e}^{-r{t}_j}.$$

$$Hf(x)=\underset{\epsilon \to 0}{\lim }\frac{1}{\pi }{\int }_{|x-y|\ge \epsilon }\frac{f(y)}{y-x}dy.$$

$$P_{m}f(x):={\int }_a^b\sum _{k=0}^m\frac{p_k(x)p_k(y)}{h_k}f(y)w(y)dy,$$

$$u_j:=\sum _{i=1}^n\sum _{k=0}^m\frac{p_k(x_j)p_k(x_i)}{h_k}f(x_i)w(x_i),\text{for}j=1,\dots ,n.$$

$$\sum _{k=0}^m\frac{p_k(x)p_k(y)}{h_k}=\frac{1}{h_m}\frac{p_{m+1}(x)p_m(y)-p_m(x)p_{m+1}(y)}{x-y},$$

$$u_j=\frac{1}{h_m}\left(f(x_j)+\sum _{i=1,i\ne j}^m\frac{p_{m+1}(x_j)p_m(x_i)-p_m(x_j)p_{m+1}(x_i)}{x_j-x_i}f(x_i)w(x_i)\right),$$

$$\sum _{k=1}^{\infty }2(\nu +k)J_{\nu +k}(w)J_{v+k}(z)=\frac{wz}{w-z}\left(J_{\nu +1}(w)J_{\nu }(z)-J_{\nu }(w)J_{\nu +1}(z)\right),$$

$$v_j=\sum _{i=1}^n\frac{{\alpha }_i}{x_i-y_j},\text{for}j=1,\dots ,m,$$

$$u_j:=\sum _{i=1,i\ne j}^n\frac{{\alpha }_i}{x_i-x_j},\text{for}j=1,\dots ,n.$$

$$\frac{\left|{\tilde{u}}_j-u_j\right|}{{\sum }_{i=1}^n|{\alpha }_i|}\le \epsilon ,\text{for}j=1,\dots ,n$$

$$N_{\delta }:=\sum _{j=1}^n\#\{x_i:|x_j-x_i|<\delta (b-a)\}.$$

$$x_j=2^{-n}\left(\frac{j-1}{n-1}\right),\text{and}{x}_{n+j}=1+2^{-n}\left(\frac{j-n}{n-1}\right),\text{for}j=1,\dots ,n.$$

$$n\log ({\delta }^{-1})+N_{\delta }=\Omega (n^2).$$

$${\tilde{u}}_j^{+}\approx \sum _{i=1}^{j-1}\frac{{\alpha }_i}{x_i-x_j},\text{for}j=1,\dots ,n,$$

$${\tilde{u}}_j^{-}\approx \sum _{i=j+1}^n\frac{{\alpha }_i}{x_i-x_j},\text{for}j=1,\dots ,n.$$

$${\tilde{u}}_j\approx \sum _{i=1,i\ne j}^n\frac{{\alpha }_i}{x_i-x_j},\text{for}j=1,\dots ,n.$$

$${\tilde{u}}_j^{+}\approx \sum _{i=1}^{j-1}\frac{{\alpha }_i}{x_i-x_j},\text{for}j=1,\dots ,n.$$

$$\frac{1}{r}\approx \sum _{i=1}^m{w}_i{e}^{-r{t}_i}\text{for}r\in [\delta (b-a),b-a],$$

$${\tilde{u}}_j^{+}\approx \sum _{i=1}^{j_0}\frac{{\alpha }_i}{x_i-x_j}+\sum _{i=j_0+1}^{j-1}\frac{{\alpha }_i}{x_i-x_j},$$

$${\tilde{u}}_j^{+}\approx -\sum _{i=1}^{j_0}\sum _{k=1}^m{w}_k{\alpha }_i{e}^{-(x_j-x_i)t_k}+\sum _{i=j_0+1}^{j-1}\frac{{\alpha }_i}{x_i-x_j}.$$

$$g_k(j_0)=\sum _{i=1}^{j_0}{\alpha }_i{e}^{-(x_{j_0}-x_i)t_k},\text{for}k=1,\dots ,m,$$

$${\tilde{u}}_j^{+}\approx -\sum _{k=1}^m{w}_k{g}_k(j_0)e^{-(x_j-x_{j_0})t_k}+\sum _{i=j_0+1}^{j-1}\frac{{\alpha }_i}{x_i-x_j}.$$

$$g_k(j_0+1)={\alpha }_{j_0}+e^{-(x_{j_0+1}-x_{j_0})t_k}{g}_k(j_0),\text{for}k=1,\dots ,m.$$

$$\left|\frac{1}{r}-\sum _{j=1}^m{w}_j{e}^{-r{t}_j}\right|\le \epsilon \text{for}r\in [\delta (b-a),b-a].$$

$$P(x_j)=f(x_j),\text{for}j=1,\dots ,n.$$

Full text

[TABLE]

where $\alpha_{i}$ are real numbers, and $x_{i}$ are points in a compact interval of $\mathbb{R}$. This expression can be viewed as representing the electrostatic potential generated by charges on a line in $\mathbb{R}^{3}$. While fast algorithms for computing the electrostatic potential of general distributions of charges in $\mathbb{R}^{3}$ exist, in a number of situations in computational physics it is useful to have a simple and extremely fast method for evaluating the potential of charges on a line; we present such a method in this paper, and report numerical results for several examples.

N.F.M. was supported in part by NSF DMS-1903015

V.R. was supported in part by AFOSR FA9550-16-1-0175 and ONR N00014-14-1-0797.

1. Introduction and motivation

1.1. Introduction

In this paper, we describe a simple fast algorithm for evaluating expressions of the form

[TABLE]

where $\alpha_{i}$ are real numbers, and $x_{i}$ are points in a compact interval of $\mathbb{R}$. This expression can be viewed as representing the electrostatic potential generated by charges on a line in $\mathbb{R}^{3}$. We remark that fast algorithms for computing the electrostatic potential generated by general distributions of charges in $\mathbb{R}^{3}$ exist, see for example the Fast Multipole Method [9] whose relation to the method presented in this paper is discussed in §1.2. However, in a number of situations in computational physics it is useful to have a simple and extremely fast method for evaluating the potential of charges on a line; we present such a method in this paper. Under mild assumptions the presented method involves $\mathcal{O}(n\log n)$ operations and has a small constant. The method is based on writing the potential $1/r$ as

[TABLE]

We show that there exists a small set of quadrature nodes $t_{1},\ldots,t_{m}$ and weights $w_{1},\ldots,w_{m}$ such that for a large range of values of $r$ we have

[TABLE]

see Lemma 4.5, which is a quantitative version of (2). Numerically the nodes $t_{1},\ldots,t_{m}$ and weights $w_{1},\ldots,w_{m}$ are computed using a procedure for constructing generalized Gaussian quadratures, see §5.2. An advantage of representing $1/r$ as a sum of exponentials is that the translation operator

[TABLE]

can be computed by taking an inner product of the weights $(w_{1},\ldots,w_{m})$ with a diagonal transformation of the vector $(e^{-rt_{1}},\ldots,e^{-rt_{m}})$. Indeed, we have

[TABLE]

The algorithm described in §3 leverages the existence of this diagonal translation operator to efficiently evaluate (1).

1.2. Relation to past work

We emphasize that fast algorithms for computing the potential generated by arbitrary distributions of charges in $\mathbb{R}^{3}$ exist. An example of such an algorithm is the Fast Multipole Method that was introduced by [9] and has been extended by several authors including [7, 10, 16]. In this paper, we present a simple scheme for the special case where the charges are on a line, which occurs in a number of numerical calcuations, see 1.3. The presented scheme has a much smaller runtime constant compared to general methods, and is based on the diagonal form (4) of the translation operator (3). The idea of using the diagonal form of this translation operator to accelerate numerical computations has been studied by several authors; in particular, the diagonal form is used in algorithms by Dutt, Gu and Rokhlin [6], and Yavin and Rokhlin [22] and was subsequently studied in detail by Beylkin and Monzón [1, 2].

The current paper improves upon these past works by taking advantage of robust generalized Gaussian quadrature codes [4] that were not previously available; these codes construct a quadrature rule that is exact for functions in the linear span of a given Chebyshev system, and can be viewed as a constructive version of Lemma 4.2 of Kreĭn [13]. The resulting fast algorithm presented in §3 simplifies past approaches, and has a small runtime constant; in particular, its computational cost is similar to the computational cost of $5$-$10$ Fast Fourier Transforms on data of a similar length, see 5.

1.3. Motivation

Expressions of the form (1) appear in a number of situations in computational physics. In particular, such expressions arise in connection with the Hilbert Transform

[TABLE]

For example, the computation of the projection $P_{m}f$ of a function $f$ onto the first $m+1$ functions in a family of orthogonal polynomials can be reduced to an expression of the form (1) by using the Christoffel–Darboux formula, which is related to the Hilbert transform; we detail the reduction of $P_{m}f$ to an expression of the form (1) in the following.

Let $\{p_{k}\}_{k=0}^{\infty}$ be a family of monic polynomials that are orthogonal with respect to the weight $w(x)\geq 0$ on $(a,b)\subseteq\mathbb{R}$. Consider the projection operator

[TABLE]

where $h_{k}:=\int_{a}^{b}p_{k}(x)^{2}w(x)dx$. Let $x_{1},\ldots,x_{n}$ and $w_{1},\ldots,w_{n}$ be the $n>m/2$ point Gaussian quadrature nodes and weights associated with $\{p_{k}\}_{k=0}^{\infty}$, and set

[TABLE]

By construction the polynomial that interpolates the values $u_{1},\ldots,u_{n}$ at the points $x_{1},\ldots,x_{n}$ will accurately approximate $P_{m}f$ on $(a,b)$ when $f$ is sufficiently smooth, see for example §7.4.6 of Dahlquist and Björck [5]. Directly evaluating (5) would require $\Omega(n^{2})$ operations. In contrast, the algorithm of this paper together with the Christoffel–Darboux Formula can be used to evaluate (5) in $\mathcal{O}(n\log n)$ operations. The Christoffel-Darboux formula states that

[TABLE]

see §18.2(v) of [17]. Using (6) to rewrite (5) yields

[TABLE]

where we have used the fact that the diagonal term of the double summation is equal to $f(x_{j})/h_{m}$. The summation in (7) can be rearranged into two expressions of the form (1), and thus the method of this paper can be used to compute a representation of $P_{m}f$ in $\mathcal{O}(n\log n)$ operations.

Remark 1.1**.**

Analogs of the Christoffel–Darboux formula hold for many other families of functions; for example, if $J_{\nu}(w)$ is a Bessel function of the first kind, then we have

[TABLE]

see [21]. This formula can be used to write a projection operator related to Bessel functions in an analogous form to (7), and the algorithm of this paper can be similarly applied

Remark 1.2**.**

A simple modification of the algorithm presented in this paper can be used to evaluate more general expressions of the form

[TABLE]

where $x_{1},\ldots,x_{n}$ are source points, and $y_{1},\ldots,y_{m}$ are target points. For simplicity, this paper focuses on the case where the source and target points are the same, which is the case in the projection application described above.

2. Main result

2.1. Main result

Our principle analytical result is the following theorem, which provides precise accuracy and computational complexity guarantees for the algorithm presented in this paper, which is detailed in §3.

Theorem 2.1**.**

Let $x_{1}<\ldots<x_{n}\in[a,b]$ and $\alpha_{1},\ldots,\alpha_{n}\in\mathbb{R}$ be given. Set

[TABLE]

Given $\delta>0$ and $\varepsilon>0$, the algorithm described in §3 computes values $\tilde{u}_{j}$ such

[TABLE]

in $\mathcal{O}\left(n\log(\delta^{-1})\log(\varepsilon^{-1})+N_{\delta}\right)$ operations, where

[TABLE]

The proof of Theorem 2.1 is given in §4. Under typical conditions, the presented algorithm involves $\mathcal{O}(n\log n)$ operations. The following corollary describes a case of interest, where the points $x_{1},\ldots,x_{n}$ are Chebyshev nodes for a compact interval $[a,b]$ (we define Chebyshev nodes in §4.2).

Corollary 2.1**.**

Fix $\varepsilon=10^{-15}$, and let the points $x_{1},\ldots,x_{n}$ be Chebyshev nodes on $[a,b]$. If $\delta=1/n$, then the algorithm of §3 involves $\mathcal{O}(n\log n)$ operations.

The proof of Corollary 2.1 is given in §4.4. The following corollary states that a similar result holds for uniformly random points.

Corollary 2.2**.**

Fix $\varepsilon=10^{-15}$, and suppose that $x_{1},\ldots,x_{n}$ are sampled uniformly at random from $[a,b]$. If $\delta=1/n$, then the algorithm of §3 involves $\mathcal{O}(n\log n)$ operations with high probability.

The proof of Corollary 2.2 is immediate from standard probabilistic estimates. The following remark describes an adversarial configuration of points.

Remark 2.1**.**

Fix $\varepsilon>0$, and let $x_{1},\ldots,x_{2n}$ be a collection of points such that $x_{1},\ldots,x_{n}$ and $x_{n+1},\ldots,x_{2n}$ are evenly spaced in $[0,2^{-n}]$ and $[1-2^{-n},1]$, respectively, that is

[TABLE]

We claim that Theorem 2.1 cannot guarantee a complexity better than $\mathcal{O}(n^{2})$ for this configuration of points. Indeed, if $\delta\geq 2^{-n}$, then $N_{\delta}\geq n^{2}/2$, and if $\delta<2^{-n}$, then $\log_{2}(\delta^{-1})>n$. In either case

[TABLE]

This complexity is indicative of the performance of the algorithm for this point configuration; the reason that the algorithm performs poorly is that structures exist at two different scales. If such a configuration were encountered in practice, it would be possible to modify the algorithm of §3 to also involve two different scales to achieve evaluation in $\mathcal{O}(n\log n)$ operations.

3. Algorithm

3.1. High level summary

The algorithm involves passing over the points $x_{1},\ldots,x_{n}$ twice. First, we pass over the points in ascending order and compute

[TABLE]

and second, we pass over the points in descending order and compute

[TABLE]

Finally, we define $\tilde{u}_{j}:=\tilde{u}_{j}^{+}+\tilde{u}_{j}^{-}$ for $j=1,\ldots,n$ such that

[TABLE]

We call the computation of $\tilde{u}_{1}^{+},\ldots,\tilde{u}_{n}^{+}$ the forward pass of the algorithm, and the computation of $\tilde{u}_{1}^{-},\ldots,\tilde{u}_{n}^{+}$ the backward pass of the algorithm. The forward pass of the algorithm computes the potential generated by all points to the left of a given point, while the backward pass of the algorithm computes the potential generated by all points to the right of a given point. In §3.2 and §3.3 we give an informal and detailed description of the forward pass of the algorithm. The backward pass of the algorithm is identical except it considers the points in reverse order.

3.2. Informal description

In the following, we give an informal description of the forward pass of the algorithm that computes

[TABLE]

Assume that a small set of nodes $t_{1},\ldots,t_{m}$ and weights $w_{1},\ldots,w_{m}$ such that

[TABLE]

where $\delta>0$ is given and fixed. The existence and computation of such nodes and weights is described in §4.4 and §5.2. We divide the sum defining $u_{j}^{+}$ into two parts:

[TABLE]

where j_{0}=\max\big{\{}i:x_{j}-x_{i}>\delta(b-a)\big{\}}. By definition, the points $x_{1},\ldots,x_{j_{0}}$ are all distance at least $\delta(b-a)$ from $x_{j}$. Therefore, by (12)

[TABLE]

If we define

[TABLE]

then it is straightforward to verify that

[TABLE]

Observe that we can update $g_{k}(j_{0})$ to $g_{k}(j_{0}+1)$ using the following formula

[TABLE]

We can now summarize the algorithm for computing $\tilde{u}_{1}^{+},\ldots,\tilde{u}_{n}^{+}$. For each $j$, we compute $\tilde{u}_{j}^{+}$ by the following three steps:

Update $g_{1},\ldots,g_{m}$ as necessary 2. 2.

Use $g_{1},\ldots,g_{m}$ to evaluate the potential from $x_{i}$ such that $x_{j}-x_{i}>\delta(b-a)$ 3. 3.

Directly evaluate the potential from $x_{i}$ such that $0<x_{j}-x_{i}<\delta(b-a)$

By (16), each update of $g_{1},\ldots,g_{m}$ requires $\mathcal{O}(m)$ operations, and we must update $g_{1},\ldots,g_{m}$ at most $n$ times, so we conclude that the total cost of the first step of the algorithm is $\mathcal{O}(nm)$ operations. For each $j=1,\ldots,n$, the second and third step of the algorithm involve $\mathcal{O}(m)$ and $\mathcal{O}(\#\{x_{i}:0<x_{j}-x_{i}<\delta(b-a)\})$ operations, respectively, see (15). It follows that the total cost of the second and third step of the algorithm is $\mathcal{O}(nm+N_{\delta})$ operations, where $N_{\delta}$ is defined in (9). We conclude that $\tilde{u}_{1}^{+},\ldots,\tilde{u}_{n}^{+}$ can be computed in $\mathcal{O}(nm+N_{\delta})$ operations. In §4, we complete the proof of the computational complexity guarantees of Theorem 2.1 by showing that there exist $m=\mathcal{O}(\log(\delta^{-1})\log(\varepsilon^{-1}))$ nodes $t_{1},\ldots,t_{m}$ and weights $w_{1},\ldots,w_{m}$ that satisfy (12), where $\varepsilon>0$ is the approximation error in (12).

3.3. Detailed description

In the following, we give a detailed description of the forward pass of the algorithm that computes $\tilde{u}_{1}^{+},\ldots,\tilde{u}_{n}^{+}$. Suppose that $\delta>0$ and $\varepsilon>0$ are given and fixed. We describe the algorithm under the assumption that we are given quadrature nodes $t_{1},\ldots,t_{m}$ and weights $w_{1},\ldots,w_{m}$ such that

[TABLE]

The existence of such weights and nodes is established in §4.4, and the computation of such nodes and weights is discussed in §5.2. To simplify the description of the algorithm, we assume that $x_{0}=-\infty$ is a placeholder node that does not generate a potential.

Algorithm 3.1**.**

Input: $x_{1}<\cdots<x_{n}\in[a,b]$, $\alpha_{1},\ldots,\alpha_{n}\in\mathbb{R}$. Output: $\tilde{u}_{1}^{+},\ldots,\tilde{u}_{n}^{+}$.

1. 1:

$j_{0}=0$ and $g_{1}=\cdots=g_{m}=0$ 2. 2: 3. 3:

main loop: 4. 4:

for $j=1,\ldots,n$ 5. 5: 6. 6:

update $g_{1},\ldots,g_{m}$ and $j_{0}$: 7. 7:

while $x_{j}-x_{j_{0}+1}>\delta(b-a)$ 8. 8:

for i = 1,…,m 9. 9:

$g_{i}=g_{i}e^{-(x_{j_{0}+1}-x_{j_{0}})t_{i}}+\alpha_{i}$ 10. 10:

end for 11. 11:

$j_{0}=j_{0}+1$ 12. 12:

end while 13. 13: 14. 14:

compute potential from $x_{i}$ such that $x_{i}\leq x_{j_{0}}:$ 15. 15:

$\tilde{u}_{j}^{+}=0$ 16. 16:

for $i=1,\ldots,m$ 17. 17:

$\tilde{u}_{j}^{+}=\tilde{u}_{j}^{+}-w_{i}g_{i}e^{-(x_{j}-x_{j_{0}})t_{i}}$ 18. 18:

end for 19. 19: 20. 20:

compute potential from $x_{i}$ such that $x_{j_{0}+1}\leq x_{i}\leq x_{j-1}$ 21. 21:

for $i=j_{0}+1,\ldots,j-1$ 22. 22:

$\tilde{u}_{j}^{+}=\tilde{u}_{j}^{+}+\alpha_{i}/(x_{i}-x_{j}).$ 23. 23:

end for 24. 24:

end for

Remark 3.1**.**

In some applications, it may be necessary to evaluate an expression of the form (1) for many different weights $\alpha_{1},\ldots,\alpha_{n}$ associated with a fixed set of points $x_{1},\ldots,x_{n}$. For example, in the projection application described in §1.3 the weights $\alpha_{1},\ldots,\alpha_{n}$ correspond to a function that is being projected, while the points $x_{1},\ldots,x_{n}$ are a fixed set of quadrature nodes. In such situations, pre-computing the exponentials $e^{-(x_{j}-x_{j_{0}})t_{i}}$ used in the Algorithm 3.1 will significantly improve the runtime, see §5.1.

4. Proof of Main Result

4.1. Organization

In this section we complete the proof of Theorem 2.1; the section is organized as follows. In §4.2 we give mathematical preliminaries. In §4.3 we state and prove two technical lemmas. In §4.4 we prove Lemma 4.5, which together with the analysis in §3 establishes Theorem 2.1. In §4.5 we prove Corollary 2.1, and Corollary 2.2.

4.2. Preliminaries

Let $a<b\in\mathbb{R}$ and $n\in\mathbb{Z}_{>0}$ be fixed, and suppose that $f:[a,b]\rightarrow\mathbb{R}$, and $x_{1}<\cdots<x_{n}\in[a,b]$ are given. The interpolating polynomial $P$ of the function $f$ at $x_{1},\ldots,x_{n}$ is the unique polynomial of degree at most $n-1$ such that

[TABLE]

This interpolating polynomial $P$ can be explicitly defined by

[TABLE]

where $q_{j}$ is the nodal polynomial for $x_{j}$, that is,

[TABLE]

We say $x_{1},\ldots,x_{n}$ are Chebyshev nodes for the interval $[a,b]$ if

[TABLE]

The following lemma is a classical result in approximation theory. It says that a smooth function on a compact interval is accurately approximated by the interpolating polynomial of the function at Chebyshev nodes, see for example §4.5.2 of Dahlquist and Björck [5].

Lemma 4.1**.**

Let $f\in C^{n}([a,b])$, and $x_{1},\ldots,x_{n}$ be Chebyshev nodes for $[a,b]$. If $P$ is the interpolating polynomial for $f$ at $x_{1},\ldots,x_{n}$, then

[TABLE]

where

[TABLE]

In addition to Lemma 4.1, we require a result about the existence of generalized Gaussian quadratures for Chebyshev systems. In 1866, Gauss [8] established the existence of quadrature nodes $x_{1},\ldots,x_{n}$ and weights $w_{1},\ldots,w_{n}$ for an interval $[a,b]$ such that

[TABLE]

whenever $f(x)$ is a polynomial of degree at most $2n-1$. This result was generalized from polynomials to Chebyshev systems by Kreĭn [13]. A collection of functions $f_{0},\ldots,f_{n}$ on $[a,b]$ is a Chebyshev system if every nonzero generalized polynomial

[TABLE]

has at most $n$ distinct zeros in $[a,b]$. The following result of Kreĭn says that any function in the span of a Chebyshev system of $2n$ functions can be integrated exactly by a quadrature with $n$ nodes and $n$ weights.

Lemma 4.2** (Kreĭn [13]).**

Let $f_{0},\ldots,f_{2n-1}$ be a Chebyshev system of continuous functions on $[a,b]$, and $w:(a,b)\rightarrow\mathbb{R}$ be a continuous positive weight function. Then, there exists unique nodes $x_{1},\ldots,x_{n}$ and weights $w_{1},\ldots,w_{n}$ such that

[TABLE]

whenever $f$ is in the span of $f_{0},\ldots,f_{2n-1}$.

4.3. Technical Lemmas

In this section, we state and prove two technical lemmas that are involved in the proof of Theorem 2.1. We remark that a similar version of Lemma 4.3 appears in [18].

Lemma 4.3**.**

Fix $a>0$ and $t\in[0,\infty)$, and let $r_{1},\ldots,r_{n}$ be Chebyshev nodes for $[a,2a]$. If $P_{t}(r)$ is the interpolating polynomial for $e^{-rt}$ at $r_{1},\ldots,r_{n}$, then

[TABLE]

Proof.

We have

[TABLE]

By writing the derivative of $t^{n}e^{-ta}$ as

[TABLE]

we can deduce that the maximum of $t^{n}e^{-ta}$ occurs at $t=n/a$, that is,

[TABLE]

By (21) and the result of Lemma 4.1, we conclude that

[TABLE]

It remains to show that $2n^{n}e^{-n}\leq n!$. Since $\ln(x)$ is a increasing function, we have

[TABLE]

Exponentiating both sides of this inequality gives $en^{n}e^{-n}\leq n!$, which is a classical inequality related to Stirling’s approximation. This completes the proof. ∎

Lemma 4.4**.**

Suppose that $\varepsilon>0$ and $M>1$ are given. Then, there exists

[TABLE]

values $r_{1},\ldots,r_{m}\in[1,M]$ such that for all $r\in[1,M]$ we have

[TABLE]

for some choice of coefficients $c_{j}(r)$ that depend on $r$.

Proof.

We construct an explicit set of $m:=(\lfloor\log_{2}M\rfloor+1)(\lfloor\log_{4}\varepsilon^{-1}\rfloor+1)$ points and coefficients such that (22) holds. Set $n:=\lfloor\log_{4}\varepsilon^{-1}\rfloor+1$. We define the points $r_{1},\ldots,r_{m}$ by

[TABLE]

for $k=1,\ldots,n$ and $i=0,\ldots,\lfloor\log_{2}M\rfloor$, and define the coefficients $c_{1}(r),\ldots,c_{m}(r)$ by

[TABLE]

for $k=1,\ldots,n$ and $i=0,\ldots,\lfloor\log_{2}M\rfloor$. We claim that

[TABLE]

Indeed, fix $r\in[1,M]$, and let $i_{0}\in\{0,\ldots,\lfloor\log_{2}M\rfloor\}$ be the unique integer such that $r\in[2^{i_{0}},2^{i_{0}+1})$. By the definition of the coefficients, see (24), we have

[TABLE]

We claim that the right hand side of this equation is the interpolating polynomial $P_{t,i_{0}}(r)$ for $e^{-rt}$ at $r_{i_{0}n+k},\ldots,r_{(i_{0}+1)n}$, that is,

[TABLE]

Indeed, see (18) and (19). Since the points $r_{i_{0}n+k},\ldots,r_{(i_{0}+1)n}$ are Chebyshev nodes for the interval $[2^{i_{0}},2^{i_{0}+1}]$, and since $i_{0}$ was chosen such that $r\in[2^{i_{0}},2^{i_{0}+1})$, it follows from Lemma 4.3 that

[TABLE]

Since $n=\lfloor\log_{4}\varepsilon^{-1}\rfloor+1$ the proof is complete. ∎

Remark 4.1**.**

The proof of Lemma 4.4 has the additional consequence that the coefficients $c_{1}(r),\ldots,c_{m}(r)$ in (22) can be chosen such that they satisfy

[TABLE]

Indeed, in (24) the coefficients $c_{j}(r)$ are either equal zero or equal to the nodal polynomial, see (19), for Chebyshev nodes on an interval that contains $r$. The nodal polynomials for Chebyshev nodes on an interval $[a,b]$ are bounded by $\sqrt{2}$ on $[a,b]$, see for example [18]. The fact that $e^{-rt}$ can be approximated as a linear combination of functions $e^{-r_{1}t},\ldots,e^{-r_{m}t}$ with small coefficients means that the approximation of Lemma 4.4 can be used in finite precision environments without any unexpected catastrophic cancellation.

4.4. Completing the proof of Theorem 2.1

Previously in §3.2, we proved that the algorithm of §3 involves $\mathcal{O}(nm+N_{\delta})$ operations. To complete the proof of Theorem 2.1 it remains to show that there exists

[TABLE]

points $t_{1},\ldots,t_{m}$ and weights $w_{1},\ldots,w_{m}$ that satisfy (17); we show the existence of such nodes and weights in the following lemma, and thus complete the proof of Theorem 2.1. The computation of such nodes and weights is described in §5.2.

Lemma 4.5**.**

Fix $a<b\in\mathbb{R}$, and let $\delta>0$ and $\varepsilon>0$ be given. Then, there exists $m=\mathcal{O}(\log(\varepsilon^{-1})\log(\delta^{-1}))$ nodes $t_{1},\ldots,t_{m}$ and weights $w_{1},\ldots,w_{m}$ such that

[TABLE]

Proof.

Fix $a<b\in\mathbb{R}$, and let $\delta,\varepsilon>0$ be given. By the possibility of rescaling $r$, $w_{j}$, and $t_{j}$, we may assume that $b-a=\delta^{-1}$ such that we want to establish (25) for $r\in[1,\delta^{-1}]$. By Lemma 4.4 we can choose $2m=\mathcal{O}(\log(\varepsilon^{-1})\log(\delta^{-1}))$ points $r_{0},\ldots,r_{2m-1}\in[1,\delta^{-1}]$, and coefficients $c_{0}(r),\ldots,c_{2m-1}(r)$ depending on $r$ such that

[TABLE]

The collection of functions $e^{-r_{0}t},\ldots,e^{-r_{2m-1}t}$ form a Chebyshev system of continuous functions on the interval $[0,\log(2\varepsilon^{-1})]$, see for example [12]. Thus, by Lemma 4.2 there exists $m$ quadrature nodes $t_{1},\ldots,t_{m}$ and weights $w_{1},\ldots,w_{m}$ such that

[TABLE]

whenever $f(t)$ is in the span of $e^{-r_{0}t},\ldots,e^{-r_{2m-1}t}$. By the triangle inequality

[TABLE]

Recall that we have assumed $r\in[1,\delta^{-1}]$, in particular, $r\geq 1$ so it follows that

[TABLE]

By (26), the function $e^{-rt}$ can be approximated to error $\varepsilon/(2\log(2\varepsilon^{-1}))$ in the $L^{\infty}$-norm on $[0,\log(2\varepsilon^{-1})]$ by functions in the span of $e^{-r_{0}t},\ldots,e^{-r_{2m-1}t}$. Since our quadrature is exact for these functions, we conclude that

[TABLE]

Combining (27), (28), and (29) completes the proof. ∎

4.5. Proof of Corollary 2.1

In this section, we prove Corollary 2.1, which states that the algorithm of §3 involves $\mathcal{O}(n\log n)$ operations when $x_{1},\ldots,x_{n}$ are Chebyshev nodes, $\varepsilon=10^{-15}$, and $\delta=1/n$.

Proof of Corollary 2.1.

By rescaling the problem we may assume that $[a,b]=[-1,1]$ such that the Chebyshev nodes $x_{1},\ldots,x_{n}$ are given by

[TABLE]

By the result of Theorem 2.1, it suffices to show that $N_{\delta}=\mathcal{O}(n\log n)$, where

[TABLE]

It is straightforward to verify that the number of Chebyshev nodes within an interval of radius $1/n$ around the point $-1<x<1$ is $\mathcal{O}(1/\sqrt{1-x^{2}})$, that is,

[TABLE]

This estimate, together with the fact that the first and last Chebyshev node are distance at least $1/n^{2}$ from $1$ and $-1$, respectively, gives the estimate

[TABLE]

Let $\pi/2>\eta>0$ be a fixed parameter; direct calculation yields

[TABLE]

Combining this estimate with (30) yields $N_{\delta}=\mathcal{O}(n\log n)$ as was to be shown. ∎

5. Numerical results and implementation details

5.1. Numerical results

We report numerical results for two different point distributions: uniformly random points in $[1,10]$, and Chebyshev nodes in $[-1,1]$. In both cases, we choose the weights $\alpha_{1},\ldots,\alpha_{n}$ uniformly at random from $[0,1]$, and test the algorithm for

[TABLE]

We time two different versions of the algorithm: a standard implementation, and an implementation that uses precomputed exponentials. Precomputing exponentials may be advantageous in situations where the expression

[TABLE]

must be evaluated for many different weights $\alpha_{1},\ldots,\alpha_{n}$ associated with a fixed set of points $x_{1},\ldots,x_{n}$, see Remark 3.1. We find that using precomputed exponentials makes the algorithm approximately ten times faster, see Tables 1, 2, and 3. In addition to reporting timings, we report the absolute relative difference between the output of the algorithm of §3 and the output of direct evaluation; we define the absolute relative difference $\epsilon_{r}$ between the output $\tilde{u}_{j}$ of the algorithm of §3 and the output $u_{j}^{d}$ of direct calculation by

[TABLE]

Dividing by $\bar{u}_{j}$ accounts were the fact that the calculations are performed in finite precision; any remaining loss of accuracy in the numerical results is a consequence of the large number of addition and multiplication operations that are performed. All calculations are performed in double precision, and the algorithm of §3 is run with $\varepsilon=10^{-15}$. The parameter $\delta>0$ is set via an empirically determined heuristic. The numerical experiments were performed on a laptop with a Intel Core i5-8350U CPU and $7.7$ GiB of memory; the code was written in Fortran and compiled with gfortran with standard optimization flags. The results are reported in Tables 1, 2, and 3.

To put the run time of the algorithm in context, we additionally perform a time comparison to the Fast Fourier Transform (FFT), which also has complexity $\mathcal{O}(n\log n)$. Specifically, we compare the run time of the algorithm of §3 on random data using precomputed exponentials with the run time of an FFT implementation from FFTPACK [20] on random data of the same length using precomputed exponentials. We report these timings in Table 4; we find that the FFT is roughly $5$-$10$ times faster than our implementation of the algorithm of §3; we remark that no significant effort was made to optimize our implementation, and that it may be possible to improve the run time by vectorization.

5.2. Computing nodes and weights

The algorithm of §3 is described under the assumption that nodes $t_{1},\ldots,t_{m}$ and weights $w_{1},\ldots,w_{m}$ are given such that

[TABLE]

where $\varepsilon>0$ and $\delta>0$ are fixed parameters. As in the proof of Lemma 4.5 we note that by rescaling $r$ it suffices to find nodes and weights satisfying

[TABLE]

Indeed, if the nodes $t_{1},\ldots,t_{m}$ and weights $w_{1},\ldots,w_{m}$ satisfy (34), then the nodes $t_{1}/(b-a),\ldots,t_{m}/(b-a)$ and weights $w_{1}/(b-a),\ldots,w_{m}/(b-a)$ will satisfy (33). Thus, in order to implement the algorithm of §3 it suffices to tabulate nodes and weights that are valid for $r\in[1,M]$ for various values of $M$. In the implementation used in the numerical experiments in this paper, we tabulated nodes and weights valid for $r\in[1,M]$ for

[TABLE]

For example, in Tables 5 and 6 we have listed $m=33$ nodes $t_{1},\ldots,t_{33}$ and weights $w_{1},\ldots,w_{33}$ such that

[TABLE]

for all $r\in[1,1024]$.

The nodes and weights satisfying (34) can be computed by using a procedure for generating generalized Gaussian quadratures for Chebyshev systems together with the proof of Lemma 4.4. Indeed, Lemma 4.4 is constructive with the exception of the step that invokes Lemma 4.2 of Kreĭn. The procedure described in [4] is a constructive version of Lemma 4.2: given a Chebyshev system of functions, it generates the corresponding quadrature nodes and weights. We remark that generalized Gaussian quadrature generation codes are a powerful tools for numerical computation with a wide range of applications. The quadrature generation code used in this paper was an optimized version of [4] recently developed by Serkh for [19].

Acknowledgements

The authors would like to thank Jeremy Hoskins for many useful discussions. Certain commercial equipment is identified in this paper to foster understanding. Such identification does not imply recommendation or endorsement by the National Institute of Standards and Technology, nor does it imply that equipment identified is necessarily the best available for the purpose.