NAG Toolbox Chapter Introduction

C06 — Summation of Series

Scope of the Chapter

Calculating the discrete Fourier transform of a sequence of real or complex data values.

(b)

Calculating the discrete convolution or the discrete correlation of two sequences of real or complex data values using discrete Fourier transforms.

(c)

Calculating the inverse Laplace transform of a user-supplied function.

(d)

Direct summation of orthogonal series.

(e)

Acceleration of convergence of a sequence of real values.

Background to the Problems

Discrete Fourier Transforms

Complex transforms

Most of the functions in this chapter calculate the finite discrete Fourier transform (DFT) of a sequence of nn complex numbers
zjzj, for j = 0,1, … ,n − 1j=0,1,…,n-1. The direct transform is defined by

n − 1

ẑk = 1/(sqrt(n))

∑

zjexp( − i(2πjk)/n)

j = 0

z^k=1n∑j=0n-1zjexp(-i2πjkn)

(1)

for k = 0,1, … ,n − 1k=0,1,…,n-1. Note that equation (1) makes sense for all integral kk and with this extension ẑkz^k is periodic with period nn, i.e., ẑk = ẑk ± nz^k=z^k±n, and in particular ẑ − k = ẑn − kz^-k=z^n-k. Note also that the scale-factor of 1/(sqrt(n))1n may be omitted in the definition of the DFT, and replaced by 1/n1n in the definition of the inverse.

If we write zj = xj + iyjzj=xj+iyj and ẑk = ak + ibkz^k=ak+ibk, then the definition of ẑkz^k may be written in terms of sines and cosines as

n − 1

ak = 1/(sqrt(n))

∑

(xjcos((2πjk)/n) + yjsin((2πjk)/n))

j = 0

ak=1n∑j=0n-1(xjcos(2πjkn)+yjsin(2πjkn))

n − 1

bk = 1/(sqrt(n))

∑

(yjcos((2πjk)/n) − xjsin((2πjk)/n)).

j = 0

bk=1n∑j=0n-1(yjcos(2πjkn)-xjsin(2πjkn)).

The original data values zjzj may conversely be recovered from the transform ẑkz^k by an inverse discrete Fourier transform:

n − 1

zj = 1/(sqrt(n))

∑

ẑkexp( + i(2πjk)/n)

k = 0

zj=1n∑k=0n-1z^kexp(+i2πjkn)

(2)

for j = 0,1, … ,n − 1j=0,1,…,n-1. If we take the complex conjugate of (2), we find that the sequence zjz-j is the DFT of the sequence ẑkz^-k. Hence the inverse DFT of the sequence ẑkz^k may be obtained by taking the complex conjugates of the ẑkz^k; performing a DFT, and taking the complex conjugates of the result. (Note that the terms forward transform and backward transform are also used to mean the direct and inverse transforms respectively.)

The definition (1) of a one-dimensional transform can easily be extended to multidimensional transforms. For example, in two dimensions we have

n1 − 1

n2 − 1

ẑk1k2 = 1/(sqrt(
n1n2
))

∑

∑

zj1j2exp( − i(2πj1k1)/(n1))exp( − i(2πj2k2)/(n2)).

j1 = 0

j2 = 0

z^k1k2=1n1n2∑j1=0n1-1∑j2=0n2-1zj1j2exp(-i2πj1k1n1)exp(-i2πj2k2n2).

(3)

Note: definitions of the discrete Fourier transform vary. Sometimes (2) is used as the definition of the DFT, and (1) as the definition of the inverse.

Real transforms

If the original sequence is purely real valued, i.e., zj = xjzj=xj, then

n − 1

ẑk = ak + ibk = 1/(sqrt(n))

∑

xjexp( − i(2πjk)/n)

j = 0

z^k=ak+ibk=1n∑j=0n-1xjexp(-i2πjkn)

and ẑn − kz^n-k is the complex conjugate of ẑkz^k. Thus the DFT of a real sequence is a particular type of complex sequence, called a Hermitian sequence, or half-complex or conjugate symmetric, with the properties

