Computing the real roots of a Fourier series-plus-linear-polynomial: A Chebyshev companion matrix approach

Computing the real roots of a Fourier series-plus-linear-polynomial: A Chebyshev companion matrix approach

Applied Mathematics and Computation 219 (2012) 819–826 Contents lists available at SciVerse ScienceDirect Applied Mathematics and Computation journa...

326KB Sizes 0 Downloads 1 Views

Applied Mathematics and Computation 219 (2012) 819–826

Contents lists available at SciVerse ScienceDirect

Applied Mathematics and Computation journal homepage: www.elsevier.com/locate/amc

Computing the real roots of a Fourier series-plus-linear-polynomial: A Chebyshev companion matrix approach John P. Boyd a,⇑, Burhan A. Sadiq b a b

Department of Atmospheric, Oceanic & Space Science, University of Michigan, 2455 Hayward Avenue, Ann Arbor, MI 48109, United States Applied & Interdisciplinary Mathematics Program, Department of Mathematics, University of Michigan, East Hall, Ann Arbor, MI 48109, United States

a r t i c l e

i n f o

Keywords: Fourier series Trigonometric polynomial Rootfinding Secular Companion matrix Chebyshev polynomials

a b s t r a c t Fourier series often need to be generalized by appending a linear polynomial to the usual series of sines and cosines. The integral of a trigonometric polynomial is one example; another is a time series of climate data where the periodic oscillations of the diurnal and annual cycles are accompanied by a non-periodic trend (global warming). Stock market averages fluctuate about a generally upward trend. Such non-periodic variations with time are commonly called ‘‘secular trends’’. We borrow ‘‘secular’’ to label a truncated Fourier series plus a linear trend as a ‘‘linear’’, secular trigonometric polynomial. Standard Fourier rootfinding methods are wrecked by the extra, nonperiodic term. Therefore, we introduce a new algorithm for computing the zeros of a Fourier polynomial-with-secular-trend. First, we expand the linear secular trigonometric polynomial fN ðtÞ as a truncated Chebyshev series. Because of the special structure, it is easy to calculate a problem-dependent truncation M such that the error of the truncated Chebyshev series is guaranteed to be less than a userspecified tolerance. We then find the roots of the truncated Chebyshev series as the eigenvalues of the Chebyshev companion matrix. This computes all roots, but we explain why the method is not reliable for complex-valued roots unless these are close to the real axis. No a priori information is required of the user except the coefficients of the linear secular trigonometric polynomial. Numerical examples show that 13 decimal place accuracy for real roots is typical. Ó 2012 Elsevier Inc. All rights reserved.

1. Introduction Nothing seems more regular than the cycles of the moon, but earth’s satellite is in fact steadily spiraling away from us; its orbit cannot be described by a pure Fourier series, but only by a Fourier series augmented by a nonperiodic trend. The daily and yearly rhythms of weather and climate are similarly modified by nonperiodic trends. Global warming is familiar to all, but there was in fact a warming trend that peaked around 1000 A. D. when the Vikings built a short-lived settlement in a Labrador that was then balmy, and a cooling trend that chilled Europe into the Little Ice Age. Investors are, except for those who have ‘‘shorted’’ the market, very pleased by the generally upward secular trend in stock market averages. In time series analysis [4,9], a pure Fourier series is really the exception rather the rule. To compute the zero-crossings of a cyclic-with-trend signal, one must craft a rootfinder that can accommodate nonperiodic terms. There is no generally accepted term for such a generalized trigonometric series or polynomial, but nonperiodic trends in cyclic phenomena are said to have a ‘‘secular’’ trend, motivating the following. ⇑ Corresponding author. E-mail addresses: [email protected] (J.P. Boyd), [email protected] (B.A. Sadiq). URL: http://www.engin.umich.edu:/~jpboyd/ (John P. Boyd). 0096-3003/$ - see front matter Ó 2012 Elsevier Inc. All rights reserved. http://dx.doi.org/10.1016/j.amc.2012.06.045

820

J.P. Boyd, B.A. Sadiq / Applied Mathematics and Computation 219 (2012) 819–826

Definition 1 (Linear secular trigonometric polynomial). Let A denote a constant. A truncated Fourier series, augmented by a linear term, of the form

