Legendre Series Analysis and Computation via Composed Abel-Fourier Transform

1. Introduction

The computation of the coefficients of Legendre expansions is a very important problem in numerical analysis and applied mathematics with a wide range of applications including, just to mention a few, approximation theory [1], special function theory [2], spectral methods for differential equations [3] and construction of quadrature formulae [3,4]. Its importance has emerged also in connection with the computation of spectra of highly oscillatory Fredholm integral operators, which play an important role in laser engineering [5]. For its relevance, this problem attracted a significant reserach attention since the Seventies [6,7].

The difficulty of the problem lies essentially in the fact that these coefficients are represented by integrals whose integrands oscillate rapidly for large values of the index of the polynomial. Standard quadrature procedures for the calculation of N Legendre coefficients lead only to slow$\mathcal{O}\left(N^2\right)$ algorithms (see, e.g., Ref. [6,8]). The first contribution for a more efficient computation of Legendre coefficients traces back to the work of Orszag [9], where the algorithm uses a slowly converging first-order WKB expansion of the Legendre polynomials.

More efficiently, in Ref. [10] the Legendre coefficients are obtained by transforming the corresponding Chebyshev coefficients through a multipole-like expansion, which yields a fast $\mathcal{O}(NlogN)$ algorithm, though requiring a considerable and rather expensive initialization phase. In this context, various improvements have been proposed, e.g., in [11,12]. Remaining within this kind of approach, [13] describes an $\mathcal{O}(N{(logN)}^2/loglogN)$ Chebyshev-Legendre transform, which is based on the Stieltjes’ asymptotic formula for the Legendre polynomials of large degree. Mori et al. [14] employed the same asymptotic formula to produce a fast $\mathcal{O}(NlogN)$ algorithm but affected by a problem of numerical instability for large N. The connection between Legendre and Chebyshev coefficients is analyzed also in [15] in the case of piecewise smooth functions. In Ref. [16] (see also [17]) an $\mathcal{O}\left(N{(logN)}^2\right)$ algorithm is given, for N a power of two, which requires a suitably preprocessed data structure.

Another way to tackle this problem has been described in [18], where an algorithm for the rapid computation of the Legendre coefficients is presented when the analytic expression of the input function is known on a Bernstein ellipse in the complex plane. The algorithm has $\mathcal{O}(NlogN)$ complexity but requires the knowledge of the region of analyticity of the function in $\mathbb{C}$.

In this paper we present an alternative procedure. The basic idea of our method consists in exploiting the Dirichlet-Murphy integral representation of the Legendre polynomials. We prove that the coefficients of the Legendre expansion of a function $f\left(x\right)$ are connected with a subset of the Fourier coefficients (the ones with nonnegative index) of an Abel-type transform of $f\left(x\right)$.

The numerical implementation of the algorithm follows straightforwardly and is very efficient. The aforementioned Fourier coefficients (which represent the sought for Legendre coefficients) can be computed in $\mathcal{O}(NlogN)$ operations by a single Fast Fourier Transform (FFT) after the evaluation of the Abel-type integral by means of standard quadrature techniques.

The dual problem of calculating the values of (finite) Legendre expansions is analyzed in Section 3. In SubSection 3.1 we prove that the Legendre expansions can be computed by means of a (inverse) cosine transform of a sequence of coefficients, which is obtained by a suitable linear transformation of the given set of N Legendre coefficients. Further, a novel algorithm is presented in SubSection 3.2, where the Legendre expansions are proved to be the Abel-type transform of a Fourier series, the coefficients of the latter being the Legendre coefficients of the function $f\left(x\right)$. This leads to an efficient algorithm which requires the computation of one FFT and of one Abel-type integral.

2. Connection between Legendre expansions and Fourier series

The standard form of the Legendre expansion reads:

$$f\left(x\right)=\sum _{n=0}^{\infty }{c}_n{P}_n\left(x\right)(x\in [-1,1]),$$

(1)

where $P_n\left(x\right)$ are the Legendre polynomials, which can be defined by the generating function [19]:

$$\sum _{n=0}^{\infty }{P}_n\left(x\right)t^n={\left(1-2xt+t^2\right)}^{-\frac{1}{2}}\left(\right|x|⩽1,\left|t\right|<1),$$

(2)

and the coefficients ${\left\{c_n\right\}}_{n=0}^{\infty }$ are given by:

$$c_n=\left(n+\frac{1}{2}\right){\int }_{-1}^{1}f\left(x\right)P_n\left(x\right)\mathrm{d}x(n⩾0).$$

(3)

The conditions to be satisfied by $f\left(x\right)$ to guarantee the uniform convergence of the series in (1) will be discussed in the next section. However, for our current purpose of computing the Legendre coefficients $c_n$, it is sufficient to assume that $f\left(x\right)$ be absolutely integrable on the interval $[-1,1]$. For later convenience, we define the normalized Legendre coefficients as

$$a_n\doteq \frac{c_n}{2n+1}(n⩾0).$$

(4)

We can now state the following theorem.

Theorem 1. 