an − k = akbn − k = − bkb0 = 0

an-k=akbn-k=-bkb0=0

and, if nn is even, bn / 2 = 0bn/2=0.

Thus a Hermitian sequence of nn complex data values can be represented by only nn, rather than 2n2n, independent real values. This can obviously lead to economies in storage, with
two schemes
being used in this chapter. In
the first
scheme, which will be referred to as the real storage format for Hermitian sequences, the real parts akak for 0 ≤ k ≤ n / 20≤k≤n/2 are stored in normal order in the first n / 2 + 1n/2+1 locations of an array x of length nn; the corresponding nonzero imaginary parts are stored in reverse order in the remaining locations of x. To clarify, if
x is declared with bounds (0 : n − 1)(0:n-1) in your calling function, the following two tables illustrate the storage of the real and imaginary parts of ẑkz^k for the two cases: nn even and nn odd.

If nn is even then the sequence has two purely real elements and is stored as follows:

Index of x

0

1

2

… …

n / 2n/2

… …

n − 2n-2

n − 1n-1

Sequence

a0a0

a1 + ib1a1+ib1

a2 + ib2a2+ib2

… …

an / 2an/2

… …

a2 − ib2a2-ib2

a1 − ib1a1-ib1

Stored values

a0a0

a1a1

a2a2

… …

an / 2an/2

… …

b2b2

b1b1

x(k) = ak,

for ​k = 0,1, … ,n / 2, and

x(n − k) = bk,

for ​k = 1,2, … ,n / 2 − 1.

xk=ak,for ​k=0,1,…,n/2, andxn-k=bk,for ​k=1,2,…,n/2-1.

If nn is odd then the sequence has one purely real element and, letting n = 2s + 1n=2s+1, is stored as follows:

Index of x

0

1

2

… …

ss

s + 1s+1

… …

n − 2n-2

n − 1n-1

Sequence

a0a0

a1 + ib1a1+ib1

a2 + ib2a2+ib2

… …

as + ibsas+ibs

as − ibsas-ibs

… …

a2 − ib2a2-ib2

a1 − ib1a1-ib1

Stored values

a0a0

a1a1

a2a2

… …

asas

bsbs

… …

b2b2

b1b1

x(k) = ak,

for ​k = 0,1, … ,s, and

x(n − k) = bk,

for ​k = 1,2, … ,s.

xk=ak,for ​k=0,1,…,s, andxn-k=bk,for ​k=1,2,…,s.

The second storage scheme, referred to in this chapter as the complex storage format for Hermitian sequences, stores the real and imaginary parts ak,bkak,bk, for 0 ≤ k ≤ n / 20≤k≤n/2, in consecutive locations of an array x of length n + 2n+2. If x is declared with bounds (0 : n + 1)(0:n+1) in your calling function, the following two tables illustrate the storage of the real and imaginary parts of ẑkz^k for the two cases: nn even and nn odd.

If nn is even then the sequence has two purely real elements and is stored as follows:

Index of x

0

1

2

3

… …

n − 2n-2

n − 1n-1

nn

n + 1n+1

Stored values

a0a0

b0 = 0b0=0

a1a1

b1b1

… …

an / 2 − 1an/2-1

bn / 2 − 1bn/2-1

an / 2an/2

bn / 2 = 0bn/2=0

x(2 × k) = ak,

for ​k = 0,1, … ,n / 2, and

x(2 × k + 1) = bk,

for ​k = 0,1, … ,n / 2.

x2×k=ak,for ​k=0,1,…,n/2, andx2×k+1=bk,for ​k=0,1,…,n/2.

If nn is odd then the sequence has one purely real element and, letting n = 2s + 1n=2s+1, is stored as follows:

For real data that is two-dimensional or higher, the symmetry in the transform persists for the leading dimension only. So, using the notation of equation (3) for the complex two-dimensional discrete transform, we have that ẑk1k2z^k1k2 is the complex conjugate of ẑ(n1 − k1)k2z^(n1-k1)k2. It is more convenient for transformed data of two or more dimensions to be stored as a complex sequence of length (n1 / 2 + 1) × n2 × … × nd(n1/2+1)×n2×…×nd where dd is the number of dimensions. The inverse discrete Fourier transform operating on such a complex sequence (Hermitian in the leading dimension) returns a real array of full dimension (n1 × n2 × … × ndn1×n2×…×nd).

Real symmetric transforms

In many applications the sequence xjxj will not only be real, but may also possess additional symmetries which we may exploit to reduce further the computing time and storage requirements. For example, if the sequence xjxj is odd, (xj = − xn − j)(xj=-xn-j), then the discrete Fourier transform of xjxj contains only sine terms. Rather than compute the transform of an odd sequence, we define the sine transform of a real sequence by

n − 1

x̂k = sqrt(2/n)

∑

xjsin((πjk)/n),

j = 1

x^k=2n∑j=1n-1xjsin(πjkn),

which could have been computed using the Fourier transform of a real odd sequence of length 2n2n. In this case the xjxj are arbitrary, and the symmetry only becomes apparent when the sequence is extended. Similarly we define the cosine transform of a real sequence by

x̂k = sqrt(2/n)

(

n − 1

)

(1/2)x0 +

∑

xjcos((πjk)/n) + (1/2)( − 1)kxn

j = 1

x^k=2n(12x0+∑j=1n-1xjcos(πjkn)+12(-1)kxn)

which could have been computed using the Fourier transform of a real even sequence of length 2n2n.

In addition to these ‘half-wave’ symmetries described above, sequences arise in practice with ‘quarter-wave’ symmetries. We define the quarter-wave sine transform by

x̂k = 1/(sqrt(n))

(

n − 1

)

∑

xjsin((πj(2k − 1))/(2n)) + (1/2)( − 1)k − 1xn

j = 1

x^k=1n(∑j=1n-1xjsin(πj(2k-1)2n)+12(-1)k-1xn)

which could have been computed using the Fourier transform of a real sequence of length 4n4n of the form

(0,x1, … ,xn,xn − 1, … ,x1,0, − x1, … , − xn, − xn − 1, … , − x1).

(0,x1,…,xn,xn-1,…,x1,0,-x1,…,-xn,-xn-1,…,-x1).

Similarly we may define the quarter-wave cosine transform by

x̂k = 1/(sqrt(n))

(

n − 1

)

(1/2)x0 +

∑

xjcos((πj(2k − 1))/(2n))

j = 1

x^k=1n(12x0+∑j=1n-1xjcos(πj(2k-1)2n))

which could have been computed using the Fourier transform of a real sequence of length 4n4n of the form

If the function f(t)f(t) is defined over some more general interval (a,b)(a,b), then the integral transform can still be approximated by the discrete transform provided a shift is applied to move the point aa to the origin.

Convolutions and correlations

One of the most important applications of the discrete Fourier transform is to the computation of the discrete convolution or correlation of two vectors xx and yy defined (as in Brigham (1974)) by

convolution: zk = ∑ j = 0n − 1xjyk − jzk=∑j=0n-1xjyk-j

correlation: wk = ∑ j = 0n − 1xjyk + jwk=∑j=0n-1x-jyk+j

(Here xx and yy are assumed to be periodic with period nn.)

Under certain circumstances (see Brigham (1974)) these can be used as approximations to the convolution or correlation integrals defined by

Applications to solving partial differential equations (PDEs)

A further application of the fast Fourier transform, and in particular of the Fourier transforms of symmetric sequences, is in the solution of elliptic PDEs. If an equation is discretized using finite differences, then it is possible to reduce the problem of solving the resulting large system of linear equations to that of solving a number of tridiagonal systems of linear equations. This is accomplished by uncoupling the equations using Fourier transforms, where the nature of the boundary conditions determines the choice of transforms – see Section [Application to Elliptic Partial Differential Equations]. Full details of the Fourier method for the solution of PDEs may be found in Swarztrauber (1977) and Swarztrauber (1984).