fN ðtÞ  At þ

N X

aj cosðjtÞ þ

j¼0

N X bj sinðjtÞ

ð1Þ

j¼1

is a ‘‘linear secular trigonometric polynomial of degree N’’. Note that the polynomial of degree N contains a total of 2N þ 2 terms. We use ‘‘secular’’ as a synonym for ‘‘nonperiodic’’, a standard meaning for this word in astronomy and time series analysis. To find the zeros of a truncated Fourier series, there are several robust techniques as reviewed by the author in [7]. Regrettably, none of the standard rootfinding methods for pure Fourier series carry over to secular trigonometric polynomials even when the trend is linear. We therefore propose a simple new algorithm to find the zeros of a trigonometric polynomial with a linear, nonperiodic term. The first step is to expand fN ðtÞ as a Chebyshev series. The second is to find the roots by computing the eigenvalues of the Chebyshev companion matrix [7].1 This two-part algorithm is actually a very general strategy that can be applied to any transcendental function that is analytic on the chosen expansion interval. Richardson and Trefethen [14,13], for example, find the zeros of a series of sinc functions, modified by a change of coordinate, by reapproximating the sinc sum as a Chebyshev series. To minimize the workload, the interval is automatically subdivided so that for difficult problems, many small-dimensional eigenproblems are solved instead of a single large matrix eigenvalue calculation as advocated in [6] and standard in the public domain Chebfun system [15,2]. The Chebyshev coefficients of such a mapped sinc series are not known analytically, so their adaptation is the same numerical procedure as employed in Chebfun. What is special here is that the Chebyshev coefficients of the trigonometric functions are known in explicit form as Bessel functions. The explicit coefficients allow us to guarantee that if the Chebyshev expansion is truncated at a certain, easily computable degree MðfN Þ, then the error will be less than the user-chosen tolerance. In the next three sections, we discuss Chebyshev approximation, rigorously bounding the error in truncating a Chebyshev series at degree M, and the Chebyshev companion matrix method for rootfinding. Section 5 contains a couple of numerical illustrations. A familiar physical example of a secular trigonometric polynomial is a time series analysis of weather data such as the temperature at 500 hectopascals during spring. By multiplying the observations by a ‘‘window function’’ which is one on the most of the time interval, but decays smoothly to zero at both endpoints, and then removing the secular trend, it is possible to represent the variation of temperature with time as a rapidly converging Fourier series. However, the temperature cannot be accurately represented by pure Fourier series because the average temperature is steadily rising in spring; the rising trend must be added back to yield a secular trigonometric polynomial as the approximation fN ðtÞ. Thus, the question of whether the temperature fell below some value T 0 during spring is really the question: does the secular trigonometric polynomial fN ðtÞ  T 0 have real zeros? The raw data do not provide as reliable an answer as the Fourier series because the temperature can dip to a local minimum between measurements, which for radiosonde data are usually made only twice a day. To avoid inflicting irrelevant intricacies of window functions, meteorological error bars and so on, our numerical examples are contrived rather than taken from data. Because a secular trigonometric polynomial is not periodic, determining the roots on one interval of width 2p does not trivially yield the zeros on all other periods, as would be true for a non-secular (and therefore periodic) trigonometric polynomial. Our recommended strategy is to compute the roots one period at a time. To find the roots of a polynomial g N ðTÞ on T 2 ½p þ 2Lp; p þ 2Lp where L is an integer we shall dub the ‘‘period index’’ where

g N ðTÞ ¼ AT þ g 0 þ

N N X X aj cosðjTÞ þ bj sinðjTÞ; j¼1

ð2Þ

j¼1

it is only necessary to find the roots of fN ðtÞ on t 2 ½p; p where all the terms in fN are identical to g N except that the constant is modified to a0 ¼ g 0 þ 2pLA. Then T root ¼ t root þ 2pL. It is only necessary to examine a finite number of periods because of k k the following. (We make no claim that these inequalities are the best possible.) Theorem 1 (A sufficient condition for a zero-free interval for a secular trigonometric polynomial). If

ja0 j > jApj þ

N N X X jaj j þ jbj j; j¼1

