hide long namesshow long names
hide short namesshow short names
Integer type:  int32  int64  nag_int  show int32  show int32  show int64  show int64  show nag_int  show nag_int

PDF version (NAG web site, 64-bit version, 64-bit version)
Chapter Contents
Chapter Introduction
NAG Toolbox

NAG Toolbox Chapter Introduction

C06 — summation of series

Scope of the Chapter

This chapter is concerned with the following tasks.
(a) 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 seuqnce 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 n complex numbers zj , for j=0,1,,n-1. The direct transform is defined by
z^k = 1n j=0 n-1 zj exp -i 2πjk n (1)
for k=0,1,,n-1 . Note that equation (1) makes sense for all integral k and with this extension z^k  is periodic with period n, i.e., z^k = z^ k±n , and in particular z^ -k = z^ n-k . Note also that the scale-factor of 1n  may be omitted in the definition of the DFT, and replaced by 1n  in the definition of the inverse.
If we write zj = xj + i yj  and z^k = ak + i bk , then the definition of z^k  may be written in terms of sines and cosines as
ak = 1n j=0 n-1 xj cos 2πjk n + yj sin 2πjk n  
bk = 1n j= 0 n- 1 yj cos 2πjk n - xj sin 2πjk n .  
The original data values zj  may conversely be recovered from the transform z^k  by an inverse discrete Fourier transform:
zj = 1n k=0 n-1 z^k exp +i 2πjk n (2)
for j=0,1,,n-1 . If we take the complex conjugate of (2), we find that the sequence z-j  is the DFT of the sequence z^- k . Hence the inverse DFT of the sequence z^k  may be obtained by taking the complex conjugates of the z^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
z^ k1 k2 = 1 n1 n2 j1=0 n1-1 j2=0 n2-1 z j1 j2 exp -i 2 π j1 k1 n1 exp -i 2 π j2 k2 n2 . (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=xj , then
z^k = ak + i bk = 1n j=0 n-1 xj exp -i 2πjk n  
and z^ n-k  is the complex conjugate of z^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
a n-k = ak b n-k = -bk b0 = 0  
and, if n is even, b n/2 = 0 .
Thus a Hermitian sequence of n complex data values can be represented by only n, rather than 2n , independent real values. This can obviously lead to economies in storage, with two schemes being used in this chapter. In the first (deprecated) scheme, which will be referred to as the real storage format for Hermitian sequences, the real parts ak  for 0 k n/2  are stored in normal order in the first n/2 + 1  locations of an array x of length n; 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  in your calling function, the following two tables illustrate the storage of the real and imaginary parts of z^k  for the two cases: n even and n odd.
If n is even then the sequence has two purely real elements and is stored as follows:
Index of x 0 1 2 n/2 n-2 n-1
Sequence a0 a1 + ib1 a2 + ib2 a n/2 a2 - ib2 a1 - ib1
Stored values a0 a1 a2 a n/2 b2 b1
xk = ak , for ​ k= 0, 1, , n/2 , and xn-k = bk , for ​ k= 1, 2, , n/2-1 .  
If n is odd then the sequence has one purely real element and, letting n=2s+1 , is stored as follows:
Index of x 0 1 2 s s+1 n-2 n-1
Sequence a0 a1 + ib1 a2 + ib2 as + i bs as - i bs a2 - ib2 a1 - ib1
Stored values a0 a1 a2 as bs b2 b1
xk = ak , for ​ k= 0, 1, , s , and xn-k = bk , for ​ k= 1, 2, , s .  
The second (recommended) storage scheme, referred to in this chapter as the complex storage format for Hermitian sequences, stores the real and imaginary parts ak , bk , for 0 k n/2 , in consecutive locations of an array x of length n+2 . If x is declared with bounds 0:n+1  in your calling function, the following two tables illustrate the storage of the real and imaginary parts of z^k  for the two cases: n even and n odd.
If n is even then the sequence has two purely real elements and is stored as follows:
Index of x 0 1 2 3 n-2 n-1 n n+1
Stored values a0 b0=0 a1 b1 a n/2-1 b n/2-1 a n/2 b n/2 = 0
x2×k = ak , for ​ k= 0, 1, , n/2 , and x2×k+1 = bk , for ​ k= 0, 1, , n/2 .  
If n is odd then the sequence has one purely real element and, letting n=2s+1 , is stored as follows:
Index of x 0 1 2 3 n-2 n-1 n n+1
Stored values a0 b0=0 a1 b1 b s-1 as bs 0
x2×k = ak , for ​ k= 0, 1, , s , and x2×k+1 = bk , for ​ k= 0, 1, , s .  
Also, given a Hermitian sequence, the inverse (or backward) discrete transform produces a real sequence. That is,
xj = 1n a0 + 2 k=1 n/2-1 ak cos 2πjk n - bk sin 2πjk n + a n/2  
where a n/2 = 0  if n is odd.
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^k1k2 is the complex conjugate of z^n1-k1k2. 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 where d 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××nd).