Inverse Laplace Transforms

Let f(t)f(t) be a real function of tt, with f(t) = 0f(t)=0 for t < 0t<0, and be piecewise continuous and of exponential order αα, i.e.,

|f(t)| ≤ Meαt

|f(t)|≤Meαt

for large tt, where αα is the minimal such exponent.

The Laplace transform of f(t)f(t) is given by

∞

F(s) =

∫

e − stf(t)dt, t > 0

0

F(s)=∫0∞e-stf(t)dt, t>0

where F(s)F(s) is defined for Re(s) > αRe(s)>α.

The inverse transform is defined by the Bromwich integral

a + i∞

f(t) = 1/(2πi)

∫

estF(s)ds, t > 0.

a − i∞

f(t)=12πi∫a-i∞a+i∞estF(s)ds, t>0.

The integration is performed along the line s = as=a in the complex plane, where a > αa>α. This is equivalent to saying that the line s = as=a lies to the right of all singularities of F(s)F(s). For this reason, the value of αα is crucial to the correct evaluation of the inverse. It is not essential to know αα exactly, but an upper bound must be known.

The problem of determining an inverse Laplace transform may be classified according to whether (a) F(s)F(s) is known for real values only, or (b) F(s)F(s) is known in functional form and can therefore be calculated for complex values of ss. Problem (a) is very ill-defined and no functions are provided. Two methods are provided for problem (b).

This may be used to compute the sum of the series. For further reading, see Hamming (1962).

Acceleration of Convergence

This device has applications in a large number of fields, such as summation of series, calculation of integrals with oscillatory integrands (including, for example, Hankel transforms), and root-finding. The mathematical description is as follows. Given a sequence of values {sn}{sn}, for n = m, … ,m + 2ln=m,…,m+2l, then, except in certain singular cases, parameters, aa, bibi, cici may be determined such that

l

sn = a +

∑

bicin.

i = 1

sn=a+∑i=1lbicin.

If the sequence {sn}{sn} converges, then aa may be taken as an estimate of the limit. The method will also find a pseudo-limit of certain divergent sequences – see Shanks (1955) for details.

To use the method to sum a series, the terms snsn of the sequence should be the partial sums of the series, e.g., sn = ∑ k = 1ntksn=∑k=1ntk, where tktk is the kkth term of the series. The algorithm can also be used to some advantage to evaluate integrals with oscillatory integrands; one approach is to write the integral (in this case over a semi-infinite interval) as

∞

a1

a2

a3

∫

f(x)dx =

∫

f(x)dx +

∫

f(x)dx +

∫

f(x)dx + …

0

0

a1

a2

∫0∞f(x)dx=∫0a1f(x)dx+∫a1a2f(x)dx+∫a2a3f(x)dx+…

and to consider the sequence of values

a1

a2

a2

s1 =

∫

f(x)dx, s2 =

∫

f(x)dx = s1 +

∫

f(x)dx, etc.,

0

0

a1

s1=∫0a1f(x)dx, s2=∫0a2f(x)dx=s1+∫a1a2f(x)dx, etc.,

where the integrals are evaluated using standard quadrature methods. In choosing the values of the akak, it is worth bearing in mind that nag_sum_accelerate (c06ba) converges much more rapidly for sequences whose values oscillate about a limit. The akak should thus be chosen to be (close to) the zeros of f(x)f(x), so that successive contributions to the integral are of opposite sign. As an example, consider the case where f(x) = M(x)sinxf(x)=M(x)sin⁡x and M(x) > 0M(x)>0: convergence will be much improved if ak = kπak=kπ rather than ak = 2kπak=2kπ.

Recommendations on Choice and Use of Available Functions

One-dimensional Fourier Transforms