ð3Þ

j¼1

then the linear secular polynomial fN ðtÞ  At þ

PN

j¼0 aj

cosðjtÞ þ

PN

j¼1 bj

sinðjtÞ has no roots on t 2 ½p; p.

1 I. J. Good, who was one of many independent discovers of the Chebyshev analogue of the Frobenius companion matrix, dubbed his invention the ‘‘colleague’’ matrix, but we see no reason to introduce a new term for the same concept expressed in a different polynomial basis.

821

J.P. Boyd, B.A. Sadiq / Applied Mathematics and Computation 219 (2012) 819–826

Proof. The right-hand side of the inequality is an upper bound on all the spatially-varying terms on the interval. If the magnitude of the constant a0 exceeds this bound, then it is impossible for the varying terms to drag the sum of the polynomial to zero. QED Because the constant grows with L, it follows that all zeros of the polynomial must lie on a finite interval as quantified by the following h

Theorem 2 (An upper bound on the range of the zeros of a secular trigonometric polynomial). Let

Lzerofree  smallest positive integer such that Lzerofree >

ja0 j þ jApj þ

Then there are no real zeros of the linear secular polynomial fN ðtÞ  At þ

PN

j¼1 jaj j

þ

PN

j¼1 jbj j

2pjAj PN

j¼0 aj

cosðjtÞ þ

jtj P pð2Lzerofree  1Þ:

PN

ð4Þ

:

j¼1 bj

sinðjtÞ with

ð5Þ

Proof. Apply the change of coordinate t ¼ t 0 þ 2pL to map t 2 ½p þ 2pL; p þ 2pL to t0 2 ½p; p. When L is sufficiently large to satisfy the inequality, then the previous theorem applies. QED h

2. Chebyshev approximation It is convenient, since the canonical interval for Chebyshev polynomials is x 2 ½1; 1, to define x ¼ t=p. Then the goal is to find the roots on x 2 ½1; 1 of

fN ðtÞ  Bx þ

N N X X aj cosðjpxÞ þ bj sinðjpxÞ; j¼0

ð6Þ

j¼1

where B ¼ Ap. The Chebyshev representation of the constant Bx follows trivially from the identity T 1 ðxÞ  x. The expansion of the trigonometric functions follows from (Section 10.12 of [12]; also [1])

cosðjpxÞ ¼ J 0 ðjpÞ þ 2

1 X ð1Þk J 2k ðjpÞT 2k ðxÞ;

ð7Þ

k¼1

sinðjpxÞ ¼ 2

1 X ð1Þk J 2kþ1 ðjpÞT 2kþ1 ðxÞ;

ð8Þ

k¼0

We then approximate the trigonometric polynomial by a truncated Chebyshev series where the Chebyshev truncation will be denoted by M:

fN 

M X g k T k ðxÞ  gðx; MÞ;

ð9Þ

k¼0

g 0 ¼ a0 þ

N X aj J 0 ðjpÞ;

ð10Þ

j¼1

g1 ¼ B þ 2

N X bj J 1 ðjpÞ;

ð11Þ

j¼1

g 2k ¼ ð1Þk 2

N X aj J 2k ðjpÞ;

k ¼ 1; 2 . . . ;

ð12Þ

j¼1

g 2kþ1 ¼ ð1Þk 2

N X bj J2kþ1 ðjpÞ;

k ¼ 1; 2 . . . ;

ð13Þ

j¼1

3. Choosing the Chebyshev truncation M The most conservative and safe estimate of error is to look at the worst-case,

fNworst ðtÞ 

1 fcosðNptÞ þ sinðNptÞg 2

ð14Þ

and to demand that the M-th term in the Chebyshev series for fNworst ðtÞ, J M ðN pÞ, must be very small. Denoting a small, userchosen error tolerance by tolerance , this requires

822

J.P. Boyd, B.A. Sadiq / Applied Mathematics and Computation 219 (2012) 819–826

jJ M ðNpÞj 6 tolerance  1

! M ¼ roundðNpÞ þ s:

ð15Þ