Real symmetric transforms

In many applications the sequence xj  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 xj  is odd, xj = -x n-j , then the discrete Fourier transform of xj  contains only sine terms. Rather than compute the transform of an odd sequence, we define the sine transform of a real sequence by
x^k = 2n j=1 n-1 xj sin πjkn ,  
which could have been computed using the Fourier transform of a real odd sequence of length 2n . In this case the xj  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 = 2n 12 x0 + j=1 n-1 xj cos πjkn + 12 -1k xn  
which could have been computed using the Fourier transform of a real even sequence of length 2n .
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 = 1n j=1 n-1 xj sin π j 2k-1 2n + 12 -1 k-1 xn  
which could have been computed using the Fourier transform of a real sequence of length 4n  of the form
0,x1,,xn,xn-1 ,,x1,0,-x1,,-xn, -x n-1 ,, -x 1 .  
Similarly we may define the quarter-wave cosine transform by
x^k = 1n 12 x0 + j= 1 n- 1 xj cos π j 2k- 1 2n  
which could have been computed using the Fourier transform of a real sequence of length 4n  of the form
x0,x1,, x n-1 ,0, -x n-1 ,,-x0,-x1,, -x n-1 ,0, x n-1 ,,x1 .  

Fourier integral transforms

The usual application of the discrete Fourier transform is that of obtaining an approximation of the Fourier integral transform
F s = - f t exp -i 2 π s t dt  
when ft  is negligible outside some region 0,c . Dividing the region into n equal intervals we have
F s cn j=0 n-1 fj exp -i 2 π s j c n  
and so
Fk cn j= 0 n- 1 fj exp -i 2 π jk n  
for k=0,1,,n-1 , where fj = f jc/n  and Fk = F k/c .
Hence the discrete Fourier transform gives an approximation to the Fourier integral transform in the region s=0  to s=n/c .
If the function ft  is defined over some more general interval a,b , then the integral transform can still be approximated by the discrete transform provided a shift is applied to move the point a 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 x and y defined (as in Brigham (1974)) by (Here x and y are assumed to be periodic with period n.)
Under certain circumstances (see Brigham (1974)) these can be used as approximations to the convolution or correlation integrals defined by
z s = - x t y s-t dt  
and
w s = - x- t y s+t dt ,   - < s < .  
For more general advice on the use of Fourier transforms, see Hamming (1962); more detailed information on the fast Fourier transform algorithm can be found in Gentleman and Sande (1966) and Brigham (1974).

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 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 ft  be a real function of t, with ft=0  for t<0 , and be piecewise continuous and of exponential order α, i.e.,
ft M eαt  
for large t, where α is the minimal such exponent.
The Laplace transform of ft  is given by
F s = 0 e-st f t dt ,   t>0  
where Fs  is defined for Res > α .
The inverse transform is defined by the Bromwich integral
f t = 12πi a-i a+i est F s ds ,   t>0 .  
The integration is performed along the line s=a  in the complex plane, where a>α . This is equivalent to saying that the line s=a  lies to the right of all singularities of Fs . 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) Fs  is known for real values only, or (b) Fs  is known in functional form and can therefore be calculated for complex values of s. Problem (a) is very ill-defined and no functions are provided. Two methods are provided for problem (b).

Direct Summation of Orthogonal Series

For any series of functions ϕi  which satisfy a recurrence
ϕr+1 x + αr x ϕr x + βr x ϕr-1 x =0  
the sum
r= 0 n ar ϕr x  
is given by
r=0 n ar ϕr x = b0 x ϕ0 x + b1 x ϕ1 x + α0 x ϕ0 x  
where
br x + αr x br+ 1 x + βr+ 1 x br+ 2 x = ar b n+ 1 x = b n+ 2 x = 0 .  
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 , for n=m,,m+2l, then, except in certain singular cases, arguments, a, bi , ci  may be determined such that
sn = a + i=1 l bi cin .  
If the sequence sn  converges, then a 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 sn  of the sequence should be the partial sums of the series, e.g., sn = k=1 n tk , where tk  is the kth 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
0 f x dx = 0a1 f x dx + a1 a2 f x dx + a2 a3 f x dx +  
and to consider the sequence of values
s1 = 0 a1 f x dx ,   s2 = 0 a2 f x dx = s1 + a1 a2 f x dx , etc.,  
where the integrals are evaluated using standard quadrature methods. In choosing the values of the ak , it is worth bearing in mind that nag_sum_accelerate (c06ba) converges much more rapidly for sequences whose values oscillate about a limit. The ak  should thus be chosen to be (close to) the zeros of fx , so that successive contributions to the integral are of opposite sign. As an example, consider the case where f x = M x sinx  and Mx>0 : convergence will be much improved if ak = kπ  rather than ak = 2kπ .