The normalized Legendre coefficients ${\left\{a_n\right\}}_{n=0}^{\infty }$ coincide with the Fourier coefficients, with $n⩾0$, of an Abel-type transform of $f\left(x\right)$, that is:

$$a_n={\int }_{-\pi }^{\pi }\widehat{f}\left(y\right)e^{\mathrm{i}ny}\mathrm{d}y(n⩾0),$$

(5)

where the $2\pi$-periodic function $\widehat{f}\left(y\right)$ is defined by

$$\widehat{f}\left(y\right)=\frac{1}{2\pi \mathrm{i}}\epsilon \left(y\right)e^{\mathrm{i}\frac{y}{2}}{\int }_{cosy}^1\frac{f\left(x\right)}{{\left[2(x-cosy)\right]}^{\frac{1}{2}}}\mathrm{d}x(y\in \mathbb{R}),$$

(6)

$\epsilon \left(y\right)$ being the sign function.

Proof. 

We consider the Dirichlet-Murphy integral representation of the Legendre polynomials [20](Ch. III, §5.4):

$$P_n(cosu)=-\frac{\mathrm{i}}{\pi }{\int }_u^{(2\pi -u)}\frac{e^{\mathrm{i}(n+\frac{1}{2})y}}{{\left[2(cosu-cosy)\right]}^{\frac{1}{2}}}\mathrm{d}y.$$

(7)

Plugging (7) into equality (3) (after the change of variable $x\to cosu$), we have:

$$2\pi \mathrm{i}{a}_n={\int }_0^{\pi }\mathrm{d}uf(cosu)sinu{\int }_u^{(2\pi -u)}\frac{e^{\mathrm{i}(n+\frac{1}{2})y}}{{\left[2(cosu-cosy)\right]}^{\frac{1}{2}}}\mathrm{d}y.$$

(8)

Interchanging the order of integration in (8) we obtain:

$$\begin{array}{cc}& 2\pi \mathrm{i}{a}_n={\int }_0^{\pi }\mathrm{d}y{e}^{\mathrm{i}(n+\frac{1}{2})y}{\int }_0^{y}f(cosu)\frac{sinu}{{\left[2(cosu-cosy)\right]}^{\frac{1}{2}}}\mathrm{d}u\hfill \\ & +{\int }_{\pi }^{2\pi }\mathrm{d}y{e}^{\mathrm{i}(n+\frac{1}{2})y}{\int }_0^{(2\pi -y)}f(cosu)\frac{sinu}{{\left[2(cosu-cosy)\right]}^{\frac{1}{2}}}\mathrm{d}u.\hfill \end{array}$$

(9)

Next, if we make the change of variables: $y\to y-2\pi$ and $u\to -u$, the second integral on the r.h.s. of (9) becomes:

$$e^{\mathrm{i}\pi }{\int }_{-\pi }^0\mathrm{d}y{e}^{\mathrm{i}(n+\frac{1}{2})y}{\int }_0^{y}f(cosu)\frac{sinu}{{\left[2(cosu-cosy)\right]}^{\frac{1}{2}}}\mathrm{d}u.$$

(10)

Then, (9) becomes

$$\begin{array}{cc}& 2\pi \mathrm{i}{a}_n={\int }_0^{\pi }\mathrm{d}y{e}^{\mathrm{i}(n+\frac{1}{2})y}{\int }_0^{y}f(cosu)\frac{sinu}{{\left[2(cosu-cosy)\right]}^{\frac{1}{2}}}\mathrm{d}u\hfill \\ & +e^{\mathrm{i}\pi }{\int }_{-\pi }^0\mathrm{d}y{e}^{\mathrm{i}(n+\frac{1}{2})y}{\int }_0^{y}f(cosu)\frac{sinu}{{\left[2(cosu-cosy)\right]}^{\frac{1}{2}}}\mathrm{d}u,\hfill \end{array}$$

(11)

which, after the change of variable $cosu\to x$ into the integrals on the r.h.s., yields:

$$a_n={\int }_{-\pi }^{\pi }\widehat{f}\left(y\right)e^{\mathrm{i}ny}\mathrm{d}y(n⩾0),$$

(12)

with $\widehat{f}\left(y\right)$ given by (6). □

It is easy to check from (6) that $\widehat{f}\left(y\right)$ satisfies the following symmetry relation:

$$\widehat{f}\left(y\right)=-e^{\mathrm{i}y}\widehat{f}(-y),$$

(13)

which, in view of (5), induces the following symmetry on the coefficients $a_n$:

$$a_n=-a_{-n-1}(n\in \mathbb{Z}).$$

(14)

2.1. Numerical issues

Let us now consider the actual problem of computing the first N Legendre coefficients ${\left\{c_n\right\}}_{n=0}^{N-1}$ (or, equivalently, Nnormalized coefficients ${\left\{a_n\right\}}_{n=0}^{N-1}$) associated with the function $f\left(x\right)$ (see (1) and (3)).