For fixed argument, Bessel functions oscillate for small degree M but decay very rapidly for large degree. WKB theory applied to the Bessel differential equation shows that the transition point or ‘‘turning point’’ is when degree and argument are equal. At the turning point, J M ðMÞ  0:47=M1=3 for large M, which is not very small. It follows that jJ M ðjpÞj  1 if and only M is sufficiently large compared to jp. Therefore, jJ M ðN pÞj 6 tolerance  1 requires that

M ¼ roundðNpÞ þ s

ð16Þ

for some sufficiently large ‘‘shift’’ s where s is a positive integer. An important question is: how large must the shift be? Since the cost of computing the eigenvalues of the Chebyshev companion matrix is proportional to the cube of matrix size, it is desirable to use the smallest possible Chebyshev truncation M consistent with a small Chebyshev error. When the argument z ¼ N p is sufficiently small compared to the Bessel degree M, the Bessel function can be approximated by the lowest term in its power series. Unfortunately, the ratio of the second term to the first is Oðz2 =4MÞ, which is OðMÞ and therefore large when z  OðMÞ. Thus, alas, the power series is useless, not even correct in magnitude, in the desired regime where the Bessel argument and degree are of the same magnitude. Peter Debye, physics Nobel laureate, derived asymptotic approximations which apply uniformly around the turning points. Unfortunately, these are rather complicated and do not generate simple criteria were choosing Mðtolerance ; NÞ. So the best we can do is to offer an algorithm and a chart. Because jT n ðxÞj 6 1 for x 2 ½1; 1, the error in truncating a Chebyshev series up to and including the term of degree Mis rigorously bounded by the sum of the absolute values of all neglected coefficients. This makes it trivial to prove the following. Theorem 3. The error in truncating the Chebyshev series up to and including the term of degree M for a linear secular trigonometric polynomial fN ðxÞ,

fN ðtÞ  Bx þ

N N X X aj cosðjpxÞ þ bj sinðjpxÞ j¼0

ð17Þ

j¼1

is rigorously bounded by

    M 1 X X   g k T n ðxÞ 6 Ebound ðM; fg k gÞ  jg k j; fN ðpxÞ    n¼Mþ1 k¼0

ð18Þ

where the coefficients required in the bound are

g 2k ¼ ð1Þk 2

N X aj J 2k ðjpÞ; j¼1

N X g 2kþ1 ¼ ð1Þk 2 bj J 2kþ1 ðjpÞ; k ¼ 1; 2 . . .

ð19Þ

j¼1

Since in this instance, the Chebyshev coefficients are given analytically to arbitrary order, it is straightforward to compute a truncation M with a guaranteed error bound. First, choose an error tolerance tolerance and an initial degree M. Compute all the Chebyshev coefficients g k up to and including g M ; this requires ðM þ 1ÞðN þ 1Þ evaluations of Bessel functions2 Next, accumulate the error bound in reverse order:

Ebound ðM  1Þ ¼ g M ;

Ebound ðmÞ ¼ Ebound ðm þ 1Þ þ jg mþ1 j;

m ¼ M  2;

M  3; . . .

ð20Þ

and continue until Ebound ðmÞ > tolerance , implying M ¼ m þ 1 is the smallest truncation such that the bound is smaller than the target error tolerance. Lastly, compute the roots of the linear secular trigonometric polynomial as the eigenvalues of the Chebyshev companion matrix for the Chebyshev series for the polynomial, truncated at Chebyshev degree M. The rate-determining step is the Oð10M 3 Þ cost of computing eigenvalues through the QZ algorithm, which requires no initialization and is the default eigensolver in Matlab. (Note that cheaper OðM 2 Þ but less accessible algorithms are available [3].) The error tolerance tolerance is the user’s free choice, dependent on the needs of the problem. Although Chebyshev interpolation and summation are very well-conditioned, some slow growth in the cumulative effects of sums and matrix–vector multiplies is inevitable, so we suggest choosing j logðtolerance Þj no larger than three-quarters of logðmachineepsilonÞ, that is, tolerance no smaller than 1012 when operating in Matlab, for example. To help in estimating M, the upper limit in (20), we have plotted in Fig. 1 the contours of the error bound for the ‘‘worstcase’’ function fNworst ðxÞ  ð1=2ÞðcosðN pÞ þ sinðN pÞÞ for which the error bound is