Recommendations on Choice and Use of Available Functions

The fast Fourier transform algorithm ceases to be ‘fast’ if applied to values of n which cannot be expressed as a product of small prime factors. All the FFT functions in this chapter are particularly efficient if the only prime factors of n are 2, 3 or 5.

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.

Real and Hermitian data

nag_sum_fft_realherm_1d (c06pa) transforms a single sequence of real data onto (and in-place) a representation of the transformed Hermitian sequence using the complex storage scheme described in Real transforms. nag_sum_fft_realherm_1d (c06pa) also performs the inverse transform using the representation of Hermitian data and transforming back to a real data sequence.
Alternatively, the two-dimensional function nag_sum_fft_real_2d (c06pv) can be used (on setting the second dimension to 1) to transform a sequence of real data onto an Hermitian sequence whose first half is stored in a separate Complex array. The second half need not be stored since these are the complex conjugate of the first half in reverse order. nag_sum_fft_hermitian_2d (c06pw) performs the inverse operation, transforming the the Hermitian sequence (half-)stored in a Complex array onto a separate real array.

Complex data

nag_sum_fft_complex_1d (c06pc) transforms a single complex sequence in-place; it also performs the inverse transform. nag_sum_fft_complex_1d_multi_col (c06ps) transforms multiple complex sequences, each stored sequentially; it also performs the inverse transform on multiple complex sequences. This function is designed to perform several transforms in a single call, all with the same value of n.
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.

Half- and Quarter-wave Transforms

Four functions are provided for computing fast Fourier transforms (FFTs) of real symmetric sequences. nag_sum_fft_sine (c06re) computes multiple Fourier sine transforms, nag_sum_fft_cosine (c06rf) computes multiple Fourier cosine transforms, nag_sum_fft_qtrsine (c06rg) computes multiple quarter-wave Fourier sine transforms, and nag_sum_fft_qtrcosine (c06rh) computes multiple quarter-wave Fourier cosine transforms.

Application to Elliptic Partial Differential Equations

As described in Applications to solving partial differential equations (PDEs), Fourier transforms may be used in the solution of elliptic PDEs.
nag_sum_fft_sine (c06re) may be used to solve equations where the solution is specified along the boundary.
nag_sum_fft_cosine (c06rf) may be used to solve equations where the derivative of the solution is specified along the boundary.
nag_sum_fft_qtrsine (c06rg) may be used to solve equations where the solution is specified on the lower boundary, and the derivative of the solution is specified on the upper boundary.
nag_sum_fft_qtrcosine (c06rh) may be used to solve equations where the derivative of the solution is specified on the lower boundary, and the solution is specified on the upper boundary.
For equations with periodic boundary conditions the full-range Fourier transforms computed by nag_sum_fft_realherm_1d (c06pa) are appropriate.

Multidimensional Fourier Transforms

The following functions compute multidimensional discrete Fourier transforms of real, Hermitian and complex data stored in complex arrays:
  double Hermitian complex
2 dimensions nag_sum_fft_real_2d (c06pv) nag_sum_fft_hermitian_2d (c06pw) nag_sum_fft_complex_2d (c06pu)
3 dimensions nag_sum_fft_real_3d (c06py) nag_sum_fft_hermitian_3d (c06pz) nag_sum_fft_complex_3d (c06px)
any number of dimensions     nag_sum_fft_complex_multid (c06pj)
The Hermitian data, either transformed from or being transformed to real data, is compacted (due to symmetry) along its first dimension when stored in Complex arrays; thus approximately half the full Hermitian data is stored.
nag_sum_fft_complex_2d (c06pu) and nag_sum_fft_complex_3d (c06px) should be used in preference to nag_sum_fft_complex_multid (c06pj) for two- and three-dimensional transforms, as they are easier to use and are likely to be more efficient.
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.
nag_sum_fft_real_2d (c06pv) performs the forward two-dimensionsal transform while nag_sum_fft_hermitian_2d (c06pw) performs the inverse of this transform.
nag_sum_fft_real_3d (c06py) performs the forward three-dimensional transform while nag_sum_fft_hermitian_3d (c06pz) performs the inverse of this transform.
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-k1k2 are the complex conjugate of z^k1k2 for k1=0,1,,n1/2, and k2=0,1,,n2-1.