The numerical implementation of the algorithm suggested by Theorem 1 requires the computation of the Fourier coefficients (5) of the Abel-type integral function $\widehat{f}\left(y\right)$ defined in (6). For what concerns the computation of $\widehat{f}\left(y\right)$, the integrand presents a weak algebraic singularity at the end point of the domain of integration, which can be effectively handled by means of a proper nonlinear change of variable. This technique, along with the use of a standard quadrature formula (e.g., the Gauss-Legendre one), allows obtaining high accuracy with a small number of nodes [21]. Quadrature formulae suited to end-point singular integrands can also be used [22,23].

It is ought to stress that the function $\widehat{f}\left(y\right)$ does not depend on the order n of the Legendre polynomial. Therefore, the computational cost for evaluating $\widehat{f}\left(y\right)$ is independent of the number N of Legendre coefficients to be computed. This means that the number of knots which are necessary to compute the integral in $\widehat{f}\left(y\right)$ (with a prescribed fixed accuracy $ϵ$) does not depend on N, but depends only upon the oscillatory characteristics of the function $f\left(x\right)$ (i.e., the number of samples needed to determine $f\left(x\right)$ uniquely) and on the accuracy $ϵ$. The oscillatory contributions due to the Legendre polynomials or, in other words, the dependence on the order n of the Legendre polynomial, are accounted for only in the Fourier formula (5) and are neatly separated from the oscillatory contributions ascribable to the input function $f\left(x\right)$. The consequence of this decoupling is the low computational complexity of the whole algorithm. In fact, formula (5) makes it possible to take full advantage of the computational efficiency of the Fast Fourier Transform both in terms of speed of computation and of accuracy [24,25]. The first N Fourier coefficients $a_n$ of $\widehat{f}\left(y\right)$ (which coincide with the normalized Legendre coefficients of $f\left(x\right)$) can thus be computed by a single FFT from the values (samples) of $\widehat{f}\left(y\right)$ at N distinct points in $\mathcal{O}(NlogN)$ time. As we have seen above, the number of operations for the calculation of each sample of $\widehat{f}\left(y\right)$ is independent of N and, consequently, the (asymptotic) computational complexity of the entire algorithm for the computation of the Legendre coefficients coincides with that of the FFT, i.e., $\mathcal{O}(NlogN)$. Moreover, in view of the symmetry relation (13), the number of calls to the procedure for the computation of the Abel integral is halved. Of course, if $\mathcal{O}\left(N\right)$ samples are needed to determine uniquely the function $f\left(x\right)$ and we want the error $ϵ$ to be as small as possible, then all the $\mathcal{O}\left(N\right)$ samples of $f\left(x\right)$ must be used for the numerical computation of the Abel integral and, consequently, the total number of operations in the algorithm will be $\mathcal{O}\left(N^2\right)$. Therefore, following logical lines similar to those adopted by Alpert and Rokhlin in Ref. [1] (and, more generally, in methods like, e.g., the Fast Multipole Methods [26]), the reduction of complexity from $\mathcal{O}\left(N^2\right)$ to $\mathcal{O}(NlogN)$ is obtained by accepting a lower precision in the computation (in our case, a lower precision in the computation of the Abel integral).

We want to remark two additional features of the algorithm just presented. First, the algorithm does not require the knowledge of the input function $f\left(x\right)$ at a specific and prescribed set of knots (see, for instance, Ref. [10] where the input function is supposed to be known at Chebyshev knots). The function $f\left(x\right)$ enters the algorithm only through its Abel-type integral. This allows a significant flexibility in choosing the quadrature scheme which is more suitable for the distribution of knots at which the input function is known (e.g., uniform grid, Chebyshev grid, set of measurements of $f\left(x\right)$ non-uniformly distributed).

The final remark concerns the robustness of the algorithm against the noise. In the case the values of $f\left(x\right)$ are only approximately known, e.g., when the input samples of $f\left(x\right)$ represent measurements affected by error, the effects of the noise in the calculation of the coefficients $c_n$ are damped by the algorithm. In fact, the noisy input samples of $f\left(x\right)$ are smoothed by the Abel integral operator [27] before being analyzed (encoded as samples of $\widehat{f}\left(y\right)$) by the oscillatory Fourier integral (5).

The procedure described above has been implemented in double precision arithmetics using the standard open source GNU Scientific Library (GSL) [28]. Feasibility and accuracy of the algorithm have been verified by direct comparison of the obtained numerical results with the true Legendre coefficients for a variety of functions. Hereafter we report some of them.

Table 1 displays the absolute error $E_n\doteq |c_n^{\mathrm{True}}-c_n^{\mathrm{Computed}}|$ in the computation of Legendre coefficients for four functions with different degree of regularity: the discontinuous function $f_1\left(x\right)=\gamma \epsilon (x-x_0)$ (with $\epsilon (·)$ being the sign function and $\gamma$ a constant), whose Legendre coefficients are $c_0^{\left(1\right)}=-\gamma x_0$, $c_n^{\left(1\right)}=\gamma \left[P_{n-1}\left(x_0\right)-P_{n+1}\left(x_0\right)\right]$ for $n⩾1$; the $C^0[-1,1]$ function: $f_2\left(x\right)={\left|x\right|}^{3/2}$, whose Legendre coefficients are [3](p. 78):