Eworst ðN; MÞ 

1 X

J n ðNpÞ:

ð21Þ

n¼Mþ1

The vertical axis is not M itself, but rather the shift s, defined as the excess above the ‘‘turning point’’ M ¼ N p. 2 If desired, the Bessel functions can be evaluated very efficiently by computing J M ðjpÞ and J M1 ðjpÞ, followed by the stable backwards recurrence Jn1 ðjpÞ ¼ ð2n=ðjpÞÞJ n ðjpÞ  Jnþ1 ðjpÞ [1], pg. 384.

J.P. Boyd, B.A. Sadiq / Applied Mathematics and Computation 219 (2012) 819–826

823

4. Chebyshev companion matrix The roots of the linear secular trigonometric polynomial are the eigenvalues, multiplied by p, of the Chebyshev companion matrix as described by the following. Theorem 4 (Chebyshev companion matrix). Let g M ðxÞ denote a polynomial of degree M written in ‘‘Chebyshev form’’ as

gðx; MÞ ¼

M X g j T j ðxÞ:

ð22Þ

j¼0

~ Then all roots of g M ðxÞ, both on and off the canonical expansion interval, x 2 ½1; 1, are eigenvalues of the M  M matrix ~ A whose elements are

8 d ; j ¼ 1; k ¼ 1; 2; . . . ; M; > < 2;k   1 Ajk ¼ 2 dj;kþ1 þ dj;k1 ; j ¼ 2; . . . ; ðM  1Þ; k ¼ 1; 2; . . . ; M; > g : ð1Þ gj1 þ ð1=2Þdk;M1 ; j ¼ M; k ¼ 1; 2; . . . ; M;

ð23Þ

M

where djk is the usual Kronecker delta function such that djk ¼ 0 if j – k while djj ¼ 1 for all j [10,11]. For a quintic polynomial, for example,

  0    ð1=2Þ    0    0    g0  ð1Þ 2g

5

1

0

0

ð1=2Þ

ð1=2Þ

0

0

ð1=2Þ

g1 ð1Þ 2g 5

g2 ð1Þ 2g

     0 0    ð1=2Þ 0 :   0 ð1=2Þ   g3 g4   ð1Þ 2g þ ð1=2Þ ð1Þ 2g 0

5

5

0

ð24Þ

5

As carefully analyzed by Day and Romero [10], the the companion matrix is well conditioned after being ’’balanced’’, but may sometimes yield inaccurate eigenvalues if balancing is skipped. Balancing is automatic in the Matlab eig function. Unfortunately, balancing is not incorporated into the Linear Algebra:Eigenvalues function in Maple, and we have observed inaccurate eigenvalues and sometimes no eigenvalues at all when N > 40. The Chebyshev approximation is only accurate on x 2 ½1; 1, and the error grows exponentially off this interval. Therefore, the only reliable roots are those roots xk that lie on this interval (or very close to it). The algorithm described here robustly generates all real roots of the secular trigonmetric polynomial on t 2 ½p; p:

troot ¼ p xk ; k

xk 2 ½1; 1:

ð25Þ

It is easy to polish the eigenvalues by one or two Newton’s iterations applied to the original secular trigonmetric polynomial: ðmÞ;root

ðmþ1Þ;root

tk

ðmÞ;root

¼ tk



fN ðtk

Þ

ðmÞ;root

dfN =dtðtk

Þ

;

ð26Þ

where m is the iteration number.

P Fig. 1. Isolines of the base-ten logarithm of Eworst ðN; MÞ  1 n¼Mþ1 jJ n ðN pÞj, which is a rigorous upper bound on the error in approximating ðcosðN pÞ þ sinðN pÞÞ=2 by a Chebyshev series truncated at M ¼ roundðN pÞ þ s where s is the vertical axis. The thick contour denotes a bound at 1016 , roughly machine epsilon in Matlab and also in double precision (‘‘binary 64’’) in the IEEE 754 floating point arithmetic standard.

824

J.P. Boyd, B.A. Sadiq / Applied Mathematics and Computation 219 (2012) 819–826