The choice of function is determined first of all by whether the data values constitute a real, Hermitian or general complex sequence. It is wasteful of time and storage to use an inappropriate function.

The choice is next determined by how you prefer to store complex data. Real storage format functions store general complex sequences and Hermitian sequences in real arrays. In the case of general complex sequences, the real and imaginary parts are stored separately in two real arrays. Since Hermitian sequences contain some symmetries, these can be stored in a compact form in a single real array. Alternatively, complex storage format functions store the corresponding sequence as a complex array for general sequences, and with real and imaginary parts in contiguous locations of a real array for Hermitian sequences.

Two groups, each of three functions, are provided in real storage format and three groups of two functions are provided in complex storage format.

Group 2 and Group 3
functions are all designed to perform several transforms in a single call, all with the same value of nn. They are likely to be much faster than the Group 1 functions on modern acrchitectures. They do however require more working storage.
It is therefore recommended that, for real storage format, Group 2 functions be used in preference to Group 1 functions, even when only one transform is to be performed.
Group 2 and Group 3 functions differ in the way sequences are stored: Group 2 functions store sequences as rows of a two-dimensional array while Group 3 functions store sequences as columns of a two-dimensional array.
Group 2 and Group 3
functions impose no practical restrictions on the value of nn; however, the fast Fourier transform algorithm ceases to be ‘fast’ if applied to values of nn which cannot be expressed as a product of small prime factors. All the above functions are particularly efficient if the only prime factors of nn are 22, 33 or 55.

If extensive use is to be made of these functions and you are concerned about efficiency, you are advised to conduct your own timing tests.

To compute inverse (backward) discrete Fourier transforms the
real storage format
functions should be used in conjunction with the complex conjugate of a Hermitian or general sequence of complex data values.
In the case of complex storage format functions, there is a direction parameter which determines the direction of the transform; a call to such a function in the forward direction followed by a call in the backward direction reproduces the original data.

Real data

The transform of multidimensional real data is stored as a complex sequence that is Hermitian in its leading dimension. The inverse transform takes such a complex sequence and computes the real transformed sequence. Consequently, separate functions are provided for performing forward and inverse transforms.

The complex sequences computed by nag_sum_fft_real_2d (c06pv) and nag_sum_fft_real_3d (c06py) contain roughly half of the Fourier coefficients; the remainder can be reconstructed by conjugation of those computed. For example, the Fourier coefficients of the two-dimensional transform ẑ(n1 − k1)k2z^(n1-k1)k2 are the complex conjugate of ẑk1k2z^k1k2 for k1 = 0,1, … ,n1 / 2k1=0,1,…,n1/2, and k2 = 0,1, … ,n2 − 1k2=0,1,…,n2-1.

Typically the inversion of a Laplace transform becomes harder as tt increases so that all numerical methods tend to have a limit on the range of tt for which the inverse f(t)f(t) can be computed. nag_sum_invlaplace_crump (c06la) is useful for small and moderate values of tt.

It is often convenient or necessary to scale a problem so that αα is close to 00. For this purpose it is useful to remember that the inverse of F(s + k)F(s+k) is exp( − kt)f(t)exp(-kt)f(t). The method used by nag_sum_invlaplace_crump (c06la) is not so satisfactory when f(t)f(t) is close to zero, in which case a term may be added to F(s)F(s), e.g., k / s + F(s)k/s+F(s) has the inverse k + f(t)k+f(t).

Singularities in the inverse function f(t)f(t) generally cause numerical methods to perform less well. The positions of singularities can often be identified by examination of F(s)F(s). If F(s)F(s) contains a term of the form exp( − ks) / sexp(-ks)/s then a finite discontinuity may be expected in the inverse at t = kt=k. nag_sum_invlaplace_crump (c06la), for example, is capable of estimating a discontinuous inverse but, as the approximation used is continuous, Gibbs' phenomena (overshoots around the discontinuity) result. If possible, such singularities of F(s)F(s) should be removed before computing the inverse.