Convolution and Correlation

nag_sum_convcorr_real (c06fk) computes either the discrete convolution or the discrete correlation of two real vectors.
nag_sum_convcorr_complex (c06pk) computes either the discrete convolution or the discrete correlation of two complex vectors.

Inverse Laplace Transforms

Two methods are provided: Weeks' method (nag_sum_invlaplace_weeks (c06lb)) and Crump's method (nag_sum_invlaplace_crump (c06la)). Both require the function Fs  to be evaluated for complex values of s. If in doubt which method to use, try Weeks' method (nag_sum_invlaplace_weeks (c06lb)) first; when it is suitable, it is usually much faster.
Typically the inversion of a Laplace transform becomes harder as t increases so that all numerical methods tend to have a limit on the range of t for which the inverse ft  can be computed. nag_sum_invlaplace_crump (c06la) is useful for small and moderate values of t.
It is often convenient or necessary to scale a problem so that α is close to 0. For this purpose it is useful to remember that the inverse of Fs+k  is exp-kt ft . The method used by nag_sum_invlaplace_crump (c06la) is not so satisfactory when ft  is close to zero, in which case a term may be added to Fs , e.g., k/s + Fs  has the inverse k+ft .
Singularities in the inverse function ft  generally cause numerical methods to perform less well. The positions of singularities can often be identified by examination of Fs . If Fs  contains a term of the form exp-ks/s  then a finite discontinuity may be expected in the inverse at t=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 Fs  should be removed before computing the inverse.

Direct Summation of Orthogonal Series

The only function available is nag_sum_chebyshev (c06dc), which sums a finite Chebyshev series
j=0 n cj Tj x ,   j=0 n cj T 2j x   or   j=0 n cj T 2j+1 x  
depending on the choice of a argument.

Acceleration of Convergence

The only function available is nag_sum_accelerate (c06ba).

Decision Trees

Tree 1: Fourier Transform of Discrete Complex Data

Is the data one-dimensional?   Multiple vectors?   Stored as rows?   nag_sum_fft_complex_1d_multi_row (c06pr)
yesyesyes
  no   no   no
Stored as columns?   nag_sum_fft_complex_1d_multi_col (c06ps)
yes
nag_sum_fft_complex_1d (c06pc)
Is the data two-dimensional?   nag_sum_fft_complex_2d (c06pu)
yes
  no
Is the data three-dimensional?   nag_sum_fft_complex_3d (c06px)
yes
  no
Transform on one dimension only?   nag_sum_fft_complex_multid_1d (c06pf)
yes
  no
Transform on all dimensions?   nag_sum_fft_complex_multid (c06pj)
yes

Tree 2: Fourier Transform of Real Data or Data in Complex Hermitian Form Resulting from the Transform of Real Data

Quarter-wave sine (inverse) transform   nag_sum_fft_qtrsine (c06rg)
yes
  no
Quarter-wave cosine (inverse) transform   nag_sum_fft_qtrcosine (c06rh)
yes
  no
Sine (inverse) transform   nag_sum_fft_sine (c06re)
yes
  no
Cosine (inverse) transform   nag_sum_fft_cosine (c06rf)
yes
  no
Is the data three-dimensional?   Forward transform on real data   nag_sum_fft_real_3d (c06py)
yesyes
  no   no
Inverse transform on Hermitian data   nag_sum_fft_hermitian_3d (c06pz)
yes
Is the data two-dimensional?   Forward transform on real data   nag_sum_fft_real_2d (c06pv)
yesyes
  no   no
Inverse transform on Hermitian data   nag_sum_fft_hermitian_2d (c06pw)
yes
Is the data multi one-dimensional?   Sequences stored by row   nag_sum_fft_realherm_1d_multi_row (c06pp)
yesyes
  no   no
Sequences stored by column   nag_sum_fft_realherm_1d_multi_col (c06pq)
yes
nag_sum_fft_realherm_1d (c06pa)

Functionality Index

Acceleration of convergence nag_sum_accelerate (c06ba)
Convolution or Correlation, 
    complex vectors nag_sum_convcorr_complex (c06pk)
    real vectors, 
        time-saving nag_sum_convcorr_real (c06fk)