4.1. Root doublets and multiple zeros Borrowing a term from spectroscopy, a pair of closely spaced roots are a ’’root doublet’’. When the roots are separated by a small distance , a roundoff-induced perturbation to the polynomial as small as 2 can modify two closely spaced real roots into a pair of zeros with OðÞ small imaginary parts, or vice versa. It therefore is a good idea to accept roots slightly off the real axis, even if interest is entirely confined to real-valued roots, and replace Newton’s method by Cauchy’s rootfinder, which is based on a local approximation of the polynomial by a quadratic Taylor expansion, able to accommodate closely spaced roots. Since these procedures are thoroughly summarized in [8], we shall discuss them no further here. 5. Numerical illustrations 5.1. Example one The problem is to find the roots of the sixth degree secular trigonometric polynomial

f6 ðtÞ ¼ t=p þ

6 6 X X an cosðntÞ þ bn sinðntÞ; n¼0

ð27Þ

n¼1

a0 ¼ 3;

a1 ¼ 1;

a2 ¼ 3;

a3 ¼ 1;

a4 ¼ 3;

b1 ¼ 3;

b2 ¼ 1;

b3 ¼ 3;

b4 ¼ 1;

b5 ¼ 3;

a5 ¼ 1;

a6 ¼ 3;

ð28Þ

b6 ¼ 1:

ð29Þ

Fig. 2 shows that the six real roots are all computed to better than 14 decimal places by our strategy of approximating the linear secular trigonometric polynomial by a Chebyshev series, and then finding the roots of the latter by the companion matrix method. A nonsecular trigonometric polynomial of the same degree would also have six complex-valued roots, but the Chebyshev expansion is accurate only on or very close to the target interval, t 2 ½p; p, so our method provides no reliable information about complex-valued zeros. 5.2. Polynomials with random coefficients We also experimented with secular trigonometric polynomials of the form

fN ðtÞ ¼ t þ

N X

expðnlÞarand cosðntÞ þ n

n¼0

N X

rand

expðnlÞbn

sinðntÞ;

ð30Þ

n¼1

rand

where arand ; b are random numbers uniformly distributed on the interval x 2 ½1; 1 and l is a positive constant, the ‘‘asymptotic rate of geometric convergence’’ [5]. As explained in Chapter 2 of [5], the Fourier series of functions that are periodic and analytic on the real axis have coefficients that decay exponentially with degree; the parameter l allows us to look for effects caused by differing rates of exponential convergence. Fig. 3 shows that there are in fact no systematic trends with either polynomial degree N or with the geometric rate of convergence l: All real roots are computed to at least 13 decimal places. Fig. 4 is similar to Fig. 2 except that it shows the polynomial, error bound and error for a random secular polynomial of degree 20 with l ¼ 0. Ebound(M))

|Error| in Root

−13

0

10

−8

10

10

−14

10

−16

10

−15

0 10 20 30 40 50 M

10

1

2 3 4 5 Root of f(t)

6

f6(t) 10 0 −10 −20 −3

−2

−1

0 t

1

2

3

Fig. 2. The secular trigonometric polynomial of Example 1 is graphed at the bottom. Ebound ðMÞ is plotted in the upper left; the intersection of the curve with the chosen tolerance 1012 [horizontal dotted line] shows that a Chebyshev truncation of M ¼ 43 is sufficient. The upper right-hand graph shows that the six real roots are all computed with an error of less than 1014 as the zeros of the Chebyshev series approximation to f6 ðtÞ.

825

J.P. Boyd, B.A. Sadiq / Applied Mathematics and Computation 219 (2012) 819–826

Max Error in the set of real roots of fN(t)

−12

Error

10

μ=0.1 μ=.2 μ=.3 μ=.4 μ=.5

−13

10

−14

10

0

10

20

30

40

50

N

N=20, μ=0 0

10

log(Ebound(M)) εtolerance −10

10

0

10 20 30 40 50 60 70 80 90 100 110 120 M

L1 Error in the value of the Root

Fig. 3. Each marker denotes the maximum error in any of the real roots for a given pair of values of the trigonometric polynomial degree N and geometric rate of convergence l. (The Fourier coefficients, though otherwise random, are bounded in magnitude by expðnlÞ where n is the degree of the Fourier term.).