$$c_n^{\left(2\right)}=\left\{\begin{array}{cc}{(\alpha +1)}^{-1}\hfill & \mathrm{if}\mathrm{n}=0(\alpha =\frac{3}{2}),\hfill \\ \frac{{\displaystyle \phantom{(}(2n+1)\alpha (\alpha -2)\cdots (\alpha -n+2)}}{{\displaystyle \phantom{(}(\alpha +1)(\alpha +3)\cdots (\alpha +n+1)}}\hfill & \mathrm{if}\mathrm{n}\mathrm{is}\mathrm{even},\hfill \\ 0\hfill & \mathrm{if}\mathrm{n}\mathrm{is}\mathrm{odd};\hfill \end{array}\right.$$

(15)

the function $f_3\left(x\right)={(1-2xt+t^2)}^{-1/2}$ with $c_n^{\left(3\right)}=t^n$ (with a singlarity for x outside the interval $[-1:1]$, see (2)); the $\mathbb{C}$-analytic function: $f_4\left(x\right)=e^{\beta x}{J}_0\left(\beta \sqrt{1-x^2}\right)$ with $c_n^{\left(4\right)}={\beta }^n/n!$ ($J_0$ denoting the Bessel function of the first kind and zeroth order).

The first observation from Table 1 is that the accuracy depends on the degree of smoothness of the function because of the computation of the Abel transform of $f\left(x\right)$ that, in order to obtain a given target precision, requires more quadrature knots for low regular functions. For what concerns the smooth (within the interval $[-1,1]$) functions $f_3$ and $f_4$, which give similar results, the accuracy increases significantly to reach values comparable with the $ϵ$-machine.

The algorithm sources of error are: the approximation of integral (5) by the Discrete Fourier Transform; the approximation of the Abel transform (6) by quadratures and, finally, the finite accuracy of the floating-point arithmetic. The first source of error is ruled by N. This is clearly visible from Figure 1(A), where the maximum error $E_{max}\left(N\right)\doteq {max}_{n\in [0,N-1]}{E}_n$ is plotted versus N. For the $C^0$-function $f_2\left(x\right)$ (asterisks), we see that for low values of N the maximum error is mainly ascribable to the error due to the approximation of the Fourier integral with the DFT and to the approximation of the Abel integral: it decreases like $N^{-2}$ up to $N\sim {10}^3$, and then starts increasing like $\mathcal{O}\left(N\right)$, as in the case of smooth functions. From the comparison with the plot of maximum error for smooth functions (triangles), we can also say that the number of points N used in the DFT need not to be unnecessary large unless the function $f\left(x\right)$ itself has low regularity.

The computation of N Legendre coefficients $c_n$ by ordinary quadrature produces an $\mathcal{O}\left(N^2\right)$ algorithm. The increment of performances of the algorithm presented here has been verified by evaluating the speed of computation at (nearly) same precision. The results of these tests are illustrated in Figure 1(B), where the execution time T (in seconds) is plotted versus the number N of computed Legendre coefficients. The results show an $\mathcal{O}\left(N\right)$ execution time (at least in the N-interval we considered, up to $N\sim 5·{10}^5$) independently of the regularity properties of the input function $f\left(x\right)$, confirming the expected great increase of computational speed (asymptotically, the awaited improvement ratio being proportional to $N/logN$). Such an increase of performances will become even more crucial for the efficient evaluation of multivariate Legendre transform [5] and of spherical harmonic expansions [30], which will be the subject of a forthcoming paper.

3. Computation of Legendre expansions

Let us now move on to consider the problem of evaluating the function $f\left(x\right)$ from the set of Legendre coefficients $c_n$. In general, in the complex plane the domain of convergence of the Legendre expansion (1) is the maximal ellipse with foci at $\pm 1$ within which the function f is analytic [31]. This rather restricted class of functions can be greatly enlarged when we limit ourselves to consider convergence issues on the segment $x\in [-1,1]$. In Ref. [32] W. H. Young showed the tight connection between Legendre series and Fourier series when $x\in (-1,1)$. In particular, he proved that in any internal closed interval of $(-1,1)$ the series of Legendre (1) behaves in respect to convergence, divergence, or oscillation, uniform or otherwise, precisely like the Fourier series of $f(cosu)$ where $x=cosu$, provided the terms $c_n$ of series (1) are such that

$$n^{-\frac{1}{2}}{c}_n\to 0\mathrm{for}n\to \infty ,$$

(16)

the latter being the natural necessary condition for the convergence of the series (1). More explicitly, e.g., E. W. Hobson [33] provided a test of convergence stating that if the condition

$${\int }_{-1}^1{(1-x^2)}^{-\frac{1}{4}}\left|f\left(x\right)\right|\mathrm{d}x<\infty$$

(17)

is satisfied (the latter condition being equivalent to require that $f(cosu)\sqrt{sinu}$ be integrable on $u\in [0,\pi ]$), then the series (1) is convergent for any $x\in (-1,1)$ near which f is of bounded variation. Since the points $x=\pm 1$ are singular points of Legendre’s equation, the convergence at the end points of the interval requires somehow more stringent (than in the Fourier case) conditions such as, e.g., $f\left(x\right)$ being of bounded variation over the entire interval $[-1,1]$.

This parallelism between Legendre and Fourier expansions have been already exploited in the previous section where we proved that the (normalized) Legendre coefficients of $f\left(x\right)$ coincide with the Fourier coefficients of the function $\widehat{f}\left(y\right)$. In the next two subsections we continue our analysis along this line, and present two algorithms for the computation of values of Legendre expansions.

3.1. An algorithm for the evaluation of Legendre series

In the theorem that follows we will prove that a Legendre expansion can be computed by calculating the (inverse) cosine transform of a sequence which is obtained by a suitable linear transformation of the set of Legendre coefficients $c_n$.

Theorem 2. 

Let $f\left(x\right)$ be the function given in (1). Assume that $f\left(x\right)$ is such that the Legendre expansion (1) and the Fourier expansion of $f(cosu)$ converge uniformly on $x\in [-1,1]$ and $u\in [0,\pi ]$, respectively. Then $f(cosu)$, $u\in [0,\pi ]$, can be written as follows:

$$f(cosu)=\frac{{\phi }_0}{2}+\sum _{m=1}^{+\infty }{\phi }_{m}cosmu,$$

(18)

the coefficients ${\phi }_m$ ($m=0,1,2,\dots$) being given by:

$${\phi }_m=\sum _{n=0}^{+\infty }{k}_{m,n}{c}_n(m=0,1,2,\dots ),$$

(19)

$c_n$ denoting the Legendre coefficients associated with $f\left(x\right)$ (see (3)) and

$$k_{m,n}=\frac{2}{\pi }{\int }_0^{\pi }{P}_n(cosu)cos\left(mu\right)\mathrm{d}u(m,n=0,1,2,\dots ).$$

(20)

Proof. 

Since $f\left(x\right)$ is compactly supported on $x\in [-1,1]$, the function $f(cosu)$ ($u\in [0,\pi ]$) is a $2\pi$-periodic function of u. Then, the latter can be represented as the Fourier series,

$$f(cosu)=\frac{1}{2\pi }\sum _{n=-\infty }^{+\infty }{f}_n{e}^{-\mathrm{i}nu}(u\in [0,\pi ]).$$

(21)

In view of the uniform convergence of the two series, from (1) and (21) we have:

$$\frac{1}{2\pi }\sum _{n=-\infty }^{+\infty }{f}_n{e}^{-\mathrm{i}nu}=\sum _{n=0}^{+\infty }{c}_n{P}_n(cosu),$$

(22)

so that, multiplying both sides of (22) by $e^{\mathrm{i}mu}$ and integrating on $u\in [-\pi ,\pi ]$, we have:

$$\sum _{n=-\infty }^{+\infty }{f}_n\left[\frac{1}{2\pi }{\int }_{-\pi }^{\pi }{e}^{\mathrm{i}(m-n)u}\mathrm{d}u\right]=\sum _{n=0}^{+\infty }{c}_n{\int }_{-\pi }^{\pi }{P}_n(cosu)e^{\mathrm{i}mu}\mathrm{d}u.$$

(23)

The term in the parenthesis on the l.h.s. of (23) is the integral representation of the Kronecker delta ${\delta }_{m,n}$. Then, from (23) we have:

$${\phi }_m\doteq \frac{f_m}{\pi }=\sum _{n=0}^{+\infty }\left[\frac{2}{\pi }{\int }_0^{\pi }{P}_n(cosu)cos\left(mu\right)\mathrm{d}u\right]·c_n=\sum _{n=0}^{+\infty }{k}_{m,n}{c}_n,$$

(24)

with $k_{m,n}$ defined in (20). Finally, expansion (18) follows readily by observing that ${\phi }_{-m}={\phi }_m$. □

It is worth remarking that the transformation $K\doteq {\left\{k_{m,n}\right\}}_{m,n=0}^{\infty }$, which connects the Legendre coefficients of $f\left(x\right)$ with the Fourier coefficients of $f(cosu)$, is independent of $f\left(x\right)$, can be computed with arbitrary precision, and can be regarded as known. In the next proposition we list some known facts on the structure of the matrix K, we give the analytic expression of its entries as smooth function of the indexes (see Ref. [10]), and provide an iterative procedure for computing its elements.

Proposition 1.(i) The coefficients $k_{m,n}$ (defined in (20)) can be written as follows:

$$k_{m,n}=\frac{2}{\pi }{\int }_{-1}^1{P}_n\left(x\right)T_n\left(x\right)\frac{\mathrm{d}x}{\sqrt{1-x^2}},$$

(25)

$T_n(·)$ denoting the Chebyshev polynomials of the first kind.

(ii) The matrix K is upper triangular, and its entries $k_{m,n}$ can be computed as follows:

$$k_{m,n}=\left\{\begin{array}{cc}\frac{1}{\pi }{\left[\Lambda \left(\frac{n}{2}\right)\right]}^2\hfill & \mathrm{for}m=0\mathrm{and}n⩾m\mathrm{with}n\mathrm{even},\hfill \\ \frac{2}{\pi }\Lambda \left(\frac{n-m}{2}\right)\Lambda \left(\frac{n+m}{2}\right)\hfill & \mathrm{if}0<m⩽n\mathrm{and}(n+m)\mathrm{is}\mathrm{even},\hfill \\ 0\hfill & \mathrm{otherwise},\hfill \end{array}\right.$$

(26)

where

$$\Lambda \left(z\right)\doteq \frac{\Gamma (z+\frac{1}{2})}{\Gamma (z+1)},$$

(27)

$\Gamma (·)$ denoting the Euler gamma function.

(iii) The following recurrence relation holds for $n=1,2,\dots$ and $0⩽m<n$:

$$k_{m,n}=\frac{(n-1)(n-m-1)}{n(m+n)}{k}_{m,n-2}+\frac{m(2n-1)}{n(m+n)}{k}_{m-1,n-1}.$$

(28)

(iv) The diagonal elements can be written as follows:

$$k_{m,m}=2\prod _{j=1}^m\left(1-\frac{1}{2j}\right)(m=1,2,\dots ).$$

(29)

(v) The elements $k_{0,2n}$ of the first row can be written as follows:

$$k_{0,2n}=2\prod _{j=1}^n{\left(1-\frac{1}{2j}\right)}^2(n=1,2,\dots ).$$

(30)

Proof. 

Formula (25) follows from definition (20) by making the substitution $cosu\to x$ into the integral and recalling the definition of the Chebyshev polynomials of the first kind: $T_m\left(x\right)\doteq cos(marccosx)$, ($m⩾0$). Representation (ii) of the coefficients $k_{m,n}$ is well-known, and it is given, e.g., by Alpert and Rokhlin in [10](formula (20)) (where the matrix K is there named M). The recurrence relation (28) follows from representation (26) and from the following three-term recurrence relations holding, respectively, for the Legendre polynomials $P_n\left(x\right)$ and for Chebyshev polynomials $T_n\left(x\right)$:

$$\begin{array}{cc}\hfill & n{P}_n\left(x\right)=(2n-1)x{P}_{n-1}\left(x\right)-(n-1)P_{n-2}\left(x\right)(n=2,3,\dots ),\hfill \end{array}$$

(31)

$$\begin{array}{cc}\hfill & P_0\left(x\right)=1\mathrm{and}{P}_1\left(x\right)=x,\hfill \end{array}$$

(32)

and

$$\begin{array}{cc}\hfill 1-1& T_n\left(x\right)=2x{T}_{n-1}\left(x\right)-T_{n-2}\left(x\right)(n=2,3,\dots ),\hfill \end{array}$$

(33)

$$\begin{array}{cc}\hfill & T_0\left(x\right)=1\mathrm{and}{T}_1\left(x\right)=x.\hfill \end{array}$$

(34)

Finally, point (iv) is obtained by putting $n=m$ into (29), and, similarly, point (v) follows by putting $m=0$ into the relation (29). □

The algorithmic prescription provided by Theorem 2 for the computation of the values of a Legendre expansion is indeed very simple. It consists of a two-step procedure:

• From the N-vector $\mathbf{c}={(c_0,c_1,\dots ,c_{N-1})}^T$ of Legendre coefficients compute the N-vector $\phi ={({\phi }_0,{\phi }_1,\dots ,{\phi }_{N-1})}^T$ by $\phi =K_N\mathbf{c}$, where the (known) upper triangular matrix $K_N={\left\{k_{m,n}\right\}}_{m,n=0}^{N-1}$ is defined in (20) (see also Proposition 1);
• Compute the cosine transform of length N of the vector $\phi$ to obtain values of $f(cosu)$ at N selected points of the interval $u\in [0,\pi ]$.

The second step of the algorithm presents no difficulties. It is fast since it can be implemented by a single FFT of length N, which requires $\mathcal{O}(NlogN)$ operations. The standard FFT procedure provides the output on a uniform grid ${\left\{u_n\right\}}_{n=0}^{N-1}$ of $[0,\pi ]$, so that, according to formula (18), the final result consists of the set of function values ${\left\{f(cos{u}_n)\right\}}_{n=0}^{N-1}$ at Chebyshev points.

The first step of the algorithm presents instead some aspects which deserve a few comments. The crux is that the matrix $K_N$ is dense and the direct calculation of the matrix-vector product $K_N\mathbf{c}$ requires $\mathcal{O}\left(N^2\right)$ operations (precisely, $\frac{N(N+2)}{4}$ multiplications). Therefore, this step can become a bottleneck for the algorithm when N becomes large. In some practical cases the question is not critical, e.g., when the triangular structure of $K_N$ can be exploited on some computer architectures (e.g., vector machines, parallel GPU) to drop significantly the complexity of the computation. Nevertheless, the reduction of the algorithmic complexity remains an issue. Approximate methods are the only tools to handle effectively the fast computation of the matrix-vector product when the matrix is not structured (i.e., when its entries depend only on $\mathcal{O}\left(N\right)$ parameters). The key to many of these fast methods is renouncing to exactness and accepting the result of the computation within an a-priori fixed level of accuracy. This allows for using approximations of the kernel K which, combined with a suitable subdivision of the matrix into panels, can reduce enormously the cost of computing. Among these algorithms it is worth citing the celebrated Fast Multipole Methods by L. Greengard and V. Rokhlin for the fast evaluation of N-body interactions [12,26,34], which are able to reduce the number of operations to $\mathcal{O}(NlogN)$ (or even $\mathcal{O}\left(N\right)$, depending on the matrix structure).

The specific problem we are dealing with here, that is, the fast computation of the matrix-vector product $K_N\mathbf{c}$ (with ${\left\{k_{m,n}\right\}}_{m,n=0}^{N-1}$ given by formula (25)) has been solved brilliantly by B. Alpert and V. Rokhlin in Ref. [10]. In that paper the matrix $K_N$ is first properly divided in square submatrices, and then the computation associated with each submatrix is performed efficiently approximating each entry of the submatrix by its finite Chebyshev expansion. The computational cost brought by a $\nu \times \nu$ submatrix is then reduced from order $\mathcal{O}\left({\nu }^2\right)$ of the naive computation to $\mathcal{O}(\nu log\frac{1}{ϵ})$, $ϵ$ being the fixed precision required in the Chebyshev interpolation of the submatrix entries. This basic building block allows thus constructing an algorithm for the fast computation of the matrix-vector product $K_N\mathbf{c}$ with $\mathcal{O}(NlogN)$ order of complexity.

3.2. A novel algorithm

In this subsection we present a novel algorithm for the inversion of the Legendre transform, which somewhat represents the natural dual algorithm of the one presented in Section 2 for the computation of the forward Legendre transform. In the next theorem we will show that a Legendre expansion can be represented as the Abel-type transform of a suitable Fourier series, the coefficients of the latter being the Legendre coefficients of $f\left(x\right)$ (cf. Theorem 1).

Theorem 3. 

Let $f\left(x\right)$ be the function given in (1). Assume that $f\left(x\right)$ is such that the Legendre expansion (1) converges uniformly on $x\in [-1,1]$. Then $f(cosu)$, $u\in [0,\pi ]$, can be written as follows:

$$f(cosu)={\int }_u^{\pi }\frac{\chi \left(t\right)}{{\left[2(cosu-cost)\right]}^{\frac{1}{2}}}\mathrm{d}t(u\in [0,\pi ]),$$

(35)

where the function $\chi \left(t\right)$ is defined by

$$\chi \left(t\right)\doteq \frac{2}{\pi }{\sum }_{n=0}^{\infty }{c}_{n}sin(n+\frac{1}{2})t(t\in [0,\pi ]),$$

(36)

$c_n$ denoting the Legendre coefficients associated with $f\left(x\right)$ (see (3)).

Proof. 

First we show that the following integral representation holds for the Legendre polynomials $P_n(cosu)$:

$$P_n(cosu)=\frac{2}{\pi }{\int }_u^{\pi }\frac{sin(n+\frac{1}{2})t}{{\left[2(cosu-cost)\right]}^{\frac{1}{2}}}\mathrm{d}t(u\in [0,\pi ]).$$

(37)

Formula (37) can be easily obtained from the Dirichlet-Murphy representation (7) by first splitting the interval of integration into two subintervals: $[u,2\pi -u]=[u,\pi ]\cup [\pi ,2\pi -u]$, and then changing the integration variable $t\to 2\pi -t$ in the second integral, i.e., we have:

$$\begin{array}{cc}& P_n(cosu)=-\frac{\mathrm{i}}{\pi }\left({\int }_u^{\pi }+{\int }_{\pi }^{2\pi -u}\right)\frac{e^{\mathrm{i}(n+\frac{1}{2})t}}{{\left[2(cosu-cost)\right]}^{\frac{1}{2}}}\mathrm{d}t\hfill \\ & =-\frac{\mathrm{i}}{\pi }\left\{{\int }_u^{\pi }\frac{e^{\mathrm{i}(n+\frac{1}{2})t}}{{\left[2(cosu-cost)\right]}^{\frac{1}{2}}}\mathrm{d}t-{\int }_u^{\pi }\frac{e^{-\mathrm{i}(n+\frac{1}{2})t}}{{\left[2(cosu-cost)\right]}^{\frac{1}{2}}}\mathrm{d}t\right\}\hfill \\ & =\frac{2}{\pi }{\int }_u^{\pi }\frac{sin(n+\frac{1}{2})t}{{\left[2(cosu-cost)\right]}^{\frac{1}{2}}}\mathrm{d}t.\hfill \end{array}$$

(38)

Now, we can plug formula (37) into Legendre expansion (1) and then, interchanging sum and integral (which is legitimate in view of the assumption of uniform convergence for expansion (1)), we get:

$$f(cosu)=\frac{2}{\pi }\sum _{n=0}^{\infty }{c}_n{\int }_u^{\pi }\frac{sin(n+\frac{1}{2})t}{{\left[2(cosu-cost)\right]}^{\frac{1}{2}}}\mathrm{d}t={\int }_u^{\pi }\frac{\left[\frac{2}{\pi }{\sum }_{n=0}^{\infty }{c}_{n}sin(n+\frac{1}{2})t\right]}{{\left[2(cosu-cost)\right]}^{\frac{1}{2}}}\mathrm{d}t,$$

(39)

which is indeed formula (35) with $\chi \left(t\right)$ defined by (36). □

Theorem 3 suggests the following algorithm for the computation of the values of a Legendre expansion:

• From the input N-vector $\mathbf{c}={(c_0,c_1,\dots ,c_{N-1})}^T$ of Legendre coefficients compute the N-term approximation ${\chi }_N\left(t\right)$ of the function $\chi \left(t\right)$ (see (36)):

  $${\chi }_N\left(t\right)\doteq \frac{2}{\pi }{\sum }_{n=0}^{N-1}{c}_{n}sin(n+\frac{1}{2})t.$$

  (40)
• Compute the values of the N-term Legendre expansion

  $$f_N(cosu)=\sum _{n=0}^{N-1}{c}_n{P}_n(cosu)(u\in [0,\pi ]),$$

  (41)

  at points $u_n$ ($n=0,1,\dots$) by the Abel-type integral (see (35)):

  $$f_N(cos{u}_n)={\int }_{u_n}^{\pi }\frac{{\chi }_N\left(t\right)}{{\left[2(cos{u}_n-cost)\right]}^{\frac{1}{2}}}\mathrm{d}t(u_n\in [0,\pi ]).$$

  (42)

The first step is fast: the set of values ${\left\{{\chi }_N\left(t_n\right)\right\}}_{n=0}^{N-1}$ of the function ${\chi }_N\left(t\right)$ at N equidistant nodes of the interval $[0,\pi ]$ can be calculated in $\mathcal{O}(NlogN)$ time by a single FFT.

For what concerns the second step, it should first be noted that the distribution of the nodes $u_n$ is not constrained and, therefore, we are free to calculate $f_N(cosu)$ at any set of points of the interval $u\in [0,\pi ]$ by means of formula (42).

The numerical computation of the Abel-type integral in (42), which presents a weak singularity at the left end-point (see also Section 2), is constrained by the fact that the integrand is known only on a uniform grid of the interval $t\in [0,\pi ]$ (see Step (1) of the algorithm) and, consequently, only quadrature formulae based on equally spaced knots can be used.

In general, the computation of $f_N(cos{u}_n)$ at N points ${\left\{u_n\right\}}_{n=0}^{N-1}$ requires $\mathcal{O}\left(N^2\right)$ operations since the second step of the algorithm requires $\mathcal{O}\left(N\right)$ operations when all the N known values of the function ${\chi }_N$ are used for the computation of the integral in (42). However, it is ought to note that for this computation it is not always necessary to implement a full N-knot quadrature procedure, in particular when N is large. The integral in (42) can usually be computed with an a-priori fixed precision by using much less (than N) quadrature knots, the number of knots depending on the implemented quadrature scheme and on the oscillatory characteristics of the integrand, which, ultimately, depend on the oscillatory characteristics of $f\left(x\right)$. An instance of a suitable quadrature scheme is given in Refs. [22,23], where quadratures rules are proposed, based on corrections to the trapezoidal rule, for integrands with end-point singularities of various types, including singularities of the form $x^{\gamma }$, $\left|\gamma \right|<1$. The convergence rate of these quadrature formulae is proved to be at least $\kappa$ (that is, the approximation error is of the order of $N^{-\kappa }$) where $\kappa$ denotes the order of the corrected quadrature rule [22,23]. Thus, in these cases and for N sufficiently large, the algorithm we propose for the evaluation of the Legendre expansion becomes convenient in terms of complexity since the number of operations of the second step is a constant independent of N and, for large values of N, the complexity of the entire algorithm will be dominated by the $\mathcal{O}(NlogN)$ complexity of the FFT used in Step (1).

4. Conclusions

In summary, we have presented a new and fast, i.e., $\mathcal{O}(NlogN)$, algorithm for the computation of the coefficients of Legendre expansions. The algorithm is very simple to implement and computes the Legendre coefficients with just a single FFT of an Abel-type transform of the input function. This algorithm calls up a number of natural generalizations to other polynomials systems, for instance, associated Legendre polynomials and the related spherical harmonic expansions.

The dual problem of evaluating a Legendre expansion has been also considered, and two algorithms have been proposed. In the first algorithm the values of the expansion at appropriately chosen points of the interval $x\in [-1,1]$ are computed by simply performing a single (inverse) FFT of a sequence of coefficients, which are obtained by a suitable linear transformation of the given set of N Legendre coefficients. The rapid evaluation of the expansion values requires using approximate methods for the fast calculation of the matrix-vector product. The second algorithm for the calculation of values of Legendre expansions is novel and somehow inverts the logical steps of the algorithm we have proposed for the computation of the Legendre coefficients: the N-term Legendre expansion can be evaluated efficiently (for N sufficiently large) at any point of the interval $x\in [-1,1]$ by computing the Abel-type transform of a suitable Fourier series.