Discrete Fourier Transform, 
    multidimensional, 
        complex sequence, 
            complex storage nag_sum_fft_complex_multid (c06pj)
            real storage nag_sum_fft_complex_multid_sep (c06fj)
    multiple half- and quarter-wave transforms, 
        Fourier cosine transforms, 
            simple use nag_sum_fft_real_cosine_simple (c06rb)
        Fourier cosine transforms, simple use nag_sum_fft_cosine (c06rf)
        Fourier sine transforms, 
            simple use nag_sum_fft_real_sine_simple (c06ra)
        Fourier sine transforms, simple use nag_sum_fft_sine (c06re)
        quarter-wave cosine transforms, 
            simple use nag_sum_fft_real_qtrcosine_simple (c06rd)
        quarter-wave cosine transforms, simple use nag_sum_fft_qtrcosine (c06rh)
        quarter-wave sine transforms, 
            simple use nag_sum_fft_real_qtrsine_simple (c06rc)
        quarter-wave sine transforms, simple use nag_sum_fft_qtrsine (c06rg)
    one-dimensional, 
        multiple transforms, 
            complex sequence, 
                complex storage by columns nag_sum_fft_complex_1d_multi_col (c06ps)
                complex storage by rows nag_sum_fft_complex_1d_multi_row (c06pr)
            Hermitian/real sequence, 
                complex storage by columns nag_sum_fft_realherm_1d_multi_col (c06pq)
                complex storage by rows nag_sum_fft_realherm_1d_multi_row (c06pp)
            Hermitian sequence, 
                real storage by rows nag_sum_fft_hermitian_1d_multi_rfmt (c06fq)
            real sequence, 
                real storage by rows nag_sum_fft_real_1d_multi_rfmt (c06fp)
        multi-variable, 
            complex sequence, 
                complex storage nag_sum_fft_complex_multid_1d (c06pf)
                real storage nag_sum_fft_complex_multid_1d_sep (c06ff)
        single transforms, 
            complex sequence, 
                time-saving, 
                    complex storage nag_sum_fft_complex_1d (c06pc)
                    real storage nag_sum_fft_complex_1d_sep (c06fc)
            Hermitian/real sequence, 
                time-saving, 
                    complex storage nag_sum_fft_realherm_1d (c06pa)
            Hermitian sequence, 
                time-saving, 
                    real storage nag_sum_fft_hermitian_1d_rfmt (c06fb)
            real sequence, 
                time-saving, 
                    real storage nag_sum_fft_real_1d_rfmt (c06fa)
    three-dimensional, 
        complex sequence, 
            complex storage nag_sum_fft_complex_3d (c06px)
            real storage nag_sum_fft_complex_3d_sep (c06fx)
        Hermitian/real sequence, 
            complex-to-real nag_sum_fft_hermitian_3d (c06pz)
            real-to-complex nag_sum_fft_real_3d (c06py)
    two-dimensional, 
        complex sequence, 
            complex storage nag_sum_fft_complex_2d (c06pu)
        Hermitian/real sequence, 
            complex-to-real nag_sum_fft_hermitian_2d (c06pw)
            real-to-complex nag_sum_fft_real_2d (c06pv)
Inverse Laplace Transform, 
    Crump's method nag_sum_invlaplace_crump (c06la)
    Weeks' method, 
        compute coefficients of solution nag_sum_invlaplace_weeks (c06lb)
        evaluate solution nag_sum_invlaplace_weeks_eval (c06lc)
Summation of Chebyshev series nag_sum_chebyshev (c06dc)

References

Brigham E O (1974) The Fast Fourier Transform Prentice–Hall
Davies S B and Martin B (1979) Numerical inversion of the Laplace transform: A survey and comparison of methods J. Comput. Phys. 33 1–32
Fox L and Parker I B (1968) Chebyshev Polynomials in Numerical Analysis Oxford University Press
Gentleman W S and Sande G (1966) Fast Fourier transforms for fun and profit Proc. Joint Computer Conference, AFIPS 29 563–578
Hamming R W (1962) Numerical Methods for Scientists and Engineers McGraw–Hill
Shanks D (1955) Nonlinear transformations of divergent and slowly convergent sequences J. Math. Phys. 34 1–42
Swarztrauber P N (1977) The methods of cyclic reduction, Fourier analysis and the FACR algorithm for the discrete solution of Poisson's equation on a rectangle SIAM Rev. 19(3) 490–501
Swarztrauber P N (1984) Fast Poisson solvers Studies in Numerical Analysis (ed G H Golub) Mathematical Association of America
Swarztrauber P N (1986) Symmetric FFT's Math. Comput. 47(175) 323–346
Wynn P (1956) On a device for computing the emSn transformation Math. Tables Aids Comput. 10 91–96

PDF version (NAG web site, 64-bit version, 64-bit version)
Chapter Contents
Chapter Introduction
NAG Toolbox

© The Numerical Algorithms Group Ltd, Oxford, UK. 2009–2015