N=20, μ=0

−12

10

−13

10

−14

10

−15

10

−16

10

1 2 3 4 5 6 7 8 9 1011121314151617181920 Root of f(t)

N=20, μ=0 10

fN(t)

5 0 −5 −10

−3

−2

−1

0 t

1

2

3

Fig. 4. Bottom: a linear secular trigonometric polynomial of degree 20 with coefficients chosen randomly from a uniform distribution on ½1; 1. (There is no exponential tapering; a physicist would call this ‘‘white noise.) Upper left: plot of the error truncation bound versus M; the horizontal dotted line marks the chosen error tolerance, 1012 . Upper right: errors in the real roots.

6. Summary We have developed a robust algorithm for computing the real-valued zeros of the linear secular trigonometric polynomial, that is, a trigonometric polynomial with an additional term which is linear in the coordinate t. The method is extremely accurate, typically providing 13 decimal place accuracy for all real roots independent of both the trigonometric polynomial

826

J.P. Boyd, B.A. Sadiq / Applied Mathematics and Computation 219 (2012) 819–826

degree N and the exponential rate of convergence of the Fourier coefficients. It is trivial to generalize the non–periodic term from a linear polynomial to a polynomial of higher degree. The only drawback of the scheme is that the size M of the Chebyshev companion matrix is always larger than N p. In contrast, the Fourier companion matrix whose eigenvalues are the roots of a non-secular trigonometric polynomial is only of size 2N. Thus, it is far more efficient to use the Fourier companion matrix whenever the secular, non-periodic term is absent. Acknowledgment This work was supported by the National Science Foundation through grants OCE 1059703 and DMS-1115277. We thank David Day for his lengthy and detailed review. References [1] M. Abramowitz, I.A. Stegun, Handbook of Mathematical Functions, Dover, New York, 1965. [2] Z. Battles, L.N. Trefethen, An extension of Matlab to continuous functions and operators, SIAM J. Sci. Comput. 25 (2004) 1743–1770. [3] D.A. Bini, L. Gemignani, V.Y. Pan, Fast and stable QR eigenvalue algorithms for generalized companion matrices and secular equations, Numer. Math. 100 (2005) 373–408. [4] G.E.P. Box, G.M. Jenkins, G.C. Reinsel, Time Series Analysis: Forecasting and Control, forth ed., Wiley, New York, 2008. p. 784. [5] J.P. Boyd, Chebyshev and Fourier Spectral Methods, second ed., Dover, Mineola, New York, 2001. p. 665. [6] J.P. Boyd, Computing zeros on a real interval through Chebyshev expansion polynomial rootfinding, SIAM J. Numer. Anal. 40 (2002) 1666–1682. [7] J.P. Boyd, Computing the zeros, maxima and inflection points of Chebyshev, Legendre and Fourier series: Solving transcendental equations by spectral interpolation and polynomial rootfinding, J. Eng. Math. 56 (2006) 203–219. [8] J.P. Boyd, D.H. Gally, Numerical experiments on the accuracy of the Chebyshev- Frobenius companion matrix method for finding the zeros of a truncated series of Chebyshev polynomials, J. Comput. Appl. Math. 205 (2007) 281–295. [9] P. Cowpertwait, A Introductory Time Series with R, Springer, Metcalfe, 2009. p. 272. [10] D. Day, L. Romero, Roots of polynomials expressed in terms of orthogonal polynomials, SIAM J. Numer. Anal. 43 (2005) 1969–1987. [11] I.J. Good, The colleague matrix, a Chebyshev analogue of the companion matrix, Quart. J. Math. 12 (1961) 61–68. [12] F.W.J. Olver, D.W. Lozier, R.F. Boisvert, C.W. Clark (Eds.), NIST Handbook of Mathematical Functions, Cambridge University Press, New York, 2010. [13] M. Richardson, Approximating functions with endpoint singularities, Master’s thesis, Oxford, Oxford, 2010. [14] M. Richardson, L.N. Trefethen, A sinc function analogue of Chebfun, SIAM J. Sci. Comput. 33 (2011) 2519–2535. [15] L.N. Trefethen et al., Chebfun Version 4.0, The Chebfun Development Team, 2011. .