Home / Math and Physics Files / Math / math misc

gautschi

PDF · 92 pages · 546.2 KB

Course notes by Walter Gautschi (not Phil's own work) on computational methods for orthogonal polynomials, with Matlab routines from his OPQ package. Part I covers recurrence coefficients, modified Chebyshev, Stieltjes and Lanczos algorithms, discretization, modification algorithms and Sobolev orthogonal polynomials. Part II treats Gauss-type, Kronrod and Turán quadrature, and Cauchy principal value integrals. Part III covers least squares, moment-preserving splines and slowly convergent series, with exercises.

AI-written summary; may contain errors.

Extracted text (machine-read; may contain errors)

Orthogonal Polynomials, Quadrature, and Approximation: Computational Methods and Software (in Matlab) WALTER GAUTSCHI Abstract Orthogonal polynomials, unless they are classical, requir e special techniques for their computation. One of the central proble ms is to generate the coeﬃcients in the basic three-term recurrence relation they are known to satisfy. There are two general approaches f or do- ing this: methods based on moment information, and discreti zation methods. In the former, one develops algorithms that take as input given moments, or modiﬁed moments, of the underlying measur e and produce as output the desired recurrence coeﬃcients. In the ory, these algorithms yield exact answers. In practice, owing to round ing errors, the results are potentially inaccurate depending on the num erical con- dition of the mapping from the given moments (or modiﬁed mome nts) to the recurrence coeﬃcients. A study of related condition n umbers is therefore of practical interest. In contrast to moment-b ased al- gorithms, discretization methods are basically approxima te methods: one approximates the underlying inner product by a discrete inner product and takes the recurrence coeﬃcients of the correspo nding dis- crete orthogonal polynomials to approximate those of the de sired or- thogonal polynomials. Finding discretizations that yield satisfactory rates of convergence requires a certain amount of skill and c reativity on the part of the user, although general-purpose discretiz ations are available if all else fails. Other interesting problems have as objective the computati on of new orthogonal polynomials out of old ones. If the measure of the new 1 orthogonal polynomials is the measure of the old ones multip lied by a rational function, one talks about modiﬁcation of orthogon al polyno- mials and modiﬁcation algorithms that carry out the transit ion from the old to the new orthogonal polynomials. This enters into a circle of ideas already investigated by Christoﬀel in the 1850s, bu t eﬀective algorithms have been obtained only very recently. They requ ire the computation of Cauchy integrals of orthogonal polynomials — another interesting computational problem. In the 1960s, a new type of orthogonal polynomials emerged — the so-called Sobolev orthogonal polynomials — which are ba sed on inner products involving derivatives. Although they prese nt their own computational challenges, moment-based algorithms and di scretiza- tion methods are still two of the main stocks of the trade. The com- putation of zeros of Sobolev orthogonal polynomials is of pa rticular interest in practice. An important application of orthogonal polynomials is to qu adra- ture, speciﬁcally quadrature rules of the highest algebrai c degree of exactness. Foremost among them is the Gaussian quadrature r ule and its close relatives, the Gauss–Radau and Gauss–Lobatto rul es. More recent extensions are due to Kronrod, who inserts n+1 new nodes into a given n-point Gauss formula, again optimally with respect to degre e of exactness, and to Tur´ an, who allows derivative terms to a ppear in the quadrature sum. When integrating functions having pole s outside the interval of integration, quadrature rules of polynomia l/rational degree of exactness are of interest. Poles inside the interv al of in- tegration give rise to Cauchy principal value integrals, wh ich pose computational problems of their own. Interpreting Gaussia n quadra- ture sums in terms of matrices allows interesting applicati ons to the computation of matrix functionals. In the realm of approximation, orthogonal polynomials, esp ecially discrete ones, ﬁnd use in curve ﬁtting, e.g. in the least squa res ap- proximation of discrete data. This indeed is the problem in w hich orthogonal polynomials (in substance if not in name) ﬁrst ap peared in the 1850s in work of Chebyshev. Sobolev orthogonal polynomi als also had their origin in least squares approximation, when one tr ies to ﬁt si- multaneously functions together with some of their derivat ives. Phys- ically motivated are approximations by spline functions th at preserve as many moments as possible. Interestingly, these also are r elated to orthogonal polynomials via Gauss and generalized Gauss-ty pe quadra- ture formulae. Slowly convergent series whose sum can be exp ressed 2 as a deﬁnite integral naturally invite the application of Ga uss-type quadratures to speed up their convergence. An example are se ries whose general term is expressible in terms of the Laplace tra nsform or its derivative of a known function. Such series occur prom inently in plate contact problems. 3 Table of Contents Part Ia. Orthogonal Polynomials 1. Recurrence coeﬃcients 2. Modiﬁed Chebyshev algorithm 3. Discrete Stieltjes and Lanczos algorithm 4. Discretization methods 5. Cauchy integrals of orthogonal polynomials 6. Modiﬁcation algorithms Part Ib. Sobolev Orthogonal Polynomials 7. Sobolev inner product and recurrence relation 8. Moment-based algorithm 9. Discretization algorithm 10. Zeros Exercises to Part I Part II. Quadrature 11. Gauss-type quadrature formulae •Gauss formula •Gauss–Radau formula •Gauss–Lobatto formula 12. Gauss–Kronrod quadrature 13. Gauss–Tur´ an quadrature 14. Quadrature formulae based on rational functions 15. Cauchy principal value integrals 16. Polynomials orthogonal on several intervals 17. Quadrature estimates of matrix functionals Exercises to Part II 4 Part III. Approximation 18. Polynomial least squares approximation •classical •constrained •in Sobolev spaces 19. Moment-preserving spline approximation •on the positive real line •on a compact interval 20. Slowly convergent series •generated by a Laplace transform or derivative thereof •occurring in plate contact problems Exercises to Part III References •W. Gautschi, Orthogonal polynomials: computation and appr oximation. Numerical Mathematics and Scientiﬁc Computation. Oxford U niversity Press, Oxford, 2004. (This item will be referenced by Ga04.) •A suite of Matlab routines, called OPQ, to be found at the web site http://www.cs.purdue.edu/archives/2002/wxg/codes Each routine can be downloaded individually. •G. Szeg¨ o, Orthogonal polynomials (4th edn), AMS Colloq. Pu bl. 23, Amer. Math. Soc., Providence, RI, 1975. (This item will be re ferenced by Sz75.) •M. Abramowitz and I.A. Stegun (eds), Handbook of mathematic al functions, Dover Publ., New York, 1992. (This item will be re ferenced by AS92.) •I.S. Gradshteyn and I.M. Ryzhik, Tables of integrals, serie s, and products (6th edn), Academic Press, San Diego, CA, 2000. (Th is item will be referenced by GR00.) 5 PART Ia ORTHOGONAL POLYNOMIALS 1 Recurrence coeﬃcients 1.1 Background and Notations Orthogonality is deﬁned with respect to an inner product, wh ich in turn involves a measure of integration, d λ. Anabsolutely continuous measure has the form dλ(t) =w(t)dton [a, b],−∞ ≤ a < b ≤ ∞, where wis referred to as a weight function . Usually, wis positive on ( a, b), in which case d λis said to be a positive measure and [a, b] is called the support of dλ. Adiscrete measure has the form dλN(t) =N/summationdisplay k=1wkδ(t−xk)dt, x 1< x2<· · ·< xN, where δis the Dirac delta function, and usually wk>0. The support of d λN consists of its Nsupport points x1, x2, . . ., x N. For absolutely continuous measures, we make the standing assumption that all moments µr=/integraldisplay Rtrdλ(t), r= 0,1,2, . . ., exist and are ﬁnite. The inner product of two polynomials pandqrelative to the measure d λis then well deﬁned by (p, q)dλ=/integraldisplay Rp(t)q(t)dλ(t), and the norm of a polynomial pby /bardblp/bardbldλ=/radicalbig (p, p)dλ. Orthogonal polynomials relative to the (positive) measure d λare deﬁned by πk(·) =πk(·; dλ) a polynomial of exact degree k, k= 0,1,2, . . ., 6 (πk, π/lscript)dλ/braceleftbigg = 0, k/negationslash=/lscript, >0, k=/lscript. They are uniquely deﬁned up to the leading coeﬃcient, if d λis absolutely continuous, and are called monic if the leading coeﬃcient is equal to 1. For a discrete measure d λN, there are exactly Northogonal polynomials π0, π1, . . ., π N−1.Orthonormal polynomials are deﬁned and denoted by ˜πk(·; dλ) =πk(·; dλ) /bardblπk/bardbldλ, k= 0,1,2, . . . . They satisfy (˜πk,˜π/lscript)dλ=δk,/lscript=/braceleftbigg0, k/negationslash=/lscript, 1, k=/lscript. Examples of measures resp. weight functions are shown in Tab les 1 and 2. The former displays the most important “classical” weigh t functions, the latter the best-known discrete measures. Table 1: “Classical” weight functions d λ(t) =w(t)dt name w(t) support comment Jacobi (1 −t)α(1 +t)β[−1,1] α >−1, β >−1 Laguerre tαe−t[0,∞]α >−1 Hermite |t|2αe−t2[−∞,∞]α >−1 2 Meixner-1 2πe(2φ−π)t|Γ(λ+it)|2[−∞,∞]λ >0, Pollaczek 0 < φ < π 1.2 Three-term recurrence relation For any n(< N−1 if dλ= dλN), the ﬁrst n+1 monic orthogonal polynomials satisfy a three-term recurrence relation (1.1)πk+1(t) = (t−αk)πk(t)−βkπk−1(t), k= 0,1, . . ., n −1, π−1(t) = 0, π 0(t) = 1, 7 Table 2: “Classical” discrete measures d λ(t) =/summationtextM k=0wkδ(t−k)dt name M w k comment discrete N−1 1 Chebyshev Krawtchouk N/parenleftbigN k/parenrightbig pk(1−p)N−k0< p < 1 Charlier ∞ e−aak/k! a¿0 Meixner ∞ck Γ(β)Γ(k+β) k!0< c < 1, β >0 Hahn N/parenleftbigα+k k/parenrightbig/parenleftbigβ+N−k N−k/parenrightbig α >−1, β >−1 where the recurrence coeﬃcients αk=αk(dλ),βk=βk(dλ) are real and positive, respectively. The coeﬃcient β0in (1.1) multiplies π−1= 0, and hence can be arbitrary. For later use, it is convenient to deﬁ ne (1.2) β0=β0(dλ) =/integraldisplay Rdλ(t). The proof of (1.1) is rather simple if one expands πk+1(t)−tπk(t)∈Pk in orthogonal polynomials π0, π1, . . ., π kand observes orthogonality and the obvious, but important, property ( tp, q)dλ= (p, tq)dλof the inner product. As a by-product of the proof, one ﬁnds the formulae of Darboux , (1.3)αk(dλ) =(tπk, πk)dλ (πk, πk)dλ, k= 0,1,2, . . ., βk(dλ) =(πk, πk)dλ (πk−1, πk−1)dλ, k= 1,2, . . . . The second yields (1.4) /bardblπk/bardbl2 dλ=β0β1· · ·βk. Placing the coeﬃcients αkon the diagonal, and√βkon the two side diagonals 8 of a matrix produces what is called the Jacobi matrix of the measure d λ, (1.5) J(dλ) = α0√β1 0√β1α1√β2√β2α2... ...... 0 . It is a real, symmetric matrix of inﬁnite order, in general. I ts principal minor matrix of order nwill be denoted by (1.6) Jn(dλ) =J(dλ)[1:n,1:n]. Noting that the three-term recurrence relation for the orth onormal poly- nomials is (1.7)/radicalbig βk+1˜πk+1(t) = (t−αk)˜πk(t)−√βk˜πk−1(t), k= 0,1,2, . . ., ˜π−1(t) = 0,˜π0(t) = 1/√β0, or, in matrix form, with ˜π(t) = [˜π0(t),˜π1(t), . . .,˜πn−1(t)]T, (1.8) t˜π(t) =Jn(dλ)˜π(t) +/radicalbig βn˜πn(t)en, one sees that the zeros τνof˜πn(·; dλ)are precisely the eigenvalues of Jn(dλ), and˜π(τν)corresponding eigenvectors . This is only one of many reasons why knowledge of the Jacobi matrix, i.e. of the recurrence coeﬃc ients, is of great practical interest. For classical measures as the ones in Ta bles 1 and 2, all recurrence coeﬃcients are explicitly known (cf. Ga04, Tabl es 1.1 and 1.2). In most other cases, they must be computed numerically. In theOPQpackage, routines generating recurrence coeﬃcients have t he syntaxab=name(N), where name identiﬁes the name of the orthogonal poly- nomial and Nis an input parameter specifying the number of αkand of βk desired. There may be additional input parameters. The αs and βs are stored in the N×2 arrayab: α0 β0 α1 β1 ...... αN−1βN−1N∈N. 9 For example, ab=rjacobi(N,a,b) generates the ﬁrst Nrecurrence coeﬃ- cients of the Jacobi polynomial with parameters α=a,β=b. Demo#1 The ﬁrst ten recurrence coeﬃcients for the Jacobi polynomia ls with parameters α=−1 2,β=3 2. The Matlab command, followed by the output, is shown in the bo x below. >> ab=r jacobi(10,-.5,1.5) ab = 6.666666666666666e-01 4.712388980384690e+00 1.333333333333333e-01 1.388888888888889e-01 5.714285714285714e-02 2.100000000000000e-01 3.174603174603174e-02 2.295918367346939e-01 2.020202020202020e-02 2.376543209876543e-01 1.398601398601399e-02 2.417355371900826e-01 1.025641025641026e-02 2.440828402366864e-01 7.843137254901961e-03 2.455555555555556e-01 6.191950464396285e-03 2.465397923875433e-01 5.012531328320802e-03 2.472299168975069e-01 2 Modiﬁed Chebyshev algorithm The ﬁrst 2 nmoments µ0, µ1, . . ., µ 2n−1of a measure d λuniquely determine the ﬁrst nrecurrence coeﬃcients αk(dλ) andβk(dλ),k= 0,1, . . ., n −1. How- ever, the corresponding moment map R2n/mapsto→R2n: [µk]2n−1 k=0/mapsto→[αk, βk]n−1 k=0is severely ill-conditioned when nis large. Therefore, other moment maps must be sought that are better conditioned. One that has been stud ied extensively in the literature is based on modiﬁed moments (2.1) mk=/integraldisplay Rpk(t)dλ(t), k= 0,1,2, . . ., where {pk},pk∈Pk, is a given system of polynomials chosen to be close in some sense to the desired polynomials {πk}. We assume that pk, likeπk, satisﬁes a three-term recurrence relation (2.2)pk+1(t) = (t−ak)pk(t)−bkπk−1(t), k= 0,1,2, . . ., p−1(t) = 0, p 0(t) = 1, but with coeﬃcients ak∈R,bk≥0, that are known. The case ak=bk= 0 yields powers pk(t) =tk, hence ordinary moments µk, which however, as already mentioned, is not recommended. 10 Figure 1: Modiﬁed Chebyshev algorithm, schematically lk k,l mσ σσ 0= =Computing stencil l 00 0 0 0 0 0 01 0 1n 2nl 1,0, l The modiﬁed moment map (2.3) R2n/mapsto→R2n: [mk]2n−1 k=0/mapsto→[αk, βk]n−1 k=0 and related maps have been well studied from the point of view of condi- tioning (cf. Ga04, §2.1.5 and 2.1.6). The maps are often remarkably well- conditioned, especially for measures supported on a ﬁnite i nterval, but can still be ill-conditioned otherwise. An algorithm that implements the map (2.3) is the modiﬁed Chebyshev al- gorithm (cf. Ga04, §2.1.7), which improves on Chebyshev’s original algorithm based on ordinary moments. To describe it, we need the mixed moments (2.4) σk/lscript=/integraldisplay Rπk(t; dλ)p/lscript(t)dλ(t), k, /lscript ≥ −1, which by orthogonality are clearly zero if /lscript < k. Algorithm 1 Modiﬁed Chebyshev algorithm initialization: α0=a0+m1/m0, β 0=m0, σ−1,/lscript= 0, /lscript= 1,2, . . .,2n−2, σ0,/lscript=m/lscript, /lscript= 0,1, . . .,2n−1 11 continuation (if n >1): for k= 1,2, . . ., n −1 do σk/lscript=σk−1,/lscript+1−(αk−1−a/lscript)σk−1,/lscript−βk−1σk−2,/lscript +b/lscriptσk−1,/lscript−1, /lscript=k, k+ 1, . . .,2n−k−1, αk=ak+σk,k+1 σkk−σk−1,k σk−1,/lscript−1, β k=σkk σk−1,k−1 Ifak=bk= 0, Algorithm 1 reduces to Chebyshev’s original algorithm. Figure 1 depicts the trapezoidal array of the mixed moments a nd the computing stencil showing that the circled entry is compute d in terms of the four entries below. The entries in boxes are those used to com pute the αs andβs. TheOPQMatlab command that implements the modiﬁed Chebyshev al- gorithm has the form ab=chebyshev(N,mom,abm) , wheremomis the 1 ×2N array of the modiﬁed moments, and abmthe (2N−1)×2 array of the recur- rence coeﬃcients ak,bkfrom (2.2) needed in Algorithm 1: m0m1m2· · ·m2N−1 mom a0 b0 a1 b1 ...... a2N−2b2N−2 abm If the input parameter abmis omitted, the routine assumes ak=bk= 0 and implements Chebyshev’s original algorithm. Demo#2 “Elliptic” orthogonal polynomials These are orthogonal relative to the measure dλ(t) = [(1 −ω2t2)(1−t2)]−1/2dton [−1,1],0≤ω <1. To apply the modiﬁed Chebyshev algorithm, it seems natural t o employ Chebyshev moments (i.e. pk=Tk, the Chebyshev polynomial of degree k) m0=/integraldisplay1 −1dλ(t), m k=1 2k−1/integraldisplay1 −1Tk(t)dλ(t), k≥1. 12 Their computation is not entirely trivial (cf. Ga04, Exampl e 2.29), but a stable algorithm is available as OPQroutinemmell.m , which for given N generates the ﬁrst 2 Nmodiﬁed moments of d λwithω2being input via the parameter om2. The complete Matlab routine is as follows: function ab=r elliptic(N,om2) abm=rjacobi(2*N-1,-1/2); mom=mmell(N,om2); ab=chebyshev(N,mom,abm) Forom2=.999 and N=40, results produced by the routine are partially shown in the box below. ab = 0 9.682265121100620e+00 0 7.937821421385184e-01 0 1.198676724605757e-01 0 2.270401183698990e-01 0 2.410608787266061e-01 0 2.454285325203698e-01 0 2.473016530297635e-01 0 2.482587060199245e-01 ...... 0 2.499915376529289e-01 0 2.499924312667191e-01 0 2.499932210069769e-01 Clearly, βk→1 4ask→ ∞, which is consistent with the fact that d λbelongs to the Szeg¨ o class (cf. Ga04, p. 12). Convergence, in fact, i s monotone for k≥2. 3 Discrete Stieltjes and Lanczos algorithm Computing the recurrence coeﬃcients of a discrete measure i s a prerequisite for discretization methods to be discussed in the next secti on. Given the measure (3.1) d λN(t) =N/summationdisplay k=1wkδ(t−xk)dt, the problem is to compute αν(dλN),βν(dλN) for all ν≤n−1,n≤N, which will provide access to the discrete orthogonal polyno mials of degrees 13 up to n, or else, to determine the Jacobi matrix JN(dλN), which will provide access to all discrete orthogonal polynomials. There are tw o methods in use, a discrete Stieltjes procedure and the Lanczos algorithm. 3.1 Discrete Stieltjes procedure Since the inner product for the measure (3.1) is a ﬁnite sum, (3.2) ( p, q)dλN=N/summationdisplay k=1wkp(xk)q(xk), Darboux’s formulae (1.3) seem to oﬀer attractive means of co mputing the desired recurrence coeﬃcients, since all inner products ap pearing in these formulae are ﬁnite sums. The only problem is that we do not yet know the orthogonal polynomials πk=πk,Ninvolved. For this, however, we can make use of an idea already expressed by Stieltjes in 1884: combin e Darboux’s formulae with the basic three-term recurrence relation. In deed, when k= 0 we know that π0,N= 1, so that Darboux’s formula for α0(dλN) can be applied, and β0(dλN) is simply the sum of the weights wk. Now that we know α0(dλN), we can apply the recurrence relation (1.1) for k= 0 to compute π1,N(t) fort=xk,k= 1,2, . . ., N . We then have all the information at hand to reapply Darboux’s formulae for α1,Nandβ1,N, which in turn allows us to compute π2,N(t) for all t=xkfrom (1.1). In this manner we proceed until allαν,N,βν,N,ν≤n−1, are determined. If n=N, this will yield the Jacobi matrix JN(dλN). The procedure is quite eﬀective, at least when n/lessmuchN. Asnapproaches N, instabilities may develop, particularly if the support po intsxkof dλNare equally, or nearly equally, spaced. TheOPQroutine implementing Stieltjes’s procedure is called by ab= stieltjes(n,xw) , wheren≤N, andxwis anN×2 array containing the support points and weights of the inner product, x1w1 x2w2 ...... xNwN xw 14 As usual, the recurrence coeﬃcients αν,N,βν,N, 0≤ν≤n−1, are stored in then×2 arrayab. 3.2 Lanczos’s algorithm Lanczos’s algorithm is a general procedure to orthogonally tridiagonalize a given symmetric matrix A. Thus, it ﬁnds an orthogonal matrix Qand a symmetric tridiagonal matrix Tsuch that QTAQ=T. Both QandTare uniquely determined by the ﬁrst column of Q. Given the measure (3.1), it is known that an orthogonal matri xQ∈ R(N+1)×(N+1)exists, with the ﬁrst column being e1= [1,0, . . .,0]T∈RN+1, such that (see Ga04, Corollary to Theorem 3.1) (3.3) QT 1√w1√w2· · ·√wN√w1x10· · · 0√w20 x2· · · 0 ...............√wN0 0 · · · xN Q= 1√β00· · · 0√β0α0√β1· · · 0 0√β1α1· · · 0 ............... 0 0 0 · · ·αN−1 . We are thus in the situation described above, where Ais the matrix displayed on the left and Tthe matrix on the right, the desired Jacobi matrix JN(dλN) bordered by a ﬁrst column and a ﬁrst row containing β0. The computation can be arranged so that only the leading principal minor matr ix of order n+ 1 is obtained. Lanczos’s algorithm in its original form (published in 1950 ) is numeri- cally unstable, but can be stabilized using ideas of Rutisha user (1963). An algorithm and pseudocode, using a sequence of Givens rotati ons to con- struct the matrix Qin (3.3), forms the basis for the OPQMatlab code ab=lanczos(n,xw) , where the input and output parameters have the same meaning as in the routine stieltjes.m . This routine enjoys good stability properties but may be con siderably slower than Stieltjes’s procedure. 15 4 Discretization methods The basic idea is to discretize the given measure d λ, i.e. approximate it by a discrete measure (4.1) d λ(t)≈dλN(t), and then use the recurrence coeﬃcients αk(dλN),βk(dλN) of the discrete measure to approximate αk(dλ),βk(dλ). The former are computed by ei- ther Stieltjes’s procedure or the Lanczos algorithm. The eﬀ ectiveness of the method is crucially tied to the quality of the discretizatio n. We illustrate this by a simple, yet interesting, example. Example 1. Chebyshev weight function plus a constant, w(t) = (1 −t2)−1/2+con [−1,1], c > 0. It suﬃces to approximate the inner product for the weight fun ction w. This can always be done by using appropriate quadrature form ulae. In the case at hand, it is natural to treat the two parts of the weight function separately, indeed to use Gauss–Chebyshev quadrature for t he ﬁrst part and Gauss–Legendre quadrature for the second, (4.2)(p, q)w=/integraldisplay1 −1p(t)q(t)(1−t2)−1/2dt+c/integraldisplay1 −1p(t)q(t)dt ≈N/summationdisplay k=1wCh kp(xCh k)q(xCh k) +cN/summationdisplay k=1wL kp(xL k)q(xL k). Here, xCh k,wCh kare the nodes and weights of the M-point Gauss–Chebyshev quadrature rule, and xL k,wL kthose of the Gauss–Legendre quadrature rule. The discrete measure implied by (4.2) is d λNwithN= 2Mand (4.3) d λN(t) =M/summationdisplay k=1wCh kδ(t−xCh k) +cM/summationdisplay k=1wL kδ(t−xL k). What is attractive about this choice is the fact that the appr oximation in (4.2) is actually an equality whenever the product p·qis a polynomial of degree ≤2M−1. Now if we are interested in computing αk(w),βk(w) for k≤n−1, then the products p·qthat occur in Darboux’s formulae are all of 16 degree ≤2n−1. Therefore, we have equality in (4.2) if n≤M. It therefore suﬃces to take M=nin (4.2) to obtain the ﬁrst nrecurrence coeﬃcients exactly. In general, the quadrature rules will not produce exact resu lts, and N will have to be increased through a sequence of integers unti l convergence occurs. Example 1 illustrates the case of a 2-component discretizat ion. In a gen- eralmultiple-component discretization , the support [ a, b] of dλis decomposed intosintervals, (4.4) [ a, b] =s/uniondisplay j=1[aj, bj], where the intervals [ aj, bj] may or may not be disjoint. The measure d λis then discretized on each interval [ aj, bj] using either a tailor-made quadrature (as in Example 1), or a general-purpose quadrature. For the l atter, a Fej´ er quadrature rule on [ −1,1], suitably transformed to [ aj, bj], has been found useful. (The Fej´ er rule is the interpolatory quadrature fo rmula based on Chebyshev points.) If the original measure d λhas also a discrete component, this component is simply added on. Rather than go into detail s (which are discussed in Ga04, §2.2.4), we present the Matlab implementation, another illustrative example, and a demo. TheOPQroutine for the multiple-component discretization is ab=mcdis (n,eps0,quad,Nmax) , where in addition to the variables abandn, which have the usual meaning, there are three other parameters, eps0: a prescribed accuracy tolerance, quad: the name of a quadrature routine carrying out the discretization on each subinterval if tailor-made (oth erwise,quadgp.m , a general-purpose quadrature routine can be used), Nmax: a maximal allowable value for the discretization parameter N. The decomposition (4.4) is input via themc×2 array AB=a1b1 a2b2 ...... amcbmc, wheremcis the number of components (the sin (4.4)). A discrete component which may possibly be present in d λis input via the array 17 DM=x1y1 x2y2 ...... xmpymp, with the ﬁrst column containing the support points, and the s econd column the associated weights. The number of support points is mp. Bothmcandmp are global variables. Another global variable is iq, which has to be set equal to 1 if the user provides his or her own quadrature routine, an d equal to 0 otherwise. Example 2. The normalized Jacobi weight function plus a discrete mea sure. This is the measure dλ(t) = (βJ 0)−1(1−t)α(1 +t)βdt+p/summationdisplay j=1wjδ(t−xj)dton [−1,1], where βJ 0=/integraldisplay1 −1(1−t)α(1 +t)βdt, α > −1, β > −1. Here, one single component suﬃces to do the discretization, and the obvious choice of quadrature rule is the Gauss–Jacobi N-point quadrature formula to which the discrete component is added on. Similarly as in Exa mple 1, taking N=nyields the ﬁrst nrecurrence coeﬃcients αk(dλ),βk(dλ),k≤n−1, exactly. The global parameters in Matlab are here mc=1,mp=p,iq=1, and AB=−11DM=x1w1 x2w2 ...... xpwp Demo#3 Logistic density function, dλ(t) =e−t (1 +e−t)2dt, t∈R. The discretization is conveniently eﬀected by the quadratu re rule /integraldisplay Rp(t)dλ(t) =/integraldisplay∞ 0p(−t) (1 +e−t)2e−tdt+/integraldisplay∞ 0p(t) (1 +e−t)2e−tdt ≈n/summationdisplay ν=1λL νp(−τL ν) +p(τL ν) (1 +e−τLν)2, 18 where τL k,λL kare the nodes and weights of the N-point Gauss–Laguerre quadrature formula. This no longer produces exact results f orN=n, but converges rapidly as N→ ∞. The exact answers happen to be known, αk(dλ) = 0 by symmetry , β0(dλ) = 1, βk(dλ) =k4π2 4k2−1, k≥1. Numerical results produced by mcdis.m withN=40,eps0=103×eps, along with errors (absolute errors for αk, relative errors for βk) are shown in the box below. The two entries in the last row are the maximum erro rs taken over 0 ≤n≤39. n β n errα errβ 0 1.0000000000(0) 7.18(–17) 3.33(–16) 1 3.2898681337(0) 1.29(–16) 2.70(–16) 6 8.9447603523(1) 4.52(–16) 1.43(–15) 15 5.5578278399(2) 2.14(–14) 0.00(+00) 39 3.7535340252(3) 6.24(–14) 4.48(–15) 6.24(–14) 8.75(–15) 5 Cauchy integrals of orthogonal polynomials 5.1 The Jacobi continued fraction TheJacobi continued fraction associated with the measure d λis (5.1) J=J(t; dλ) =β0 t−α0−β1 t−α1−β2 t−α2−· · ·, where αk=αk(dλ),βk=βk(dλ). From the theory of continued fractions it is readily seen that the nth convergent of Jis (5.2)β0 z−α0−β1 z−α1−· · ·βn−1 z−αn−1=σn(z; dλ) πn(z; dλ), n= 1,2,3, . . ., where πnis the monic orthogonal polynomial of degree n, andσna polynomial of degree n−1 satisfying the same basic three-term recurrence relation as πn, but with diﬀerent starting values, (5.3)σk+1(z) = (z−αk)σk(z)−βkσk−1(z), k= 1,2,3, . . ., σ0(z) = 0, σ 1(z) =β0. 19 Recall that β0=/integraltext Rdλ(t). If we deﬁne σ−1=−1, then (5.3) holds also for k= 0. We have, moreover, (5.4) σn(z) =/integraldisplay Rπn(z)−πn(t) z−tdλ(t), n= 0,1,2, . . . , as can be seen by showing that the integral on the right also sa tisﬁes (5.3). If we deﬁne (5.5) F(z) =F(z; dλ) =/integraldisplay Rdλ(t) z−t to be the Cauchy transform of the measure d λ, and more generally consider (5.6) ρn(z) =ρn(z; dλ) =/integraldisplay Rπn(t) z−tdλ(t), theCauchy integral of the orthogonal polynomial πn, we can give (5.4) the form (5.7) σn(z) =πn(z)F(z)−ρn(z), and hence (5.8)σn(z) πn(z)=F(z)−ρn(z) πn(z). An important result from the theory of the moment problem tel ls us that, whenever the moment problem for d λis determined, then (5.9) lim n→∞σn(z) πn(z)=F(z) for z∈C\[a, b], where [ a, b] is the support of the measure d λ. If [a, b] is a ﬁnite interval, then the moment problem is always determined, and (5.9) is known a sMarkov’s theorem . Note from (5.7) that, since σ−1=−1, we have (5.10) ρ−1(z) = 1, and the sequence {ρn}∞ n=−1satisﬁes the same three-term recurrence relation as{πn}∞ n=−1. As a consequence of (5.8) and (5.9), however, it behaves qui te diﬀerently at inﬁnity, (5.11) lim n→∞ρn(z) πn(z)= 0, 20 which implies that {ρn(z)}is the minimal solution of the three-term recur- rence relation having the initial value (5.10). It is well kn own that a minimal solution of a three-term recurrence relation is uniquely de termined by its starting value, and, moreover, that (5.12)ρn(z) ρn−1(z)=βn z−αn−βn+1 z−αn+1−βn+2 z−αn+2−· · ·, i.e. the successive ratios of the minimal solution are the su ccessive tailsof the Jacobi continued fraction (Pincherele’s theorem). In p articular, by (5.12) forn= 0, and (5.2) and (5.9), (5.13) ρ0(z) =F(z), i.e.ρ0is the Cauchy transform of the measure. We remark that (5.6) is meaningful also for real z=xin (a, b), if the integral is interpreted as a Cauchy principal value integral (5.14) ρn(x) =/integraldisplay R−πn(t; dλ) x−tdλ(t), x∈(a, b), and the sequence {ρn(x)}satisﬁes the basic three-term recurrence relation with initial values (5.15) ρ−1(x) = 1, ρ 0(x) =/integraldisplay R−dλ(t) x−t, but is no longer minimal. 5.2 Continued fraction algorithm This is an algorithm for computing the minimal solution ρn(z),z∈C\[a, b], of the basic three-term recurrence relation. Denote the rat io in (5.12) by (5.16) rn−1=ρn(z) ρn−1(z). Then, clearly, (5.17) rn−1=βn z−αn−rn. 21 If, for some ν≥N, we knew rν, we could apply (5.17) for r=ν, ν−1, . . .,0, and then obtain (5.18) ρn(z) =rn−1ρn−1(z), n= 0,1, . . ., N. Thecontinued fraction algorithm is precisely this algorithm, except that rνis replaced by 0. All quantities generated then depend on ν, which is indicated by a superscript. Algorithm 2 Continued fraction algorithm backward phase; ν≥N: r[ν] ν= 0, r[ν] n−1=βn z−αn−r[ν] n, n=ν, ν−1, . . .,0 forward phase: ρ[ν] −1(z) = 1, ρ[ν] n(z) =r[ν] n−1ρ[ν] n−1(z), n= 0,1, . . ., N It can be shown that, as a consequence of the minimality of {ρn(z)}(cf. Ga04, pp. 114–115), (5.19) lim ν→∞ρ[ν] n(z) =ρn(z), n= 0,1, . . ., N, ifz∈C\[a, b]. Convergence is faster the larger dist( z,[a, b]). To compute ρn(z), it suﬃces to apply Algorithm 2 for a sequence of increasing values of νuntil convergence is achieved to within the desired accuracy. TheOPQcommand implementing this algorithm is [rho,r,nu]=cauchy(N,ab,z,eps0,nu0,numax) where the meanings of the output variables rho,rand input variable abare as shown below. ρ0(z) ρ1(z) ... ρN(z)r0(z) r1(z) ... rN(z)α0 β0 α1 β1 ...... αnumax βnumax rho r ab The input variable eps0 is an error tolerance, the variable nu0a suitable starting values of νin Algorithm 2, which is incremented in steps of, say 5, until the algorithm converges to the accuracy eps0. If convergence does not occur within ν≤numax , an error message is issued. 22 6 Modiﬁcation algorithms By “modiﬁcation” of a measure d λ, we mean here multiplication of d λby a rational function rwhich is positive on the support [ a, b] of dλ. The modiﬁed measure thus is (6.1) d ˆλ(t) =r(t)dλ(t), rrational and r >0 on [a, b]. We are interested in determining the recurrence coeﬃcients ˆαk,ˆβkfor dˆλin terms of the recurrence coeﬃcients αk,βkof dλ. An algorithm that carries out the transition from αk,βkto ˆαk,ˆβkis called a modiﬁcation algorithm . While the passage from the orthogonal polynomials relative to dλto those relative to d ˆλis classical (at least in the case when ris a polynomial), the transition in terms of recurrence coeﬃcients is more recent . It was ﬁrst treated for linear factors in 1971 by Galant. Example 3. Linear factor r(t) =s(t−c),c∈R\[a, b],s=±1. Here, sis a sign factor to make r(t)>0 on ( a, b). Galant’s approach is to determine the Jacobi matrix of d ˆλfrom the Jacobi matrix of d λby means of one step of the symmetric, shifted LR algorithm: by t he choice of s, the matrix s[Jn+1(dλ)−cI] is symmetric positive deﬁnite, hence admits a Choleski decomposition s[Jn+1(dλ)−cI] =LLT, where Lis lower triangular. The Jacobi matrix Jn(dˆλ) is now obtained by reversing the order of the product on the right, adding back t he shift c, and then discarding the last row and column,1 Jn(dˆλ) =/parenleftbig LTL+cI/parenrightbig [1:n,1:n]. Since the matrices involved are tridiagonal, the procedure can be imple- mented by simple nonlinear recurrence relations. These can also be obtained more systematically via Christoﬀel’s theorem and its gener alizations. 1See, e.g. W. Gautschi, “The interplay between classical ana lysis and (numerical) linear algebra—a tribute to Gene H. Golub”, Electron. Trans. Numer . Anal. 13 (2002), 119– 147, where it is also shown how a quadratic factor ( t−c1)(t−c2) can be dealt with by one step of the QR algorithm; see in particular §3.2 and 3.3 23 6.1 Generalized Christoﬀel’s theorem We write (6.2) d ˆλ(t) =u(t) v(t)dλ(t), u(t) =±/lscript/productdisplay λ=1(t−uλ), v(t) =m/productdisplay µ=1(t−vµ), where uλandvµare real numbers outside the support of d λ. The sign of u(t) is chosen so that d ˆλis a positive measure. Christoﬀel’s original theorem (1858) relates to the case v(t) = 1, i.e. m= 0. The generalization to arbitrary vis due to Uvarov (1969). It has a diﬀerent form depending on whether m≤norm > n . In the ﬁrst case, it states that (6.3)u(t)πn(t; dˆλ) = const × /vextendsingle/vextendsingle/vextendsingle/vextendsingle/vextendsingle/vextendsingle/vextendsingle/vextendsingle/vextendsingle/vextendsingle/vextendsingle/vextendsingle/vextendsingle/vextendsingleπn−m(t)· · · πn−1(t)πn(t)· · · πn+/lscript(t) πn−m(u1)· · ·πn−1(u1)πn(u1)· · ·πn+/lscript(u1) · · · · · · · · · · · · · · · · · · πn−m(u/lscript)· · ·πn−1(u/lscript)πn(u/lscript)· · ·πn+/lscript(u/lscript) ρn−m(v1)· · ·ρn−1(v1)ρn(v1)· · ·ρn+/lscript(v1) · · · · · · · · · · · · · · · · · · ρn−m(vm)· · ·ρn−1(vm)ρn(vm)· · ·ρn+/lscript(vm)/vextendsingle/vextendsingle/vextendsingle/vextendsingle/vextendsingle/vextendsingle/vextendsingle/vextendsingle/vextendsingle/vextendsingle/vextendsingle/vextendsingle/vextendsingle/vextendsingle, where ρk(z) =/integraldisplay Rπk(t; dλ) z−tdλ(t), k= 0,1,2, . . . , are the Cauchy integrals of the orthogonal polynomials πk. They occur only ifm >0. To get monic polynomials, the constant in (6.3) must be tak en to be the reciprocal of the (signed) cofactor of the element πn+/lscript(t). Ifm > n , the generalized Christoﬀel theorem has the form (6.4)u(t)πn(t; dˆλ) = const × /vextendsingle/vextendsingle/vextendsingle/vextendsingle/vextendsingle/vextendsingle/vextendsingle/vextendsingle/vextendsingle/vextendsingle/vextendsingle/vextendsingle/vextendsingle/vextendsingle0 0 · · · 0 π0(t)· · · πn+/lscript(t) 0 0 · · · 0 π0(u1)· · ·πn+/lscript(u1) · · · · · · · · · · · · · · · · · · · · · 0 0 · · · 0 π0(u/lscript)· · ·πn+/lscript(u/lscript) 1v1· · ·vm−n−1 1 ρ0(v1)· · ·ρn+/lscript(v1) · · · · · · · · · · · · · · · · · · · · · 1vm· · ·vm−n−1 m ρ0(vm)· · ·ρn+/lscript(vm)/vextendsingle/vextendsingle/vextendsingle/vextendsingle/vextendsingle/vextendsingle/vextendsingle/vextendsingle/vextendsingle/vextendsingle/vextendsingle/vextendsingle/vextendsingle/vextendsingle. Both versions of the theorem remain valid for complex uλ,vµif orthogo- nality is understood in the sense of formal orthogonality . 24 6.2 Linear factors Generalizing Example 3 to arbitrary complex shifts, we let (6.5) d ˆλ(t) = (t−z)dλ(t), z∈C\[a, b]. Using Christoﬀel’s theorem, letting ˆ πn(·) =πn(·; dˆλ), we have (6.6) ( t−z)ˆπn(t) =/vextendsingle/vextendsingle/vextendsingle/vextendsingleπn(t)πn+1(t) πn(z)πn+1(z)/vextendsingle/vextendsingle/vextendsingle/vextendsingle −πn(z)=πn+1(t)−rnπn(t), where (6.7) rn=πn+1(z) πn(z). Following Verlinden (1999), we write ( t−z)tˆπn(t) in two diﬀerent ways: in the ﬁrst, we use the three-term recurrence relation for πkto obtain (t−z)tˆπk(t) =tπk+1(t)−rk·tπk(t) =πk+2(t) + (αk+1−rk)πk+1(t) + (βk+1−rkαk)πk(t)−rkβkπk−1(t); in the second, we use the three-term recurrence relation dir ectly on ˆ πk, and then apply (6.6), to write (t−z)tˆπk(t) = (t−z)[ˆπk+1+ ˆαkˆπk(t) +ˆβkˆπk−1(t)] =πk+2(t) + (ˆαk−rk+1)πk+1(t) + (ˆβk−rkˆαk)πk(t)−rk−1ˆβkπk−1(t). Since orthogonal polynomials are linearly independent, th e coeﬃcients in the two expressions obtained must be the same. This yields ˆαk−rk+1=αk+1−rk, r k−1ˆβk=rkβk, hence the following algorithm. Algorithm 3 Modiﬁcation by a linear factor t−z initialization: r0=z−α0, r 1=z−α1−β1/r0, ˆα0=α1+r1−r0,ˆβ0=−r0β0. 25 continuation (if n >1): for k= 1,2, . . ., n −1 do rk+1=z−αk+1−βk+1/rk, ˆαk=αk+1+rk+1−rk, ˆβk=βkrk/rk−1. Note that this requires αn,βnin addition to the usual nrecurrence coeﬃcients αk,βkfork≤n−1. Algorithm 3 has been found to be numerically stable. TheOPQMatlab command implementing Algorithm 3 is ab=chri1(N,ab0,z) whereab0is an (N+1)×2 array containing the recurrence coeﬃcients αk,βk, k= 0,1, . . .,N. 6.3 Quadratic factor We consider (real) quadratic factors ( t−x)2+y2= (t−z)(t−z),z=x+iy, y >0. Christoﬀel’s theorem is now applied with u1=z,u2=zto express (t−z)(t−z)ˆπn(t) as a linar combination of πn,πn+1, and πn+2, (6.8) ( t−z)(t−z)ˆπn(t) =πn+2(t) +snπn+1(t) +tnπn(t), where (6.9) sn=−/parenleftbigg r/prime n+1+r/prime/prime n+1 r/prime/primenr/prime n/parenrightbigg , t n=r/prime/prime n+1 r/prime/primen|rn|2. Here we use the notation (6.10) r/prime n= Rern(z), r/prime/prime n= Imrn(z),|rn|2=|rn(z)|2, n= 0,1,2, . . ., where rn(z) continues to be the quantity deﬁned in (6.7). The same techn ique used in §6.2 can be applied to (6.8): express ( t−z)(t−z)tˆπk(t) in two diﬀerent ways as a linear combination of πk+3, πk+2, . . ., π k−1and compare the respective coeﬃcients. The result gives rise to 26 Algorithm 4 Modiﬁcation by a quadratic factor ( t−z)(t−z),z=x+ iy initialization: r0=z−α0, r1=z−α1−β1/r0, r2=z−α2−β2/r1, ˆα0=α2+r/prime 2+r/prime/prime 2 r/prime/prime 1r/prime 1−/parenleftbigg r/prime 1+r/prime/prime 1 r/prime/prime 0r/prime 0/parenrightbigg , ˆβ0=β0(β1+|r0|2). continuation (if n >1): for k= 1,2, . . ., n −1 do rk+2=z−αk+2−βk+2/rk+1, ˆαk=αk+2+r/prime k+2+r/prime/prime k+2 r/prime/prime k+1r/prime k+1−/parenleftbigg r/prime k+1+r/prime/prime k+1 r/prime/prime kr/prime k/parenrightbigg , ˆβk=βkr/prime/prime k+1r/prime/prime k−1 [r/prime/prime k]2/vextendsingle/vextendsingle/vextendsingle/vextendsinglerk rk−1/vextendsingle/vextendsingle/vextendsingle/vextendsingle2 . Note that this requires αk,βkforkup to n+ 1. Algorithm 4 is also quite stable, numerically. TheOPQroutine for Algorithm 4 is ab=chri2(N,ab0,x,y) with obvious meanings of the variables involved. Since any real polynomial can be factored into a product of re al linear and quadratic factors of the type considered, Algorithms 3 and 4 can be applied repeatedly to deal with modiﬁcation by an arbitrary polynom ial which is positive on the support [ a, b]. 6.4 Linear divisor In analogy to (6.5), we consider (6.11) d ˆλ(t) =dλ(t) t−z, z∈C\[a, b]. Now the generalized Christoﬀel theorem (with /lscript= 0,m= 1) comes into play, giving (6.12) ˆ πn(t) =/vextendsingle/vextendsingle/vextendsingle/vextendsingleπn−1(t)πn(t) ρn−1(z)ρn(z)/vextendsingle/vextendsingle/vextendsingle/vextendsingle −ρn−1(z)=πn(t)−rn−1πn−1(t), 27 where now (6.13) rn=ρn+1(z) ρn(z). Similarly as in §6.2 and 6.3, we express tˆπk(t) in two diﬀerent ways as a linear combination of πk+1, πk, . . ., π k−2and compare coeﬃcients. By convention, ˆβ0=/integraldisplay Rdˆλ(t) =/integraldisplay Rdλ(t) t−z=−ρ0(z). The result is: Algorithm 5 Modiﬁcation by a linear divisor initialization: ˆα0=α0+r0,ˆβ0=−ρ0(z). continuation (if n >1): for k= 1,2, . . ., n −1 do ˆαk=αk+rk−rk−1, ˆβk=βk−1rk−1/rk−2. Note that here no coeﬃcient αk,βkbeyond k≤n−1 is needed, not even βn−1. The ratios rkof Cauchy integrals that appear in Algorithm 5 can be precomputed by Algorithm 2, where only the backward phase is relevant, convergence being tested on the r[ν] k. Once converged, the algorithm also provides ρ0(z) =r[∞] −1. Aszapproaches the support interval [ a, b], the strength of minimality of the Cauchy integrals {ρk(z)}weakens and ceases altogether when z=x∈ [a, b]. For zvery close to [ a, b], Algorithm 2 therefore converges very slowly. On the other hand, since minimality is very weak, one can gene rateρkwith impunity, if nis not too large, by forward application of the basic three-t erm recurrence relation, using the initial values ρ−1(z) = 1 and ρ0(z). All of this is implemented in the OPQroutine [ab,nu]=chri4(N,ab0,z,eps0,nu0,numax,rho0,iopt) where all variables except rhoandiopthave the same meaning as before. The parameter rhoisρ0(z), whereas ioptcontrols the method of computa- tion for rk: Algorithm 2 if iopt=1, and forward recursion otherwise. 28 6.5 Quadratic divisor We now consider (6.14) d ˆλ(t) =dλ(t) (t−z)(t−z)=dλ(t) (t−x)2+y2, z=x+ iy, x∈R, y > 0. Here we have (6.15) ˆ α0=/integraldisplay Rtdλ(t)/|t−z|2 /integraldisplay Rdλ(t)/|t−z|2=x+yReρ0(z) Imρ0(z),ˆβ0=−1 yImρ0(z). We are in the case /lscript= 0,m= 2 of the generalized Christoﬀel theorems (6.3) and (6.4), which give (6.16) ˆπn(t) =/vextendsingle/vextendsingle/vextendsingle/vextendsingle/vextendsingle/vextendsingleπn−2(t)πn−1(t)πn(t) ρn−2(z)ρn−1(z)ρn(z) ρn−2(z)ρn−1(z)ρn(z)/vextendsingle/vextendsingle/vextendsingle/vextendsingle/vextendsingle/vextendsingle /vextendsingle/vextendsingle/vextendsingle/vextendsingleρn−2(z)ρn−1(z) ρn−2(z)ρn−1(z)/vextendsingle/vextendsingle/vextendsingle/vextendsingle, n≥2; ˆπ1(t) =/vextendsingle/vextendsingle/vextendsingle/vextendsingle/vextendsingle/vextendsingle0π0(t)π1(t) 1ρ0(z)ρ1(z) 1ρ0(z)ρ1(z)/vextendsingle/vextendsingle/vextendsingle/vextendsingle/vextendsingle/vextendsingle /vextendsingle/vextendsingle/vextendsingle/vextendsingle1ρ0(z) 1ρ0(z)/vextendsingle/vextendsingle/vextendsingle/vextendsingle. This becomes (6.17) ˆ πn(t) =πn(t) +snπn−1(t) +tnπn−2(t), n≥1, where (6.18) sn=−/parenleftbigg r/prime n−1+r/prime/prime n−1 r/prime/prime n−2r/prime n−2/parenrightbigg , n≥1;tn=r/prime/prime n−1 r/prime/prime n−2|rn−2|2, n≥2, withrnas deﬁned in (6.13) and notations as in (6.10). Exactly the sa me procedure used to obtain Algorithm 5 yields Algorithm 6 Modiﬁcation by a quadratic divisor initialization: ˆα0=x+ρ/prime 0y/ρ/prime/prime 0,ˆβ0=−ρ/prime/prime 0/y, ˆα1=α1−s2+s1,ˆβ1=β1+s1(α0−ˆα1)−t2, ˆα2=α2−s3+s2,ˆβ2=β2+s2(α1−ˆα2)−t3+t2. 29 continuation (if n >3): for k= 3,4, . . ., n −1 do ˆαk=αk−sk+1+sk,ˆβk=βk−2tk/tk−1. TheOPQroutine for Algorithm 6 is [ab,nu]=chri5(N,ab0,z,eps0,nu0,numax,rho0,iopt) where the input and output variables have the same meaning as in the routine chri4.m . Just like Algorithms 3 and 4, also Algorithms 5 and 6 can be app lied repeatedly to deal with more general polynomial divisors. 30 PART Ib SOBOLEV ORTHOGONAL POLYNOMIALS 7 Sobolev inner product and recurrence re- lation In contrast to the orthogonal polynomials considered so far , the inner product here involves not only function values, but also successive derivative values, all being endowed with their own measures. Thus, (7.1)(p, q)S=/integraldisplay Rp(t)q(t)dλ0(t) +/integraldisplay Rp/prime(t)q/prime(t)dλ1(t) +· · ·+/integraldisplay Rp(s)(t)q(s)(t)dλs(t), s≥1. If all the measures d λσare positive, as we assume, the inner product (7.1) has associated with it a sequence of (monic) polynomials πk(·;S),k= 0,1,2, . . ., orthogonal in the sense (7.2) ( πk, π/lscript)S/braceleftbigg = 0, k/negationslash=/lscript, >0, k=/lscript. These are called Sobolev orthogonal polynomials . We cannot expect them to satisfy a three-term recurrence relation, since the inner p roduct no longer has the shift property ( tp, q) = (p, tq). However, like any sequence of monic polynomials of degrees 0 ,1,2, . . . , they must satisfy an extended recurrence relation of the type (7.3) πk+1(t) =tπk(t)−k/summationdisplay j=0βk jπk−j(t), k= 0,1,2, . . . . Associated with it is the upper Hessenberg matrix of recurrence coeﬃcients (7.4) Hn= β0 0β1 1β2 2· · ·βn−2 n−2βn−1 n−1 1β1 0β2 1· · ·βn−2 n−3βn−1 n−2 0 1 β2 0· · ·βn−2 n−4βn−1 n−3 · · · · · · · · · · · · · · · · · · 0 0 0 · · ·βn−2 0βn−1 1 0 0 0 · · · 1βn−1 0 . 31 In the case s= 0 (of ordinary orthogonal polynomials) there holds βk j= 0 forj >1, and the matrix Hnis tridiagonal. If symmetrized by a (real) diagonal similarity transformation, it becomes the Jacobi matrix Jn(dλ0). When s >0, however, symmetrization of Hnis no longer possible, since Hn may well have complex eigenvalues (see Example 5). 8 Moment-based algorithm There are now s+1 sets of modiﬁed moments, one set for each measure d λσ, (8.1) m(σ) k=/integraldisplay Rpk(t)dλσ, k= 0,1,2, . . .;σ= 0,1, . . ., s. The ﬁrst 2 nmodiﬁed moments of all the sets will uniquely determine the matrix Hnin (7.4), i.e. there is a well-determined map (8.2) [ m(σ) k]2n−1 k=0, σ= 0,1, . . ., s /mapsto→Hn, called modiﬁed moment map for Sobolev orthogonal polynomials. In the case where the polynomials pkin (8.1) satisfy a three-term recurrence relation with known coeﬃcients, and for s= 1, an algorithm has been developed that implements the map (8.2). It very much resembles the modiﬁed Chebyshev algorithm for ordinary orthogonal polynomials, but is tech nically much more elaborate.2The algorithm, however, is implemented in the OPQroutine B=chebyshev sob(N,mom,abm) which produces the N×Nupper triangular matrix Bof recurrence coeﬃcients, withβk j, 0≤j≤k, 0≤k≤N–1, occupying the position ( j+ 1, k+ 1) in the matrix. The input parameter momis the 2 ×(2N) array of modiﬁed moments m(σ) k,k= 0,1, . . .,2N–1;σ= 0,1, of the two measures d λ0and dλ1, andabm the (2N−1)×2 array of coeﬃcients ak,bk,k= 0,1, . . .,2N–2, deﬁning the polynomials pk. Example 3. Althammer’s polynomials (1962) These are the Sobolev polynomials relative to the measures d λ0(t) = d t, dλ1(t) =γdton [−1,1],γ >0. 2See W. Gautschi and M. Zhang,“Computing orthogonal polynom ials in Sobolev spaces”, Numer. Math. 71 (1995), 159–183. 32 A natural choice of modiﬁed moments are the Legendre moments , i.e. pk(t) is the monic Legendre polynomial of degree k. By orthogonality of the Legendre polynomials, all modiﬁed moments m(0) kandm(1) kare zero for k >0, while m(0) 0= 2 and m(1) 0= 2γ. The following Matlab routine, therefore, can be used to generate the Althammer polynomials. mom=zeros(2,2*N); mom(1,1)=2; mom(2,1)=2*g; abm=rjacobi(2*N-1); B=chebyshev sob(N,mom,abm); 9 Discretization algorithm Taking the inner product of both sides of (7.3) with πk−jgives 0 = (πk+1, πk−j)S= (tπk, πk−j)S−βk j(πk−j, πk−j)S, j= 0,1, . . ., k, hence (9.1) βk j=(tπk, πk−j)S (πk−j, πk−j)S, j= 0,1, . . ., k ;k= 0,1, . . ., n −1. These are the analogues of Darboux’s formulae for ordinary o rthogonal poly- nomials, and like these, can be combined with the recurrence relation (7.3) to successively build up the recurrence coeﬃcients βk jin the manner of Stielt- jes’s procedure. The technical details, of course, are more involved, since we must generate not only the polynomials πk, but also their derivatives, in order to be able to compute the Sobolev inner products in (9.1 ). This all is implemented, for arbitrary s≥1, in the Matlab routine stieltjes sob.m . The basic assumption in the design of this routine is the avai lability, for each measure d λσ, of an nσ-point quadrature rule (9.2)/integraldisplay Rp(t) dλσ(t) =nσ/summationdisplay k=1w(σ) kp(x(σ) k), p∈P2(n−σ)−1, σ= 0,1, . . ., s, that is exact for polynomials pof degree ≤2(n−σ)−1. These are typically Gaussian quadrature rules, possibly with discrete compone nts (present in dλσ) added on. The information is supplied to the routine via the 1×(s+1) array 33 nd= [n0, n1, . . ., n s] and themd×(2s+ 2) array xw=x(0) 1· · ·x(s) 1w(0) 1· · ·w(s) 1 x(0) 2· · ·x(s) 2w(0) 2· · ·w(s) 2 ............ x(0) md· · ·x(s) mdw(0) md· · ·w(s) md wheremd=max(nd) . In each column of xwthe entries after x(σ) nσresp.w(σ) nσ(if any) are not used by the routine. Two more input parameters ar e needed; the ﬁrst is a0, the coeﬃcient α0(dλ0), which allows us to initialize the matrix of recurrence coeﬃcients, β0 0=(t,1)S (1,1)S=(t,1)dλ0 (1,1)dλ0=α0(dλ0). The other, same, is a logical variable set equal to 1 if all quadrature rules h ave the same set of nodes, and equal to 0 otherwise. The role of thi s parameter is to switch to a simpliﬁed, and thus faster, procedure if same=1. A call to the routine, therefore, has the form B=stieltjes sob(N,s,nd,xw,a0,same) Example 4. Althammer’s polynomials, revisited. Here, the obvious choice of the quadrature rule for d λ0and d λ1is the n-point Gauss–Legendre rule. This gives rise to the followin g routine: s=1; nd=[N N]; a0=0; same=1; ab=rjacobi(N); zw=gauss(N,ab); xw=[zw(:,1) zw(:,1) ... zw(:,2) g*zw(:,2)]; B=stieltjes sob(N,s,nd,xw,a0,same); The results are identical with those obtained in Example 3. 34 10 Zeros If we let πT(t) = [π0(t), π1(t), . . ., π n−1(t)], where πkare the Sobolev orthog- onal polynomials, then the recurrence relation (7.3) can be written in matrix form as (10.1) tπT(t) =πT(t)Hn+πn(t)eT n in terms of the matrix Hnin (7.4). This immediately shows that the zeros τν ofπnare the eigenvalues of HnandπT(τν) corresponding left eigenvectors. Naturally, there is no guarantee that the eigenvalues are re al; some may well be complex. Also, if nis large, there is a good chance that some of the eigenvalues are ill-conditioned. TheOPQroutine for the zeros of πnis z=sobzeros(n,N,B) whereBis theN×Nmatrix returned by chebyshev sob.m orstieltjes sob.m , andzthen-vector of the zeros of πn, 1≤n≤N. Example 5. Sobolev orthogonal polynomials with only a few real zeros (Mei- jer, 1994). The Sobolev inner product in question is (10.2) (u, v)S=/integraldisplay3 −1u(t)v(t) dt+γ/integraldisplay1 −1u/prime(t)v/prime(t) dt+/integraldisplay3 1u/prime(t)v/prime(t) dt, γ > 0. Meijer proved that for n(even) ≥2 and γsuﬃciently large, the polynomial πn(·;S) has exactly two real zeros, one in [ −3,−1] and the other in [1 ,3]. If n(odd) ≥3, there is exactly one real zero, located in [1 ,3], ifγis suﬃciently large. We use the routine stieltjes sob.m andsobzeros.m to illustrate this for n= 6 and γ= 44,000. (The critical value of γabove which Meijer’s theorem takes hold is about γ= 43,646.2; see Ga04, Table 2.30.) The inner product corresponds to the case s= 1 and dλ0(t) = dton [−1,3],dλ1(t) =/braceleftbigg γdtift∈[−1,1], dtift∈(1,3]. Thus, we can write, with suitable transformations of variab les, /integraldisplay3 −1p(t) dλ0(t) = 2/integraldisplay1 −1p(2x+1) dx,/integraldisplay3 −1p(t) dλ1(t) =/integraldisplay1 −1[γp(x)+p(x+2)] d x 35 and apply n-point Gauss–Legendre quadrature to the integrals on the ri ght. This will produce the matrix Hnexactly. The parameters in the routine stieltjes sob.m have to be chosen as follows: nd= [n,2n],xw= 2τG 1+ 1 τG 1 2λG 1γλG 1 ............ 2τG n+ 1 τG n 2λG nγλG n τG 1+ 2 λG 1 ...... τG n+ 2 λG n ∈R2n×4, where τG ν,λG νare the nodes and weight of the n-point Gauss–Legendre quadrature rule. Furthermore, a0=1 andsame=1. The complete program, therefore, is as follows: N=6; s=1; a0=1; same=0; g=44000; nd=[N 2*N]; ab=rjacobi(N); zw=gauss(N,ab); xw=zeros(2*N,2*(s+1)); xw(1:N,1)=2*zw(:,1)+1; xw(1:N,2)=zw(:,1); xw(1:N,3)=2*zw(:,2); xw(1:N,4)=g*zw(:,2); xw(N+1:2*N,2)=zw(:,1)+2; xw(N+1:2*N,4)=zw(:,2); B=stieltjes sob(N,s,nd,xw,a0,same); z=sobzeros(N,N,B) It produces the output z = -4.176763898909848e-01 - 1.703657992747233e-01i -4.176763898909848e-01 + 1.703657992747233e-01i 8.453761089539369e-01 - 1.538233952529940e-01i 8.453761089539369e-01 + 1.538233952529940e-01i -1.070135059563751e+00 2.598402134930250e+00 conﬁrming Meijer’s theorem for n= 6. A more detailed numerical study, also in the case of odd values of n, has been made in Ga04, Table 2.30. 36 Exercises to Part I (Stars indicate more advanced exercises.) 1. Explain why, under the assumptions made about the measure dλ, the inner product ( p, q)dλof two polynomials p,qis well deﬁned. 2. Show that monic orthogonal polynomials relative to an abs olutely con- tinuous measure are uniquely deﬁned. {Hint: Use Gram-Schmidt or- thogonalization. }Discuss the uniqueness in the case of discrete mea- sures. 3. Supply the details of the proof of (1.1). In particular, de rive (1.3) and (1.4). 4. Derive the three-term recurrence relation (1.7) for the o rthonormal polynomials. 5. (a) With ˜ πkdenoting the orthonormal polynomials relative to a mea- sure d λ, show that /integraldisplay Rt˜πk(t)˜π/lscript(t)dλ(t) =  0 if |k−/lscript|>1,/radicalbig βk+1if|k−/lscript|= 1, αkifk=/lscript, where αk=αk(dλ),βk=βk(dλ). (b) Use (a) to prove J=Jn(dλ) =/integraldisplay Rtp(t)pT(t)dλ(t), where pT(t) = [˜π0(t),˜π1(t), . . .,˜πn−1(t)]. (c) With notations as in (b), prove tp(t) =Jp(t) +/radicalbig βn˜πn(t)en, where en= [0,0, . . .,1]T∈Rn. 6∗. Symmetry of orthogonal polynomials. Let dλ(t) =w(t)dtbe symmetric on [ −a, a],a >0, i.e. w(−t) =w(t) on [−a, a]. 37 (a) Show that π2k(t; dλ) =π+ k(t2), π 2k+1(t; dλ) =tπ− k(t2), where π± kare the monic polynomials orthogonal on [0 , a2] with respect to d λ±(t) =t∓1/2w(t1/2)dt. Let πk+1(t) =tπk(t)−βkπk−1(t), k= 0,1,2, . . ., π−1(t) = 0, π0(t) = 1 be the recurrence relation for {πk(·; dλ)}, and let α± k,β± kbe the recurrence coeﬃcients for {π± k}. (b Show that β1=α+ 0 β2k=β+ k/β2k−1 β2k+1=α+ k−β2k  k= 1,2,3, . . . . (c) Derive relations similar to those in (b) which involve α+ 0andα− k, β− k. (d) Write a Matlab program that checks the numerical stabili ty of the nonlinear recursions in (b) and (c) when {πk}are the monic Legendre polynomials. 7. The recurrence relation, in Matlab, of the Chebyshev poly nomials of the second kind. (a) Using Matlab, compute Uk(x) for 1 ≤k≤Neither by means of the three-term recurrence relation Un+1(x) = 2xUn(x)−Un−1(x) forn= 0,1, . . ., N −1 (where U−1(x) = 0, U0(x) = 1), or else by putting n= 1 : Nin the explicit formula Un(cosθ) = sin(n+ 1)θ/sinθ, where x= cos θ. For selected values of xand N, determine which of the two methods, by timing each, is more eﬃcient. (b) Using Matlab, compute the single value UN(x) either by use of the three-term recurrence relation, or by direct computation b ased on 38 the trigonometric formula for UN(cosθ). For selected values of x andN, determine which of the two methods, by timing each, is more eﬃcient. 8∗. Orthogonality on two separate (symmetric) intervals. Let 0 < ξ < 1 and consider orthogonal polynomials πkrelative to the weight function w(t) =  |t|γ(1−t2)α(t2−ξ2)β, t ∈[−1, ξ]∪[ξ,1], 0, otherwise . Here, γ∈Randα >−1,β >−1. Evidently, wis a symmetric weight function. Deﬁne π± kas in Exercise 6. (a) Transform the polynomials π± korthogonal on [ ξ2,1] to orthogonal polynomials ˚ π± kon the interval [ −1,1] and obtain the respective weight function ˚ w±. (b) Express β2kandβ2k+1in terms of ˚ γ± r, the leading coeﬃcient of the orthonormal polynomial of degree rrelative to the weight function ˚w±on [−1,1].{Hint: Use βr=/bardblπr/bardbl2//bardblπr−1/bardbl2(cf. eqn (1.4)) and relate this to the leading coeﬃcients γk,γ± k, and ˚γ± k, with obvious notations. } (c) Prove that lim k→∞β2k=1 4(1−ξ)2,lim k→∞β2k+1=1 4(1 +ξ)2. {Hint: Use the result of (b) in combination with the asymptotic equivalence ˚γ± k∼2k˚γ±,˚γ±=π−1/2exp/braceleftbigg −1 2π/integraldisplay1 −1ln ˚w±(x)(1−x2)−1/2dx/bracerightbigg , ask→ ∞ (cf. Sz75, eqn (12.7.2)). You may also want to use /integraldisplay1 0ln(1−a2x2)(1−x2)−1/2dx=πln1 + (1 −a2)1/2 2, a2<1 (see GR00, eqn 4.295.29). } 39 (d) Prove that lim k→∞α± k=1 +ξ2 2,lim k→∞β± k=/parenleftbigg1−ξ2 4/parenrightbigg2 . {Hint: Express α± k,β± kin terms of ˚ α± k,˚β± k, and use the fact that the weight function ˚ w±is in the Szeg¨ o class. } (e) The recurrence coeﬃcients {βk}must satisfy the two nonlinear recursions of Exercise 6(b),(c). Each of them can be interpr eted as a pair of ﬁxed-point iterations for the even-indexed and f or the odd-indexed subsequence, the ﬁxed points being respective ly the limits in (c). Show that, asymptotically, both ﬁxed points a re “at- tractive” for the recursion in 6(b), and “repelling” for the one in 6(c). Also show that in the latter, the ﬁxed points become att rac- tive if they are switched. What are the numerical implicatio ns of all this? (f) Consider the special case γ=±1 and α=β=−1 2. In the case γ= 1, use Matlab to run the nonlinear recursion of Exercise 6(b ) and compare the results with the known answers β2k=1 4(1−ξ)21 +η2k−2 1 +η2k, k= 1,2,3, . . ., 0≤t≤1 and β1=1 2(1 +ξ2), β2k+1=1 4(1 +ξ)21 +η2k+2 1 +η2k, k= 1,2,3, . . ., where η= (1−ξ)/(1 +ξ) (see Ga04, Example 2.30). Likewise, in the case γ=−1, run the nonlinear recursion of Exercise 6(c) and compare the results with the exact answers β2=1 2(1−ξ)2, β2k=1 4(1−ξ)2, k= 2,3, . . ., and β1=ξ, β 2k+1=1 4(1 +ξ)2, k= 1,2,3, . . . . Comment on what you observe. 9. Prove the validity of Algorithm 1. 40 (a) Verify the initialization part. (b) Combine σk+1,k−1= 0 with the three-term recurrence relation for πkto prove the formula for βkin the continuation part. (c) Combine σk+1,k= 0 with the three-term recurrence relation for both, πkandpk, and use the result of (b), to prove the formula forαkin the continuation part. 10∗. Orthogonal polynomials {πk(·;w)}relative to the weight function (“hat function”) w(t) =  1 +tif−1≤t≤0, 1−tif 0≤t≤1, 0 otherwise . (a) Develop a modiﬁed Chebyshev algorithm for generating th e ﬁrst nrecurrence coeﬃcients βk(w),k= 0,1, . . ., n −1 (all αk(w) = 0; why?). Deﬁne modiﬁed moments with respect to a suitable syst em of (monic) orthogonal polynomials. (b) What changes in the routine are required if one wants {πk(·; 1− w)}, or{πk(·;w(1−w))}, or{πk(·;wp)}where p >−1? (c) Download the routine chebyshev.m , write a routine mom.m for the modiﬁed moments to be used in conjunction with chebyshev.m to implement (a), and write a Matlab driver to produce result s for selected values of n. (d) Devise a 2-component discretization scheme for computi ng the ﬁrst nrecurrence coeﬃcients βk(w),k= 0,1,2, . . ., n −1, which uses ann-point discretization of the inner product on each componen t interval and is to yield exact answers (in the absence of roun ding errors). (e) Same as (b). (f) Download the routine mcdis , write a quadrature routine qhatf necessary to implement (d), and append a script to the driver of (c) that produces results of the discretization procedure for s elected values of n. Download whatever additional routines you need. Run the procedure with irout = 1 and irout /negationslash= 1 and observe the respective timings and the maximum discrepancy between the two sets of answers. Verify that the routine “converges” aft er one 41 iteration if idelta is properly set. Compare the results with those of (a). (g) Use the routines acondG.m andrcondG.m to print the absolute and relative condition numbers of the relevant map Gn. Do any of these correlate well with the numerical results obtained in (c)? If not, why not? 11∗. Orthogonal polynomials {πk(·;w)}relative to the weight function (“ex- ponential integral”) w(t) =E1(t), E 1(t) =/integraldisplay∞ 1e−ts sdson [0,∞]. These are of interest in the theory of radiative transfer; se e S. Chan- drasekhar, Radiative transfer, Oxford Univ. Press, Oxford , 1950, Chap- ter II, §23. (a) Develop and run a multiple-component discretization ro utine for generating the ﬁrst nrecurrence coeﬃcients αk(w),βk(w),k= 0,1, . . ., n −1. Check your results for n= 20 against B. Danloy, Math. Comp. 27 (1973), 861–869, Table 3. {Hint: Decompose the interval [0 ,∞] into two subintervals [0 ,2] and [2 ,∞] (addi- tional subdivisions may be necessary to implement the devel op- ments that follow) and incorporate the behavior of E1(t) near t= 0 and t=∞to come up with appropriate discretizations. For 0≤t≤2, use the power series E1(t)−ln(1/t) =−γ−∞/summationdisplay k=1(−1)ktk kk!, where γ=.57721566490153286 . . .is Euler’s constant, and for t >2 the continued fraction (cf. AS92, eqn 5.1.22) tetE1(t) =1 1 +a1 1 +a2 1 +a3 1 +a4 1 +· · ·, a k=⌈k/2⌉/t. Evaluate the continued fraction recursively by (cf. W. Gaut schi, Math. Comp. 31 (1977), 994–999, §2) 1 1 +a1 1 +a2 1 +· · ·=∞/summationdisplay k=0tk, 42 where t0= 1, tk=ρ1ρ2· · ·ρk, k= 1,2,3, . . ., ρ0= 0, ρk=−ak(1 +ρk−1) 1 +ak(1 +ρk−1), k= 1,2,3, . . . . Download the array abjaclog(101:200,:) to obtain the recur- rence coeﬃcients abfor the logarithmic weight function ln(1 /t).} (b) Do the same for w(t) =E2(t), E 2(t) =/integraldisplay∞ 1e−ts s2dson [0,∞]. Check your results against the respective two- and three-po int Gauss quadrature formulae in Chandrasekhar, ibid., Table VI. (c) Do the same for w(t) =Em(t) on [0 , c],0< c < ∞, m= 1,2. Check your results against the respective two-point Gauss q uadra- ture formulae in Chandrasekhar, ibid., Table VII. 12. Let C=b0+a1 b1+a2 b2+a3 b3+· · ·be an inﬁnite continued fraction, and Cn=b0+a1 b1+· · ·an bn=An Bnitsnth convergent. From the theory of continued fractions, it is known that An=bnAn−1+anAn−2 Bn=bnBn−1+anBn2  n= 1,2,3, . . ., where A−1= 1, A0=b0;B−1= 0, B0= 1. Use this to prove (5.2) and (5.3). 13. Prove (5.4). 14. Show that (5.11) implies lim n→∞ρn yn= 0, where ynis any solution of the three-term recurrence relation (satisﬁed by ρnandπn) which is linearly independent of ρn. Thus, {ρn}is indeed a minimal solution. 43 15. Show that the minimal solutions of a three-term recurren ce relation form a one-dimensional manifold. 16. (a) Derive (6.9). (b) Supply the details for deriving Algorithm 4. 17. Supply the details for deriving Algorithm 5. 18. (a) Prove (6.15). (b) Prove (6.17), (6.18). (c) Supply the details for deriving Algorithm 6. 19. Show that a Sobolev inner product does not satisfy the shi ft property (tp, q) = (p, tq). 20. Prove (7.3). 21. The Sobolev inner product (7.1) is called symmetric if each measure dλσis symmetric in the sense of Problem 6. For symmetric Sobolev inner products, (a) show that πk(−t;S) = (−1)kπk(t;S); (b) show that βk 2r= 0 for r= 0,1, . . .,⌊k/2⌋. 44 PART II QUADRATURE 11 Gauss-type quadrature formulae 11.1 Gauss formula Given a positive measure d λ, then-point Gaussian quadrature formula asso- ciated with the measure d λis (11.1)/integraldisplay Rf(t)dλ(t) =n/summationdisplay ν=1λG νf(τG ν) +RG n(f), which has maximum algebraic degree of exactness 2 n−1, (11.2) RG n(f) = 0 if f∈P2n−1. It is well known that the nodes τG νare the zeros of πn(·; dλ), and hence the the eigenvalues of the Jacobi matrix Jn(dλ); cf.§1.2. Interestingly, the weights λG ν, too, can be expressed in terms of spectral data of Jn(dλ); indeed, they are (11.3) λG ν=β0v2 ν,1, where vν,1is the ﬁrst component of the normalized eigenvector vνcorre- sponding to the eigenvalue τG ν, (11.4) Jn(dλ)vν=τG νvν,vT νvν= 1, and, as usual, β0=/integraltext Rdλ(t). This is implemented in the OPQMatlab routine xw=gauss(N,ab) whereab, as in all previous routines, is the N×2 array of recurrence coef- ﬁcients for d λ, andxwtheN×2 array containing the nodes τG νin the ﬁrst column, and the weights λG νin the second. We remark, for later purposes, that the Gauss quadrature sum forf suﬃciently regular can be expressed in matrix form as (11.5)n/summationdisplay ν=1λG νf(τG ν) =β0eT 1f(Jn(dλ))e1,e1= [1,0, . . .,0]T. 45 This is an easy consequence of (11.3) and the spectral decomp osition of Jn, Jn(dλ)V=V D τ,Dτ= diag( τG 1, τG 2, . . ., τG n), where V= [v1,v2, . . .,vn]. Example 6. Zeros of Sobolev orthogonal polynomials of Gegenbauer ty pe (Groenevelt, 2002). The polynomials in question are those orthogonal with respe ct to the Sobolev inner product (u, v)S=/integraldisplay1 −1u(t)v(t)(1−t2)α−1dt+γ/integraldisplay1 −1u/prime(t)v/prime(t)(1−t2)α t2+y2dt. Groenevelt proved that in the case γ→ ∞ the Sobolev orthogonal polyno- mials of even degrees n≥4 have complex zeros if yis suﬃciently small. By symmetry, they must in fact be purely imaginary, and by the re ality of the Sobolev polynomials, must occur in conjugate complex pairs . As we illustrate this theorem, we have an opportunity to apply not only the rou tinegauss.m , but also a number of other routines, speciﬁcally the modiﬁca tion algorithm embodied in the routine chri6.m , dealing with the special quadratic divisor t2+y2in the second integral, and the routine stieltjes sob.m generating the recurrence matrix of the Sobelev orthogonal polynomial s: s=1; same=0; eps0=1e-14; numax=250; nd=[N N]; ab0=rjacobi(numax,alpha); z=complex(0,y); nu0=nu0jac(N,z,eps0); rho0=0; iopt=1; ab1=chri6(N,ab0,y,eps0,nu0,numax,rho0,iopt); zw1=gauss(N,ab1); ab=rjacobi(N,alpha-1); zw=gauss(N,ab); xw=[zw(:,1) zw1(:,1) zw(:,2) gamma*zw1(:,2)]; a0=ab(1,1); B=stieltjes sob(N,s,nd,xw,a0,same); z=sobzeros(N,N,B) Demo#4 The case N=12,α=1 2, and γ= 1 of Example 6. Applying the above routine for y=.1 and y=.09 yields the following zeros (with positive imaginary parts; the other six zeros ar e the same with opposite signs): 46 y zeros y zeros .1 .027543282225 .09 .011086169153 i .284410786673 .281480077515 .541878443180 .540697645595 .756375307278 .755863108617 .909868274113 .909697039063 .989848649239 .989830182743 The numerical results (and additional tests) suggest that G roenevelt’s theo- rem also holds for ﬁnite, not necessarily large, values of γ, and, when γ= 1, that the critical value of ybelow which there are complex zeros must be betweeen .09 and .1. 11.2 Gauss–Radau formula If there is an interval [ a,∞],−∞< a, containing the support of d λ, it may be desirable to have an ( n+ 1)-point quadrature rule of maximum degree of exactness that has aas a prescribed node, (11.6)/integraldisplay Rf(t)dλ(t) =λa 0f(a) +n/summationdisplay ν=1λa νf(τa ν) +Ra n(f). Here, Ra n(f) = 0 for all f∈P2n, andτa νare the zeros of πn(·; dλa), dλa(t) = (t−a)dλ(t). This is called the Gauss–Radau formula . There is again a symmetric, tridiagonal matrix, the Jacobi–Radau matrix (11.7) JR,a n+1(dλ) = Jn(dλ)√βnen √βneT nαR n , αR n=a−βnπn−1(a) πn(a), where en= [0,0, . . .,1]T∈Rn,βn=βn(dλ), and πk(·) =πk(·; dλ), which allows the Gauss–Radau formula to be characterized in terms of eigenvalues and eigenvectors: all nodes of (11.6), including the node a, are the eigen- values of (11.7), and the weights λa νexpressible as in (11.3) in terms of the corresponding normalized eigenvectors vνof (11.7), (11.8) λa ν=β0v2 ν,1, ν= 0,1,2, . . ., n. As in (11.5), this implies that the Gauss–Radau quadrature s um, for smooth f, can be expressed as β0eT 1f(JR,a n)e1. 47 Naturally, if the support of d λis contained in an interval [ −∞, b],b <∞, there is a companion formula to (11.6) which has the prescrib ed node b, (11.9)/integraldisplay Rf(t) dλ(t) =n/summationdisplay ν=1λb νf(τb ν) +λb n+1f(τb n+1) +Rb n(f). The eigenvalue/vector characterization also holds for (11 .9) if in the formula forαR nin (11.7), the variable a, at every occurrence, is replaced by b. The remainder terms of (11.6) and (11.9), if f∈C2n+1[a, b], have the useful property (11.10) Ra n(f)>0, Rb n(f)<0 if sgn f(2n+1)= 1 on [ a, b], with the inequalities reversed if sgn f(2n+1)=−1. For Jacobi resp. generalized Laguerre measures with parame tersα,β resp. α, the quantity αR nis explicitly known (cf. Ga04, Examples 3.4 and 3.5). For example, if a=−1 (in the case of Jacobi measures), (11.11) αR n=−1 +2n(n+α) (2n+α+β)(2n+α+β+ 1), αR n=n, whereas for a= 1, the sign of αR nmust be changed and αandβinterchanged. The respective OPQMatlab routines are xw=radau(N,ab,end0) xw=radau jacobi(N,iopt,a,b) xw=radau laguerre(N,a) In the ﬁrst, abis the (N+1)×2 array of recurrence coeﬃcients for d λ, and end0either a(for (11.6)) or b(for (11.9)). The last two routines make use of the explicit formulae for αR nin the case of Jacobi resp. Laguerre measures, the parameters being α=a,β=b. The parameter ioptchooses between the two Gauss–Radau formulae: the left-handed, if iopt=1, the right-handed otherwise. 11.3 Gauss–Lobatto formula If the support of d λis contained in the ﬁnite interval [ a, b], we may wish to prescribe two nodes, the points aandb. Maximizing the degree of exactness 48 subject to these constraints yields the Gauss–Lobatto formula (11.12)/integraldisplayb af(t)dλ(t) =λL 0f(a) +n/summationdisplay ν=1λL νf(τL ν) +λL n+1f(b) +Ra,b n(f), which we write as an ( n+2)-point formula; we have Ra,b n(f) = 0 for f∈P2n+1. The internal nodes τL νare the zeros of πn(·; dλa,b), dλa,b(t) = (t−a)(b− t)dλ(t). All nodes and weights can be expressed in terms of eigenval ues and eigenvectors exactly as in the two preceding subsections, e xcept that the matrix involved is the Jacobi–Lobatto matrix (11.13) JL n+2(dλ) = Jn+1(dλ)/radicalBig βL n+1en+1/radicalBig βL n+1eT n+1 αL n+1 , where αL n+1andβL n+1are the solution of the 2 ×2 system of linear equations (11.14) πn+1(a)πn(a) πn+1(b)πn(b)  αL n+1 βL n+1 = aπn+1(a) bπn+1(b) . For smooth f, the quadrature sum is expressible as β0eT 1f(JL n+2)e1. For f∈C2n+2[a, b] with constant sign on [ a, b], the remainder Ra,b n(f) satisﬁes (11.15) Ra,b n(f)<0 if sgn f(2n+2)= 1 on [ a, b], with the inequality reversed if sgn f(2n+2)=−1. The parameters αL n+1,βL n+1for Jacobi measures on [ −1,1] with parameters α,βanda=−b=−1 are explicitly known (cf. Ga04, Example 3.8), (11.16)αL n+1=α−β 2n+α+β+ 2, βL n+1= 4(n+α+ 1)(n+β+ 1)(n+α+β+ 1) (2n+α+β+ 1)(2 n+α+β+ 2). TheOPQMatlab routines are xw=lobatto(N,ab,endl,endr) xw=lobatto jacobi(N,a,b) 49 with the meaning of ab,a,bthe same as in the Gauss–Radau routines, and endl=a,endr=b. We remark that both Gauss–Radau and Gauss–Lobatto formulae can be generalized to include boundary points of multiplicity r >1. The inter- nal (simple) nodes and weights are still related to orthogon al polynomials, but the boundary weights require new techniques for their co mputation; see Exercises 8–9. 12 Gauss–Kronrod quadrature In an attempt to estimate the error of the n-point Gauss quadrature rule, Kronrod in 1965 had the idea of inserting n+1 additional nodes and choosing them, along with all 2 n+ 1 weights, in such a way as to achieve maximum degree of exactness. The resulting quadrature rule can be ex pected to yield much higher accuracy than the Gauss formula, so that the diﬀe rence of the two provides an estimate of the error in the Gauss formula. Th e extended formula thus can be written in the form (12.1)/integraldisplay Rf(t)dλ(t) =n/summationdisplay ν=1λK νf(τG ν) +n+1/summationdisplay µ=1λ∗K µf(τK µ) +RGK n(f), and having 3 n+ 2 free parameters λK ν,λ∗K µ,τK µat disposal, one ought to be able to achieve degree of exactness 3 n+ 1, (12.2) RGK n(f) = 0 for f∈P3n+1. A quadrature formula (12.1) that satisﬁes (12.2) is called a Gauss–Kronrod formula . The nodes τK µ, called Kronrod nodes , are the zeros of the polynomial πK n+1of degree n+ 1 which is orthogonal to all polynomials of lower degree in the sense (12.3)/integraldisplay RπK n+1(t)p(t)πn(t; dλ) dλ(t) = 0 for all p∈Pn. Note that the measure of orthogonality here is πn(t; dλ)dλ(t) and thus os- cillates on the support of d λ. Stieltjes (1894) was the ﬁrst to consider poly- nomials πK n+1of this kind (for d λ(t) = d t); a polynomial πK n+1satisfying (12.3) is therefore called a Stieltjes polynomial . Stieltjes conjectured (in the case d λ(t) = d t) that all zeros of πK n+1are real and interlace with the n 50 Gauss nodes— a highly desirable conﬁguration! This has been proved only later by Szeg¨ o (1935) not only for Legendre measures, but al so for a class of Gegenbauer measures. The study of the reality of the zeros fo r more general measures is an interesting and ongoing activity. The computation of Gauss–Kronrod formulae is a challenging problem. A solution has been given only recently by Laurie (1997), at l east in the case when a Gauss–Kronrod formula exists with real nodes and posi tive weights. Interestingly, it can be computed again in terms of eigenval ues and eigen- vectors of a symmetric tridiagonal matrix, just like the pre vious Gauss-type formulae. The relevant matrix, however, is the Jacobi–Kronrod matrix (12.4) JK 2n+1(dλ) = Jn(dλ)√βnen 0 √βneT n αn/radicalbig βn+1eT 1 0/radicalbig βn+1e1 J∗ n . Here, αn=αn(dλ),βn=βn(dλ), etc, and J∗ n(which is partially known) can be computed by Laurie’s algorithm (cf. Ga04, §3.1.2.2). Should some of the eigenvalues of (12.4) turn out to be complex, this would b e an indication that a Gauss–Kronrod formula (with real nodes) does not exis t. There are two routines in OPQ, ab=rkronrod(N,ab0) xw=kronrod(n,ab) that serve to compute Gauss–Kronrod formulae. The ﬁrst gene rates the Jacobi-Kronrod matrix of order 2 N+1, the other the nodes and weights of the quadrature formula, stored respectively in the ﬁrst and second column of the (2N+1)×2 arrayxw. The recurrence coeﬃcients of the given measure dλare input via the ⌈3N/2+1⌉×2 arrayab0. 13 Gauss–Tur´ an quadrature The idea of allowing derivatives to appear in a Gauss-type qu adrature formula is due to Tur´ an (1950). He considered the case where each nod e has the same multiplicity r≥1, that is, (13.1)/integraldisplay Rf(t) dλ(t) =n/summationdisplay ν=1[λνf(τν)+λ/prime νf/prime(τν)+· · ·+λ(r−1) νf(r−1)(τν)]+Rn(f). 51 This is clearly related to Hermite interpolation. Indeed, i f all nodes were prescribed and distinct, one could use Hermite interpolati on to obtain a formula with degree of exactness rn−1 (there are rnfree parameters). Tur´ an asked, like Gauss before him, whether one can do better by cho osing the nodes τνjudiciously. The answer is yes; more precisely, we can get de gree of exactness rn−1 +k,k >0, if and only if (13.2)/integraldisplay Rωr n(t)p(t) dλ(t) = 0 for all p∈Pk−1, where ωn(t) =/producttextn ν=1(t−τν) is the node polynomial of (13.1). We have here a new type of orthogonality: the rth power of ωn, notωn, must be orthogonal to all polynomials of degree k−1. This is called power orthogonality . It is easily seen that rmust be odd, (13.3) r= 2s+ 1, s≥0, so that (13.1) becomes (13.4)/integraldisplay Rf(t) dλ(t) =n/summationdisplay ν=12s/summationdisplay σ=0λ(σ) νf(σ)(τν) +Rn,s(f). Then in (13.2), necessarily k≤n, and k=nis optimal. The maximum possible degree of exactness, therefore, is (2 s+ 2)n−1, and is achieved if (13.5)/integraldisplay R[ωn(t)]2s+1p(t) dλ(t) = 0 for all p∈Pn−1. The polynomial ωn=πn,ssatisfying (13.5) is called s-orthogonal . It exists uniquely and has distinct simple zeros contained in the supp ort interval of dλ. The formula (13.4) is the Gauss-Tur´ an formula if its node polynomial ωn satisﬁes (13.5) and the weights λ(σ) νare obtained by Hermite interpolation. The computation of Gauss-Tur´ an formulae is not as simple as in the case of ordinary Gauss-type formulae. The basic idea, however, i s to consider the positive measure d λn,s(t) = [πn,s(t)]2sdλ(t) and to note that πn,sis the nth-degree polynomial orthogonal relative to d λn,s. The diﬃculty is that this deﬁnes πn,simplicitly, since πn,salready occurs in the measure d λn,s. Nevertheless, the diﬃculty can be surmounted, but at the exp ense of having to solve a system of nonlinear equations; for details, see Ga 04,§3.1.3.2. The procedure is embodied in the OPQroutine 52 xw=turan(n,s,eps0,ab0,hom) where the nodes are stored in the ﬁrst column of the n×(2s+2) array xw, and the successive weights in the remaining 2 s+1 columns. The input parameter eps0is an error tolerance used in the iterative solution of the no nliner sys- tem of equations, and the measure d λis speciﬁed by the (( s+1)n) ×2 input arrayab0of its recurrence coeﬃcients. Finally, hom=1 orhom/negationslash=1 depending on whether or not a certain homotopy in the variable sis used to facilitate convergence of Newton’s method for solving the system of non linear equa- tions. 14 Quadrature formulae based on rational functions All quadrature formulae considered so far were based on poly nomial degree of exactness. This is meaningful if the integrand is indeed p olynomial-like. Not infrequently, however, it happens that the integrand ha s poles outside the interval of integration. In this case, exactness for app ropriate rational functions, in addition to polynomials, is more natural. We d iscuss this for the simplest type of quadrature rule, (14.1)/integraldisplay Rg(t)dλ(t) =n/summationdisplay ν=1λνg(τν) +Rn(g). The problem, more precisely, is to determine λν,τνsuch that Rn(g) = 0 if g∈S2n, where S2nis a space of dimension 2 nconsisting of rational functions and polynomials, (14.2)S2n=Qm⊕P2n−m−1,0≤m≤2n, P2n−m−1= polynomials of degree ≤2n−m−1, Qm= rational functions with prescribed poles . Here, mis an integer of our choosing, and (14.3) Qm= span/braceleftbigg r(t) =1 1 +ζµt, µ= 1,2, . . ., m/bracerightbigg , 53 where (14.4) ζµ∈C, ζ µ/negationslash= 0,1 +ζµt/negationslash= 0 on supp(d λ). The idea is to select the poles −1/ζµof the rational functions in Qmto match the pole(s) of gclosest to the support interval of d λ. In principle, the solution of the problem is rather simple: p utωm(t) =/producttextm µ=1(1 +ζµt) and construct, if possible, the n-point (poynomial) Gauss formula (14.5)/integraldisplay Rg(t)dλ(t) ωm(t)=n/summationdisplay ν=1λG νg(τG ν), g∈P2n−1, for the modiﬁed measure d ˆλ(t) = dλ(t)/ωm(t). Then (14.6) τν=τG ν, λ ν=ωm(τG ν)λG ν, ν= 1,2, . . ., n, are the desired nodes and weights in (14.1). We said “if possible”, since in general ωmis complex-valued, and the existence of a Gauss formula for d ˆλis not guaranteed. There is no problem, however, if ωm≥0 on the support of d λ. Fortunately, in many instances of practical interest, this is indeed the case. There are a number of ways the formula (14.5) can be construct ed: a discretization method using Gauss quadrature relative to d λto do the dis- cretization; repeated application of modiﬁcation algorit hms involving linear or quadratic divisors; special techniques to handle “diﬃcu lt” poles, that is, poles very close to the support interval of d λ. Rather than going into details (which can be found in Ga04, §3.1.4), we present an example taken from solid state physics. Example 7. Generalized Fermi–Dirac integral. This is the integral Fk(η, θ) =/integraldisplay∞ 0tk/radicalbig 1 +θt/2 e−η+t+ 1dt, where η∈R,θ≥0, and kis the Boltzmann constant (=1 2,3 2, or5 2). The ordinary Fermi–Dirac integral corresponds to θ= 0. The integral is conveniently rewritten as (14.7) Fk(η, θ) =/integraldisplay∞ 0/radicalbig 1 +θt/2 e−η+e−tdλ[k](t),dλ[k](t) =tke−tdt, 54 which is of the form (14.1) with g(t) =/radicalbig 1 +θt/2/(e−η+e−t) and d λ= dλ[k] a generalized Laguerre measure. The poles of gevidently are η+µiπ,µ= ±1,±3,±5, . . ., and all are “easy”, that is, at a comfortable distance from the interval [0 ,∞]. It is natural to take meven, and to incorporate the ﬁrst m/2 pairs of conjugate complex poles. An easy computation then yields (14.8) ωm(t) =m/2/productdisplay ν=1[(1 +ξνt)2+ηνt2],2≤m(even) ≤2n, where (14.9) ξν=−η η2+ (2ν−1)2π2, η ν=(2ν−1)π η2+ (2ν−1)2π2. Once the nodes and weights τν,λνhave been obtained according to (14.6), the rational/polynomial quadrature approximation is give n by (14.10) Fk(η, θ)≈N/summationdisplay n=1λn/radicalbig 1 +θτn/2 e−η+e−τn. It is computed in the OPQroutine xw=fermi dirac(N,m,eta,theta,k,eps0,Nmax) whereeps0is an error tolerance, Nmaxa limit on the discretization parameter, and the other variables having obvious meanings. 15 Cauchy principal value integrals When there is a (simple) pole inside the support interval [ a, b] of the measure dλ, the integral must be taken in the sense of a Cauchy principal value integral (15.1) ( Cf)(x; dλ) =/integraldisplayb a−f(t) x−tdλ(t), x∈(a, b). There are two types of quadrature rules for Cauchy principal value integrals: one in which xoccurs as a node, and one in which it does not. They have essentially diﬀerent character and will be considered sepa rately. 55 15.1 Modiﬁed quadrature rule This is a quadrature rule of the form (15.2) ( Cf)(x; dλ) =c0(x)f(x) +n/summationdisplay ν=1cν(x)f(τν) +Rn(f;x). It can be made “Gaussian”, that is, Rn(f;x) = 0 for f∈P2n, by rewriting the integral in (15.1) as (15.3) ( Cf)(x; dλ) =f(x)/integraldisplay R−dλ(t) x−t−/integraldisplay Rf(x)−f(t) x−tdλ(t) and applying the n-point Gauss formula for d λto the second integral. The result is (15.4) ( Cf)(x; dλ) =ρn(x) πn(x)f(x) +n/summationdisplay ν=1λG νf(τG ν) x−τGν+Rn(f;x), where ρn(x) is the Cauchy principal value integral (5.14) and τG ν,λG νare the Gauss nodes and weights. Formula (15.4) is not devoid of numerical diﬃculties. The ma jor one occurs when xapproaches one of the Gauss nodes τG ν, in which case two terms on the right go to inﬁnity, but with opposite signs. In e ﬀect, this means that for xnear a Gaussian node severe cancellation must occur. The problem can be avoided by expanding the integral (15.1) i n Cauchy integrals ρk(x). Let pn(f;·) be the polynomial of degree ninterpolating f at the nGauss nodes τG νand at x. The quadrature sum in (15.4) is then precisely the Cauchy integral of pn, (15.5) ( Cf)(x; dλ) =/integraldisplayb a−pn(f;t) x−tdλ(t) +Rn(f;x). Expanding pnin the orthogonal polynomials πk, (15.6) pn(f;t) =n/summationdisplay k=0akπk(t), a k=1 /bardblπk/bardbl2/integraldisplayb apn(f;t)πk(t)dλ(t), and integrating, one ﬁnds (15.7) ( Cf)(x; dλ) =n/summationdisplay k=0akρk(x) +Rn(f;x), 56 where (15.8) ak=1 /bardblπk/bardbl2n/summationdisplay ν=1λG νf(τG ν)πk(τG ν), k < n ;an=n/summationdisplay ν=1f(x)−f(τG ν) (x−τGν)π/primen(τGν). The Cauchy integrals ρk(x) in (15.7) can be computed in a stable manner by forward recursion; cf. the last paragraph of §5.1. This requires ρ0(x), which is either explicitly known or can be computed by the continue d fraction algorithm; cf. §5.2. Some care must be exercised in computing the divided diﬀerence of fin the formula for an. The procedure is inplemented in the OPQroutine cpvi=cauchyPVI(N,x,f,ddf,iopt,ab,rho0) withiopt/negationslash=1, which produces the ( N+1)-term approximation (15.7) where Rn(f;x) is neglected. The input parameter ddfis a routine for computing the divided diﬀerence of fin a stable manner. It is used only if iopt/negationslash=1. The meaning of the other parameters is obvious. 15.2 Quadrature rule in the strict sense This rule, in which the node t=xis absent, is obtained by interpolating f at the nGauss nodes τG νby a polynomial pn−1(f;·) of degree n−1, f(t) =pn−1(f;t) +En−1(f;t), p n−1(f;t) =n/summationdisplay ν=1πn(t) (t−τG ν)π/prime n(τG ν)f(τG ν), where En−1is the interpolation error, which vanishes identically if f∈Pn−1. The formula to be derived, therefore, will have degree of exa ctness n−1 (which can be shown to be maximum possible). Integrating in t he sense of (15.1) yields (15.9) ( Cf)(x; dλ) =n/summationdisplay ν=1ρn(x)−ρn(τG ν) (x−τG ν)π/prime n(τG ν)f(τG ν) +R∗ n(f;x), where R∗ n(f;x) =/integraltextb a−En−1(f;t)dλ(t)/(x−t). This formula, too, suﬀers from severe cancellation errors w henxis near a Gauss node. The resolution of this problem is similar (in fac t, simpler) than 57 in§15.1: expand pn−1(f;·) in the orthogonal polynomials πkto obtain (15.10)(Cf)(x; dλ) =n−1/summationdisplay k=0a/prime kρk(x) +R∗ n(f;x), a/prime k=1 /bardblπk/bardbl2/integraldisplayb apn−1(f;t)πk(t)dλ(t). It turns out that (15.11) a/prime k=ak, k= 0,1, . . ., n −1, where ak,k < n , is given by (15.8). This is implemented in the OPQroutine cauchyPVI.m withiopt=1. 16 Polynomials orthogonal on several inter- vals We are given a ﬁnite set of intervals [ cj, dj], which may be disjoint or not, and on each interval a positive measure d λj. Let d λbe the “composite” measure (16.1) d λ(t) =/summationdisplay jχ[cj,dj](t)dλj(t), where χ[cj,dj]is the characteristic function of the interval [ cj, dj]. Assuming known the Jacobi matrices J(j)=Jn(dλj) of the component measures d λj, we now consider the problem of determining the Jacobi matrix J=Jn(dλ) of the composed measure d λ. We provide two solutions, one based on Stieltjes’s procedure, and one based on the modiﬁed Chebyshev algorithm . 16.1 Solution by Stieltjes’s procedure The main problem in applying Stieltjes’s procedure is to com pute the inner products ( tπk, πk)dλand (πk, πk)dλfork= 0,1,2, . . ., n −1. This can be done by using Gaussian quadrature on each component interval, (16.2)/integraldisplaydj cjp(t)dλj(t) =n/summationdisplay ν=1λ(j) νp(τ(j) ν), p∈P2n−1. 58 Here we use (11.5) to express the quadrature sum in terms of th e Jacobi matrix J(j), (16.3)/integraldisplaydj cjp(t)dλj(t) =β(j) 0eT 1p(J(j))e1, β(j) 0=/integraldisplaydj cjdλj(t). Then, for the inner products ( tπk, πk)dλ,k≤n−1, we get (tπk, πk)dλ=/integraldisplay Rtπ2 k(t)dλ(t) =/summationdisplay j/integraldisplaydj cjtπ2 k(t)dλj(t) =/summationdisplay jβ(j) 0eT 1J(j)[πk(J(j))]2e1 =/summationdisplay jβ(j) 0eT 1[πk(J(j))]TJ(j)πk(J(j))e1 and for ( πk, πk)dλsimilarly (in fact, simpler) (πk, πk)dλ=/summationdisplay jβ(j) 0eT 1[πk(J(j))]Tπk(J(j))e1. This can be conveniently expressed in terms of the vectors ζ(j) k:=πk(J(j))e1,e1= [1,0, . . .,0]T, which, as required in Stieltjes’s procedure, can be updated by means of the basic three-term recurrence relation. This leads to the fol lowing algorithm Algorithm 7 Stieltjes procedure for polynomials orthogonal on several in- tervals initialization : ζ(j) 0=e1,ζ(j) −1= 0 (all j), α0=/summationtext jβ(j) 0eT 1J(j)e1 /summationtext jβ(j) 0, β0=/summationdisplay jβ(j) 0. continuation (ifn >1): for k= 0,1, . . ., n −2 do ζ(j) k+1= (J(j)−αkI)ζ(j) k−βkζ(j) k−1(allj), αk+1=/summationtext jβ(j) 0ζ(j)T k+1J(j)ζ(j) k+1/summationtext jβ(j) 0ζ(j)T k+1ζ(j) k+1, βk+1=/summationtext jβ(j) 0ζ(j)T k+1ζ(j) k+1/summationtext jβ(j) 0ζ(j)T kζ(j) k. In Matlab, this is implemented in the OPQroutine 59 ab=rmultidomain sti(N,abmd) whereabmd is the array containing the ( α, β)-coeﬃcients of the measures dλ(j). Example 8. Example 1, revisited. This is the case of two identical intervals [ −1,1] and two measures d λ(j) on [−1,1], one a multiple cof the Legendre measure, the other the Cheby- shev measure. This was solved in Example 1 by a 2-component di scretization method. The solution by the 2-domain algorithm of this subse ction, in Mat- lab, looks as follows: ab1=rjacobi(N); ab1(1,2)=2*c; ab2=rjacobi(N,-.5); abmd=[ab1 ab2]; ab=rmultidomain sti(N,abmd) It produces results identical with those produced by the met hod of Example 1. 16.2 Solution by the modiﬁed Chebyshev algorithm The quadrature procedure used in §16.1 to compute inner products can equally be applied to compute the ﬁrst 2 nmodiﬁed moments of d λ, (16.4) mk=/summationdisplay j/integraldisplaydj cjpk(t)dλj(t) =/summationdisplay jβ(j) 0eT 1pk(J(j))e1. The relevant vectors are now z(j) k:=pk(J(j))e1,e1= [1,0, . . .,0]T, and the computation proceeds as in Algorithm 8 Modiﬁed moments for polynomials orthogonal on several in- tervals initialization z(j) 0=e1,z(j) −1= 0 (all j), m 0=/summationdisplay jβ(j) 0. 60 continuation : fork= 0,1, . . .,2n−2 do z(j) k+1= (J(j)−akI)z(j) k−bkz(j) k−1(allj), mk+1=/summationdisplay jβ(j) 0z(j)T k+1e1. With these moments at hand, we can apply Algorithm 1 to obtain the desired recurrence coeﬃcients. This is done in the OPQroutine ab=rmultidomain cheb(N,abmd,abmm) The input array abmdhas the same meaning as in the routine of §16.1, and abmm is a (2N×2) array of the recurrence coeﬃcients ak,bkgenerating the polynomials pk. Applied to Example 8, the Matlab program, using Legendre mom ents (pk the monic Legendre polynomials), is as follows: abm=rjacobi(2*N-1); ab1=abm(1:N,:); ab1(1,2)=2*c; ab2=rjacobi(N,-.5); abmd=[ab1 ab2]; ab=rmultidomain cheb(N,abmd,abm) It produces results identical with those obtained in §16.1, but takes about three times as long to run. 17 Quadrature estimates of matrix function- als The problem to be considered here is to ﬁnd lower and upper bou nds for the quadratic form (17.1) uTf(A)u,u∈RN,/bardblu/bardbl= 1, where A∈RN×Nis a symmetric, positive deﬁnite matrix, fa smooth func- tion (for which f(A) makes sense), and ua given vector. While this looks more like a linear algebra problem, it can actually be solved , for functions f 61 with derivatives of constant sign, by applying Gauss-type q uadrature rules. The connecting link is provided by the spectral resolution ofA, (17.2) AV=VΛ,Λ= diag( λ1, λ2, . . ., λ N),V= [v1,v2, . . .,vN], where λkare the eigenvalues of A(which for simplicity are assumed distinct), andvkthe normalized eigenvectors of A. If we put (17.3) u=N/summationdisplay k=1ρkvk=V ρ,ρ= [ρ1, ρ2, . . ., ρ N]T, and again for simplicity assume ρk/negationslash= 0, all k, then (17.4)uTf(A)u=ρTVTVf(Λ)VTV ρ=ρTf(Λ)ρ, =N/summationdisplay k=1ρ2 kf(λk) =:/integraldisplay R+f(t)dρN(t). This shows how the matrix functional is related to an integra l relative to a discrete positive measure. Now we know from (11.10) and (11. 15) how Gauss– Radau or Gauss–Lobatto rules (and for that matter also ordin ary Gauss rules, in view of RG n= [f(2n)(τ)/(2n)!]/integraltextb a[πn(t; dλ)]2dλ(t),a < τ < b ) can be applied to obtain two-sided bounds for (17.4) when some deri vative of fhas constant sign. To generate these quadrature rules, we need t he orthogonal polynomials for the measure d ρN, and for these the Jacobi matrix JN(dρN). The latter can be computed by an algorithm due to Lanczos (195 0). 17.1 Lanczos algorithm Letρkbe as in (17.4) and h0=/summationtextN k=1ρkvk(=u),/bardblh0/bardbl= 1, as in (17.3). Algorithm 9 Lanczos algorithm initialization : h0prescribed with /bardblh0/bardbl= 1,h−1=0. continuation : forj= 0,1, . . ., N −1 do αj=hT jAhj, ˜hj+1= (A−αjI)hj−γjhj−1, γj+1=/bardbl˜hj+1/bardbl, hj+1=˜hj+1/γj+1. 62 While γ0in Algorithm 9 can be arbitrary (it multiplies h−1=0), it is conve- nient to deﬁne γ0= 1. The vectors h0,h1, . . .,hNgenerated by Algorithm 9 are called Lanczos vectors . It can be shown that αkgenerated by the Lanczos algorithm is precisely αk(dρN), and γk=/radicalbig βk(dρN), fork= 0,1,2, . . ., N −1. This provides us with the Jacobi matrix JN(dρN). 17.2 Examples Example 9. Error bounds for linear algebraic systems. Consider the system (17.5) Ax=b,Asymmetric ,positive deﬁnite . Given an approximation x∗≈x=A−1bto the exact solution x, and the residual vector r=b−Ax∗, we have x−x∗=A−1b+A−1(r−b) =A−1r, thus /bardblx−x∗/bardbl2= (A−1r)TA−1r=rTA−2r, and therefore (17.6) /bardblx−x∗/bardbl2=/bardblr/bardbl2·uTf(A)u, where u=r//bardblr/bardblandf(t) =t−2. All derivatives of fare here of constant sign on R+, (17.7) f(2n)(t)>0, f(2n+1)(t)<0 for t∈R+. By (17.4), we now have (17.8) /bardblx−x∗/bardbl=/bardblr/bardbl2/integraldisplay R+t−2dρN(t). Then-point Gauss quadrature rule applied to the integral on the r ight of (17.8), by the ﬁrst inequality in (17.7), yields a lower bound of/bardblx−x∗/bardbl, without having to know the exact support interval of d ρN. If, on the other hand, we know that the support of d ρNis contained in some interval [ a, b], 0< a < b , we can get a lower bound also from the right-handed ( n+1)-point Gauss–Radau formula, and upper bounds from the left-handed ( n+ 1)-point Gauss–Radau formula on [ a, b], or from the ( n+ 2)-point Gauss–Lobatto formula on [ a, b]. 63 Exercises to Part II (Stars indicate more advanced exercises.) 1. Prove (11.5). 2. Prove that complex zeros of the Sobolev orthogonal polyno mials of Example 6 must be purely imaginary. 3∗. Circle theorems for quadrature weights (cf. P.J. Davis and P. Rabi- nowitz, J. Math. Anal. Appl. 2 (1961), 428–437). (a) Gauss–Jacobi quadrature Letw(t) = (1 −t)α(1 +t)βbe the Jacobi weight function. It is known (Sz75, eqn (15.3.10)) that the nodes τνand weights λνof then-point Gauss–Jacobi quadrature formula satisfy λν∼π nw(τν)/radicalbig 1−τ2ν, n→ ∞, forτνon any compact interval contained in ( −1,1). Thus, suitably normalized weights, plotted against the nodes, lie asympto tically on the unit circle. Use Matlab to demonstrate this graphical ly. (b) Gauss quadrature for the logarithmic weight function w(t) = tαln(1/t) on [0 ,1] (cf. Ga04, Example 2.27). Try, numerically, to ﬁnd a circle theorem in this case also, a nd experiment with diﬀerent values of the parameter α >−1. (Use theOPQroutinerjaclog.m to generate the recurrence coeﬃcients of the orthogonal polynomials for the weight function w.) (c) Gauss–Kronrod quadrature. With was in (a), the analogous result for the 2 n+ 1 nodes τν and weights λνof the (2 n+ 1)-point Gauss–Kronrod formula is expected to be λν∼π 2nw(τν)/radicalbig 1−τ2ν, n→ ∞. That this indeed is the case, when α, β∈[0,5 2), follows from Theorem 2 in F. Peherstorfer and K. Petras, Numer. Math. 95 (2003), 689–706. Use Matlab to illustrate this graphically . (d) Experiment with the Gauss–Kronrod formula for the logar ithmic weight function of (b), when α= 0. 64 4. Discrete orthogonality. Letπk(·; dλ),k= 0,1,2, . . ., be the orthogonal polynomials relative to an absolutely continuous measure. Show that for each N≥2, the ﬁrst Nof them are orthogonal with respect to the discrete inner pro duct (p, q)N=N−1/summationdisplay ν=0λG νp(τG ν)q(τG ν), where τG ν,λG νare the nodes and weights of the N-point Gauss formula for dλ. Moreover, /bardblπk/bardbl2 N=/bardblπk/bardbl2 dλfork≤n−1. 5. (a) Consider the Cauchy integral ρn(z) =ρn(z; dλ) =/integraldisplayb aπn(t; dλ) z−tdλ(t), where [ a, b] is the support of d λ. Show that ρn(z) =O(z−n−1) as z→ ∞. {Hint: Expand the integral deﬁning ρn(z) in descending powers ofz.} (b) Show that /integraldisplayb adλ(t) z−t−σn(z) πn(z)=ρn(z) πn(z)=O(z−2n−1) as z→ ∞. {Hint: Use (5.8). } (c) Consider the partial fraction decomposition σn(z) πn(z)=n/summationdisplay ν=1λν z−τGν ofσn(z)/πn(z) in (5.8). Use (b) to show that λν=λG νare the weights of the n-point Gaussian quadrature formula for d λ. In particular, show that λG ν=σn(τG ν) π/prime n(τG ν). 65 (d) Discuss what happens if z→x,x∈(a, b). 6. Characterize the nodes τb νin (11.9) as zeros of an orthogonal polyno- mial of degree n, and identify the appropriate Gauss–Radau matrix for (11.9). 7. Prove (11.10). {Hint: Use the fact that both formulae (11.6) and (11.9) are interpolatory. } 8. (a) Prove the ﬁrst formula in (11.11). {Hint: Use the relation be- tween the Jacobi polynomials Pk=P(α,β) kcustomarily deﬁned and the monic Jacobi polynomials πk=π(α,β) k, expressed by Pk(t) = 2−k/parenleftbig2k+α+β k/parenrightbig πk(t). You also need Pk(−1) = ( −1)k/parenleftbigk+β k/parenrightbig and the β-coeﬃcient for Jacobi polynomials, βJ n= 4n(n+α)(n+β)(n+ α+β)/(2n+α+β)2(2n+α+β+ 1)(2 n+α+β−1).} (b) Prove the second formula in (11.11). {Hint: With π(α) kandL(α) k denoting the monic resp. conventional generalized Laguerr e poly- nomials, use L(α) k(t) =/parenleftbig (−1)k/k!/parenrightbig π(α) k(t). You also need L(α) k(0) =/parenleftbigk+α k/parenrightbig , and βL n=n(n+α).} 9. Prove (11.16). {Hint: With notation as in the hint to Exercise 6(a), usePk(1) =/parenleftbigk+α k/parenrightbig in addition to the information provided there. } 10. The (left-handed) generalized Gauss–Radau formula is /integraldisplay∞ af(t) dλ(t) =r−1/summationdisplay ρ=0λ(ρ) 0f(ρ)(a) +n/summationdisplay ν=1λR νf(τR ν) +RR n,r(f), where r >1 is the multiplicity of the end point τ0=a, and RR n,r(f) = 0 for f∈P2n−1+r. Let d λ[r](t) = ( t−a)rdλ(t) and τ[r] ν,λ[r],ν= 1,2, . . ., n , be the nodes and weights of the n-point Gauss formula for dλ[r]. (a) Show that τR ν=τ[r] ν, λR ν=λ[r] ν (τR ν−a)r, ν= 1,2, . . ., n. 66 (b) Show hat not only the internal weights λR νare all positive (why?), but also the boundary weights λ0,λ/prime 0ifr= 2. 11. The generalized Gauss–Lobatto formula is /integraldisplayb af(t) dλ(t) =r−1/summationdisplay ρ=0λ(ρ) 0f(ρ)(a)+n/summationdisplay ν=1λL νf(τL ν)+r−1/summationdisplay ρ=0(−1)ρλ(ρ) n+1f(ρ)(b)+RL n,r(f), where r >1 is the multiplicity of the end points τ0=a,τn+1=b, andRL n,r(f) = 0 for P2n−1+2r. Let d λ[r](t) = [(t−a)(b−t)]rdλ(t) and τ[r] ν,λ[r] ν,ν= 1,2, . . ., n , be the nodes and weights of the n-point Gauss formula for d λ[r]. (a) Show that τL ν=τ[r] ν, λL ν=λ[r] ν [(τLν−a)(b−τLν)]r, ν= 1,2, . . ., n. (b) Show that not only the internal weights λL νare all positive (why?), but also the boundary weights λ0,λ/prime 0andλn+1,λ/prime n+1. (c) Show that λ(ρ) 0=λ(ρ) n+1,ρ= 0,1, . . ., r −1, if the measure d λis symmetric. 12∗. Generalized Gauss-Radau quadrature. (a) Write a Matlab routine gradau.m for generating the generalized Gauss-Radau quadrature rule of Exercise 8 for a measure d λon [a,∞], having a ﬁxed node aof multiplicity r,r >1.{Hint: To compute the boundary weights, set up an (upper triangular ) system of linear equations by applying the formula in turn wi th π2 n(t), (t−a)π2 n(t), . . .,(t−a)r−1π2 n(t), where πn(t) =/producttextn ν=1(t− τR ν).} (b) Check your routine against the known formulae with r= 2 for the Legendre and Chebyshev measures (see Ga04, Examples 3.1 0 and 3.11). Devise and implement a check that works for arbitr ary r≥2 and, in particular, for r= 1. (c) Use your routine to explore positivity of the boundary we ights and see whether you can come up with any conjectures. 67 13∗. Generalized Gauss-Lobatto quadrature. (a) Write a Matlab routine globatto.m for generating the generalized Gauss-Lobatto rule of Exercise 9 for a measure d λon [a, b], having ﬁxed nodes at aandbof multiplicity r,r >1. For simplicity, start with the case r≥2 even; then indicate the changes necessary to deal with odd values of r.{Hint: Similar to the hint in Exercise 10(a). } (b) Check your routine against the known formulae with r= 2 for the Legendre and Chebyshev measures (see Ga04, Exercises 3. 13 and 3.14). Devise and implement a check that works for arbitr ary r≥2 and, in particular, for r= 1. (c) Explore the positivity of the boundary weights λ(ρ) 0and the quan- titiesλ(ρ) n+1in the quadrature formula. 14. Show that the monic Stieltjes polynomial πK n+1in (12.3) exists uniquely. 15. (a) Let d λbe a positive measure. Use approximation theory to show that the minimum of/integraltext R|π(t)|pdλ(t), 1< p < ∞, extended over all monic polynomials πof degree nis uniquely determined. (b) Show that the minimizer of the extremal problem in (a), wh en p= 2s+ 2,s≥0 an integer, is the s-orthogonal polynomial π=πn,s.{Hint: Diﬀerentiate the integral partially with respect to the variable coeﬃcients of π.} 16. (a) Show that rin (13.2) has to be odd. (b) Show that in (13.2) with ras in (13.3), one cannot have k > n. 17. Derive (14.8) and (14.9). 18. Derive (15.4) from (15.3) and explain he meaning of Rn(f;x).{Hint: Use Exercise 5(c) and (5.8). } 19. Prove (15.5). {Hint: Use Exercise 5(c). } 20. Prove (15.8). {Hint: For k < n , use Gauss quadrature, and for k=n insert the expression for pn(f;t) from (15.4) into the formula for anin (15.6). Also use the fact that the elementary Lagrange inter polation polynomials sum up to 1. } 68 21. Derive (15.9). 22. Prove (15.11). {Hint: Use Exercise 5(a). } 23. Derive (15.9). 24. Prove (15.11). 25. (a) Prove that the Lanczos vectors are mutually orthonor mal. (b) Show that the vectors {hj}n j=0,n < N , form an orthonormal basis of the Krylov space Kn(A,h0) = span( h0,Ah0, . . .,Anh0). (c) Prove that hj=pj(A)h0, j= 0,1, . . ., N, where pjis a polynomial of degree jsatisfying the three-term recurrence relation γj+1pj+1(λ) = (λ−αj)pj(λ)−γjpj−1(λ), j= 0,1, . . ., N −1, p0(λ) = 1, p−1(λ) = 0. {Hint: Use mathematical induction. } 26. Prove that the polynomial pkof Exercise 25(c) is equal to the orthonor- mal polynomial ˜ πk(·; dρN).{Hint: Use the spectral resolution of Aand Exercises 25(a) and (c). } 69 PART III APPROXIMATION 18 Polynomial least squares approximation 18.1 Classical least squares problem We are given Ndata points ( tk, fk),k= 1,2, . . ., N , and wish to ﬁnd a polynomial ˆ πnof degree ≤n,n < N , such that a weighted average of the squared errors [ p(tk)−fk]2is as small as possible among all polynomials p of degree n, (18.1)N/summationdisplay k=1wk[ˆpn(tk)−fk]2≤N/summationdisplay k=1wk[p(tk)−fk]2for all p∈Pn. Here, wk>0 are positive weights, which allow placing more emphasis on data points that are reliable, and less emphasis on others, b y choosing them larger resp. smaller. If the quality of the data is uniformly the same, then equal weights, say wk= 1, are appropriate. The problem as formulated suggests a discrete N-point measure (18.2) d λN(t) =N/summationdisplay k=1wkδ(t−tk), δ= Dirac delta function , in terms of which the problem can be written in the compact for m (18.3) /bardblˆpn−f/bardbl2 dλN≤ /bardblp−f/bardbl2 dλNfor all p∈Pn. The polynomials πk(·) =πk(·; dλN) orthogonal (not necessarily monic) with respect to the discrete measure (18.2) provide an easy solut ion: one writes (18.4) p(t) =n/summationdisplay i=0ciπi(t), n < N, and obtains for the squared error, using the orthogonality o fπk, (18.5) E2 n=/parenleftBiggn/summationdisplay i=0ciπi−f,n/summationdisplay j=0cjπj−f/parenrightBigg =n/summationdisplay i,j=0cicj(πi, πj)−2n/summationdisplay i=0ci(f, πi) +/bardblf/bardbl2 =n/summationdisplay i=0/parenleftbigg /bardblπi/bardblci−(f, πi) /bardblπi/bardbl/parenrightbigg2 +/bardblf/bardbl2−n/summationdisplay i=0(f, πi)2 /bardblπi/bardbl2. 70 (All norms and inner products are understood to be relative t o the measure dλN.) Evidently, the minimum is attained for ci= ˆci(f), where (18.6) ˆ ci(f) =(f, πi) /bardblπi/bardbl2, i= 0,1, . . ., n, are the “Fourier coeﬃcients” of frelative to the orthogonal system π0, π1, . . ., π N−1. Thus, (18.7) ˆ pn(t) =n/summationdisplay i=0ˆci(f)πi(t; dλN). In Matlab, the procedure is implemented in the OPQroutine [phat,c]=least squares(n,f,xw,ab,d) The given function values fkare input through the N×1 arrayf, the abscissae tkand weights wkthrough the N×2 arrayxw, and the measure d λNthrough the (N+1)×2 arrayabof recurrence coeﬃcients. The 1 ×(n+1) arraydis the vector of leading coeﬃcients of the orthogonal polynomials . The procedure returns as output the N×(n+1) arrayphatof the values ˆ pν(tk), 0≤ν≤n, 1≤k≤N, and the ( n+1)×1 arraycof the Fourier coeﬃcients. Example 10. Equally weighted least squares approximation on N= 10 equally spaced points on [ −1,1]. Matlab program: N=10; k=(1:N)’; d=ones(1,N); xw(k,1)=-1+2*(k-1)/(N-1); xw(:,2)=2/N; ab=rhahn(N-1); ab(:,1)=-1+2*ab(:,1)/(N-1); ab(:,2)=(2/(N-1))^2*ab(:,2); ab(1,2)=2; [phat,c]=least squares(N-1,f,xw,ab,d); Demo#5 The program is applied to the function f(t) = ln(2+ t) on [−1,1], and selected least squares errors ˆEnare compared in the table below with maximum errors E∞ n(taken over 100 equally spaced points on [ −1,1]). n ˆEn E∞ n 0 4.88(–01) 6.37(–01) 3 2.96(–03) 3.49(–03) 6 2.07(–05) 7.06(–05) 9 1.74(–16) 3.44(–06) 71 Ifn=N−1, the least squares error ˆEN−1is zero, since the Ndata points can be interpolated exactly by a polynomial of degree ≤N−1. This is conﬁrmed in the ﬁrst tabular entry for n= 9. 18.2 Constrained least squares approximation It is sometimes desirable to impose constraints on the least squares approx- imation, for example to insist that at certain points sjthe error should be exactly zero. Thus, the polynomial p∈Pnis subject to the constraints (18.8) p(sj) =fj, j= 1,2, . . ., m ;m≤n, but otherwise is freely variable. For simplicity we assume t hat none of the sjequals one of the support points tk. (Otherwise, the procedure to be described requires some simple modiﬁcations.) In order to solve the constrained least squares problem, let (18.9) pm(f;t) =pm(f;s1, . . ., s m;t), σ m(t) =m/productdisplay j=1(t−sj) be respectively the polynomial of degree m−1 interpolating fat the points sjand the constraint polynomial of degree m. We then write (18.10) p(t) =pm(f;t) +σm(t)q(t). This clearly satisﬁes the constraints (18.8), and qis a polynomial of degree n−mthat can be freely varied. The problem is to minimize the squa red error /bardblf−pm(f;·)−σmq/bardbl2 dλN=/integraldisplay R/bracketleftbiggf(t)−pm(f;t) σm(t)−q(t)/bracketrightbigg2 σ2 m(t)dλN over all polynomials qof degree n−m. This is an unconstrained least squares problem, but for a new function f∗and a new measure d λ∗ N, (18.11) minimize : /bardblf∗−q/bardbldλ∗ N, q∈Pn−m, where (18.12) f∗(t) =f(t)−pm(f;t) σm(t),dλ∗ N(t) =σ2 m(t)dλN(t). 72 If ˆqn−mis the solution of (18.11), then (18.13) ˆ pn(t) =pm(f;t) +σm(t)ˆqn−m(t) is the solution of the constrained least squares problem. Th e function f∗, incidentally, can be given the form of a divided diﬀerence, f∗(t) = [s1, s2, . . ., s m, t]f, t ∈supp d λ∗ N, as follows from the theory of interpolation. Note also that t he discrete or- thogonal polynomials πk(·; dλ∗ N) needed to solve (18.11) can be obtained from the polynomials πk(·; dλN) bymmodiﬁcations of the measure d λNby linear factors sj. Example 11. Bessel function J0(t) for 0 ≤t≤j0,3. Here, j0,3is the third positive zero of J0. A natural constraint is to reproduce the ﬁrst three zeros of J0exactly, that is, m= 3 and s1=j0,1, s2=j0,2, s3=j0.3. Demo#6 The constrained least squares approximations of degrees n= 3,4,5 (that is, n−m= 0,1,2) using N= 51 equally spaced points on [0, j0,3] (end points included) are shown in the ﬁgure below. The soli d curve 0 1 2 3 4 5 6 7 8 9−0.500.511.5 xBessel represents the exact function, the dashdotted, dashed, an d otted curves the approximants for n= 3,4, and 5, respectively. The approximations are not particularly satisfactory and show spurious behavior near t= 0. 73 Example 12. Same as Example 11, but with two additional constraints p(0) = 1 , p/prime(0) = 0 . Demo#7 Derivative constraints, as the one in Example 12, can be inco rpo- rated similarly as before. In this example, the added constr aints are designed to remove the spurious behavior near t= 0; they also improve considerably the overall accuracy, as is shown in the next ﬁgure. 0 1 2 3 4 5 6 7 8 9−0.500.51 xBessel 18.3 Least squares approximation in Sobolev spaces The task now is to approximate simultaneously functions and some of their ﬁrst derivatives. More precisely, we want to minimize s/summationdisplay σ=0N/summationdisplay k=1w(σ) k[p(σ)(tk)−f(σ) k]2 over all polynomials p∈Pn, where f(σ) k,σ= 0,1, . . ., s , are given function and derivative values, and w(σ) k>0 appropriate weights for each derivative. These are often chosen such that w(σ) k=γσwk, γσ>0, k= 1,2, . . ., N, in terms of one set of positive weights wk. Evidently,, the problem, analo- gously to (18.3), can be written in terms of the Sobolev inner product and 74 norm (18.14) ( u, v)S=s/summationdisplay σ=0N/summationdisplay k=1w(σ) ku(σ)(tk)v(σ)(tk),/bardblu/bardblS=/radicalbig (u, u)S in the compact form (18.15) minimize : /bardblp−f/bardbl2 Sfor all p∈Pn. The solution is entirely analogous to the one provided in §18.1, (18.16) ˆ pn(t) =n/summationdisplay i=0ˆci(f)πi(t),ˆci(f) =(f, πi)S /bardblπi/bardbl2 S, where {πi}are the orthogonal polynomials of Sobolev type. In Matlab, t he procedure is [phat,c]=least squares sob(n,f,xw,B) The input parameter fis now an N×(s+ 1) array containing the Nvalues of the given function and its ﬁrst sderivatives at the points tk. The abscissae tk and the weights w(σ) kof the Sobolev inner product are input via the N×(s+1) arrayxw. (The routine determines sautomatically from the size of the array xw.) The user also has to provide the N×Nupper triangular array of the recurrence coeﬃcients for the Sobolev orthogonal polynomi als, which for s= 1 can be generated by the routine chebyshev sob.m and for arbitrary sby the routine stieltjes sob.m . The output phat is an array of dimension (n+1)×(Ns) containing the Nvalues of the derivative of order σof thenth degree approximant ˆ pnin positions ( n+1,σ:s+1:Ns) of the array phat. The Fourier coeﬃcients ˆ ciare output in the ( n+1)×1 vectorc. Example 12. The complementary error function on [0 ,2]. This is the function f(t) =et2erfct=2√πet2/integraldisplay∞ te−u2du,0≤t≤2, whose derivatives are easily calculated. Demo#8 The routine leastsquares sob.m is applied to the function fof Example 12 with s= 2 andN=5 equally spaced points tkon [0,2]. All weights are chosen to be equal, w(σ) k= 1/Nforσ= 0,1,2. The table below, in the top half, shows selected results for the Sobolev least squares error ˆEn 75 s n ˆEn E∞ n,0 E∞ n,1 E∞ n,2 2 0 1.153(+00) 4.759(–01) 1.128(+00) 2.000(+00) 2 7.356(–01) 8.812(–02) 2.860(–01) 1.411(+00) 4 1.196(–01) 1.810(–02) 5.434(–02) 1.960(–01) 9 2.178(–05) 4.710(–06) 3.011(–05) 3.159(–04) 14 3.653(–16) 1.130(–09) 1.111(–08) 1.966(–07) 0 0 2.674(–01) 4.759(–01) 1.128(+00) 2.000(+00) 2 2.245(–02) 3.865(–02) 3.612(–01) 1.590(+00) 4 1.053(–16) 3.516(–03) 5.160(–02) 4.956(–01) 9 1.053(–16) 5.409(–03) 8.124(–02) 7.959(–01) 14 1.053(–16) 5.478(–03) 8.226(–02) 8.057(–01) and the maximum errors E∞ n,0,E∞ n,1,E∞ n,2(over 100 equally spaced points on [0,2]) for the function and its ﬁrst two derivatives. In the bott om half are shown the analogous results for ordinary least squares a pproximation (s= 0). Note that the Sobolev least squares error ˆE3N−1is essentially zero, reﬂecting the fact that the Hermite interpolation polynomi al of degree 3 N−1 interpolates the data exactly. In contrast, ˆEn= 0 for n≥N−1 in the case of ordinary least squares. As expected, the table shows rather convincingly that Sobol ev least squares approximation approximates the derivatives decidedly bet ter than ordinary least squares approximation, and even the function itself w hennis suﬃciently large. 19 Moment-preserving spline approximation There are various types of approximation: those that contro l the maximum pointwise error; those that control some average error (lik e least squares error); and those, often motivated by physical considerati ons, that try to preserve the moments of the given function, or at least as man y of the ﬁrst moments as possible. It is this last type of approximation th at we now wish to study. We begin with piecewise constant approximation on the whole real lineR+, then proceed to spline approximation on R+, and end with spline approximation on a compact interval. 76 19.1 Piecewise constant approximation on R+ The piecewise constant approximants to be considered are (19.1) sn(t) =n/summationdisplay ν=1aνH(tν−t), t∈R+, where aν∈R, 0< t1< t2<· · ·< tn, and His the Heaviside function H(u) =  1 ifu≥0, 0 otherwise . The problem is, for given f∈C1(R+), to ﬁnd, if possible, the aνandtνsuch that (19.2)/integraldisplay∞ 0sn(t)tjdt=µj, j= 0,1, . . .,2n−1, where (19.3) µj=/integraldisplay∞ 0f(t)tjdt, j = 0,1, . . .,2n−1, are the moments off, assumed to exist. The solution can be formulated in terms of Gauss quadrature r elative to the measure (19.4) d λ(t) =−tf/prime(t)dtonR+. Indeed, if f(t) =o(t−2n) ast→ ∞, then the problem has a unique solution if and only if d λin (19.4) admits an n-point Gauss quadrature formula (19.5)/integraldisplay∞ 0g(t)dλ(t) =n/summationdisplay ν=1λG νg(τG ν), g∈P2n−1, satisfying 0 < τG 1< τG 2<· · ·< τG n. If that is the case, then the desired knots tνand coeﬃcients aνare given by (19.6) tν=τG ν, aν=λG ν τG ν, ν= 1,2, . . ., n. 77 A Gauss formula (19.5) always exists if f/prime<0 onR+, that is, d λ(t)≥0. For the proof, we use integration by parts, /integraldisplayT 0f(t)tjdt=1 j+ 1tj+1f(t)/vextendsingle/vextendsingle/vextendsingle/vextendsingleT 0−1 j+ 1/integraldisplayT 0f/prime(t)tj+1dt, j ≤2n−1, and let T→ ∞. The integrated part on the right goes to zero by assumption onf, and the left-hand side converges to the jth moment of f, again by assumption. Therefore, the last term on the right also conve rges, and since −tf/prime(t) = dλ(t), one ﬁnds µj=1 j+ 1/integraldisplay∞ 0tjdλ(t), j= 0,1, . . .,2n−1. This shows in particular that the ﬁrst 2 nmoments of d λexist, and therefore, if dλ≥0, also the Gauss formula (19.5). On the other hand, the approximant snhas moments /integraldisplay∞ 0sn(t)tjdt=n/summationdisplay ν=1aν/integraldisplayτν 0tjdt=1 j+ 1n/summationdisplay ν=1aνtj+1 ν, so that the ﬁrst 2 nmoments µjoffare preserved if and only if n/summationdisplay ν=1(aνtν)tj ν=/integraldisplay∞ 0tjdλ(t), j= 0,1, . . .,2n−1. This is equivalent to saying that the knots tνare the Gauss nodes in (19.5), andaνtνthe corresponding weights. Example 13. Maxwell distribution f(t) = e−t2onR+. Here, dλ(t) = 2t2e−t2dtonR+, which is a positive measure obtained (up to the factor 2) by tw ice modi- fying the half-range Hermite measure by a linear factor t. The ﬁrst n+ 2 recurrence coeﬃcients of the half-range Hermite measure ca n be computed by a discretization method. Applying to these recurrence co eﬃcients twice the routine chri1.m , with zero shift, then yields the recurrence coeﬃcients αk(dλ),βk(dλ),k≤n−1, and hence the required n-point Gauss quadrature rule (19.5) for d λ. The result for n= 5 is depicted in the ﬁgure below. 78 0 0.5 1 1.5 2 2.5 300.20.40.60.81 t1t2t3t4t5 19.2 Spline approximation on R+ The approximant snof§19.1 can be interpreted as a spline function of degree 0. We now consider spline functions sn,mof degree m >0, (19.7) sn,m(t) =n/summationdisplay ν=1aν(tν−t)m +, t∈R+, where um +is the truncated power um +=umifu≥0, and um += 0 if u <0. Given the ﬁrst 2 nmoments (19.3) of f, the problem again is to determine aν∈Rand 0 < t1< t2<· · ·< tnsuch that (19.8)/integraldisplay∞ 0sn,m(t)tjdt=µj, j= 0,1, . . .,2n−1. By a reasoning similar to the one in §19.1, but more complicated, involving mintegrations by part, one proves that for f∈Cm+1(R+) and satisfying f(µ)(t) =o(t−2n−µ) ast→ ∞,µ= 0,1, . . ., m , the problem has a unique solution if and only if the measure (19.9) d λ[m](t) =(−1)m+1 m!tm+1f(m+1)(t)dtonR+ admits an n-point Gauss quadrature formula (19.10)/integraldisplay∞ 0g(t) dλ[m](t) =n/summationdisplay ν=1λG νg(τG ν) for all g∈P2n−1 79 satisfying 0 < τG 1< τG 2<· · ·< τG n. If that is the case, the knots tνand coeﬃcients aνare given by (19.11) tν=τG ν, aν=λG ν [τG ν]m+1, ν= 1,2, . . ., n. Note that d λ[m]in (19.9) is a positive measure, for each m≥0, and hence (19.10) exists, if fis completely monotonic on R+, that is, ( −1)µf(µ)(t)>0, t∈R+, forµ= 0,1,2, . . .. Example 14. Maxwell distribution f(t) = e−t2onR+, revisited. We now have dλ[m](t) =1 m!tm+1Hm+1(t)e−t2dtonR+, where Hm+1is the Hermite polynomial of degree m+1. Here, d λ[m]ifm >0 is no longer of constant sign on R+, and hence the existence of the Gauss rule (19.10) is in doubt. Numerical exploration, using discreti zation methods, yields the situation shown in the table below, where a dash in dicates the presence of a negative Gauss node τG ν, and an asterisk the presence of a pair n m = 1 m= 2 m= 3 n m = 1 m= 2 m= 3 1 6.9(–2) 1.8(–1) 2.6(–1) 11 — 1.1(–3) 1.1(–4) 2 8.2(–2) — 2.3(–1) 12 — — * 3 — 1.1(–2) 2.5(–3) 13 7.8(–3) 6.7(–4) * 4 3.5(–2) 6.7(–3) 2.2(–3) 14 8.3(–3) 5.6(–4) 8.1(–5) 5 2.6(–2) — 1.6(–3) 15 7.7(–3) — 7.1(–5) 6 2.2(–2) 3.1(–3) * 16 — 4.9(–4) 7.8(–5) 7 — 2.4(–3) * 17 — 3.8(–4) 3.8(–5) 8 1.4(–2) — 3.4(–4) 18 5.5(–3) 3.8(–4) * 9 1.1(–2) 1.7(–3) 2.5(–4) 19 5.3(–3) — * 10 9.0(–3) 1.1(–3) — 20 5.4(–3) 3.1(–4) * of conjugate complex Gauss nodes. In all cases computed, the re were never more than one negative Gauss node, or more than one pair of com plex nodes. The numbers in the table represent the maximum errors /bardblsn,m−f/bardbl∞, the maximum being taken over 100 equally spaced points on [0 , τG n]. 80 19.3 Spline approximation on a compact interval The problem on a compact interval, say [0 ,1], is a bit more involved than the problem on R+. For one, the spline function sn,mmay now include a polynomial pof degree m, which was absent before since no moment of p exists on R+unless p≡0. Thus, the spline approximant has now the form (19.12) sn,m(t) =p(t) +n/summationdisplay ν=1aν(tν−t)m +, p∈Pm,0≤t≤1, where aν∈Rand 0 < t1< t2<· · ·< tn<1. There are two problems of interest: Problem I. Find sn,msuch that (19.13)/integraldisplay1 0sn,m(t)tjdt=µj, j= 0,1, . . .,2n+m. Since we have m+ 1 additional parameters at our disposal (the coeﬃcients ofp), we can impose m+ 1 additionl moment conditions. Problem II. Rather than matching more moments, we use the added degre e of freedom to impose m+ 1 “boundary conditions” at the end point t= 1. More precisely, we want to ﬁnd sn,msuch that (19.14)/integraldisplay1 0sn,m(t)tjdt=µj, j= 0,1, . . .,2n−1 and (19.15) s(µ) n,m(1) = f(µ)(1), µ= 0,1, . . ., m. It is still true that a solution can be given in terms of quadra ture formulae, but they are now respectively generalized Gauss–Lobatto an d generalized Gauss–Radau formulae relative to the measure (19.16) d λ[m](t) =(−1)m+1 m!f(m+1)(t)dton [0,1] (see M. Frontini, W. Gautschi, and G.V. Milovanovi´ c, Facta Univ. Ser. Math. Inform. 4 (1989), 45–56). Problem I, in fact, has a uniq ue solution if 81 and only if the generalized Gauss–Lobatto formula (19.17)/integraldisplay1 0g(t)dλ[m](t) =m/summationdisplay µ=0[λ(µ) 0g(µ)(0) + ( −1)µλ(µ) n+1g(µ)(1)] +n/summationdisplay ν=1λL νg(τL ν), g∈P2n+2m+1, exists with 0 < τL 1<· · ·< τL n<1. In this case, (19.18) tν=τL ν, aν=λL ν, ν= 1,2, . . ..n, andpis uniquely determined by (19.19) p(µ)(1) = f(µ)(1) + ( −1)mm!λ(m−µ) n+1, µ= 0,1, . . ., m. Similarly, Problem II has a unique solution if and only if the generalized Gauu–Radau formula (19.20)/integraldisplay1 0g(t)dλ[m](t) =m/summationdisplay µ=0λ(µ) 0g(µ)(0) +n/summationdisplay ν=1λR νg(τR ν), g∈P2n+m, exists with 0 < τR 1<· · ·< τR n<1. Then (19.21) tν=τR ν, aν=λR ν, ν= 1,2, . . ..n, and (trivially) (19.22) p(t) =m/summationdisplay µ=0f(µ)(1) µ!(t−1)µ. In both cases, complete monotonicity of fimplies d λ≥0 and the existence of the respective quadrature formulae. 20 Slowly convergent series Standard techniques of accelerating the convergence of slo wly convergent series are based on linear or nonlinear sequence transforma tions: the sequence of partial sums is transformed somehow into a new sequence th at converges 82 to the same limit, but a lot faster. Here we folllow another ap proach, more in the spirit of these lectures: the sum of the series is repre sented as a deﬁnite integral; a sequence of quadrature rules is then app lied to this integral which, when properly chosen, will produce a sequence of appr oximations that converges quickly to the desired sum. An easy way (and certainly not the only one) to obtain an inegr al repre- sentation presents itself when the general term of the serie s, or part thereof, is expressible in terms of the Laplace transform of a known fu nction. Several instances of this will now be described. 20.1 Series generated by a Laplace transform The series (20.1) S=∞/summationdisplay k=1ak to be considered ﬁrst has terms akthat are the Laplace transform (Lf)(s) =/integraldisplay∞ 0e−stf(t)dt of some known function fevaluated at s=k, (20.2) ak= (Lf)(k), k= 1,2,3, . . . . In this case, S=∞/summationdisplay k=1/integraldisplay∞ 0e−ktf(t)dt =/integraldisplay∞ 0∞/summationdisplay k=1e−(k−1)t·e−tf(t)dt =/integraldisplay∞ 01 1−e−te−tf(t)dt that is, (20.3) S=/integraldisplay∞ 0t 1−e−tf(t) te−tdt. There are at least three diﬀerent approaches to evaluate thi s integral nu- merically: one is Gauss–Laguerre quadrature of ( t/(1−e−t))f(t)/twith 83 dλ(t) = e−tdtonR+; another is rational/polynomial Gauss–Laguerre quadra- ture of the same function; and a third Gauss–Einstein quadra ture of the function f(t)/twith d λ(t) =tdt/(et−1) on R+. In the last method, the weight function t/(et−1) is widely used in solid state physics, where it is named after Einstein (coming from the Einstein-Bose distri bution). It is also, incidentally, the generating function of the Bernoulli pol ynomials. Example 15. The Theodorus constant S=∞/summationdisplay k=11 k3/2+k1/2= 1.860025 . . . . This is a universal constant introduced by P.J. Davis (1993) in connection with a spiral attributed to the ancient mathematician Theod orus of Cyrene. Here we note that 1 s3/2+s1/2=s−1/21 s+ 1=/parenleftbigg L1√ πt∗e−t/parenrightbigg (s), where the star stands for convolution. A simple computation yields (20.2) with f(t) =2√πF(√ t), where F(x) = e−x2/integraldisplayx 0et2dt is Dawson’s integral. Demo#9 To make f(t) regular at t= 0, we divide by√ tand write S=2√π/integraldisplay∞ 0t 1−e−tF(√ t)√ tt−1/2e−tdt =2√π/integraldisplay∞ 0F(√ t)√ tt−1/2t et−1dt. To the ﬁrst integral we apply Gauss–Laguerre quadrature wit h dλ(t) = t−1/2e−tdtonR+, or rational Gauss–Laguerre with the same d λ, and to the second integral Gauss–Einstein quadrature (modiﬁed by the factor t−1/2). The errors committed in these quadrature methods are shown i n the table below. 84 nGauss-Laguerre rational Gauss-Laguerre Gauss-Einstein 1 9.6799(–03) 1.5635(–02) 1.3610(–01) 4 5.5952(–06) 1.1893(–08) 2.1735(–04) 7 4.0004(–08) 5.9689(–16) 3.3459(–07) 10 5.9256(–10) 5.0254(–10) 15 8.2683(–12) 9.4308(–15) 20 8.9175(–14) 4.7751(–16) timing: 10.8 timing: 8.78 timing: 10.4 The clear winner is rational Gauss–Laguerre, both in terms o f accuracy and run time. Example 16. The Hardy–Littlewood function H(x) =∞/summationdisplay k=11 ksinx k, x > 0. It can be shown that ak:=1 ksinx k= (Lf(t;x))(k), where f(t;x) =1 2i[I0(2√ ixt)−I0(2√ −ixt)] andI0is the modiﬁed Bessel function. This gives rise to the two int egral representations H(x) =/integraldisplay∞ 0t 1−e−tf(t;x) te−tdt=/integraldisplay∞ 0f(t;x) tt et−1dt. Among the three quadrature methods, Gauss–Einstein perfor ms best, but all suﬀer from internal cancellation of terms in the quadrat ure sum. The problem becomes more prominent as the number nof terms increases. The ﬁgure below shows the behavior of H(x) in the range 0 ≤x≤100. 20.2 “Alternating” series generated by a Laplace trans- form These are series in which the general terms are Laplace trans forms with alternating signs of some function f, that is, series (20.1) with (20.4) ak= (−1)k−1(Lf)(k), k= 1,2,3, . . . . 85 0 10 20 30 40 50 60 70 80 90 100−0.500.511.522.533.5 xH(x) An elementary computation similar to the one carried out in §20.1 will show that (20.5) S=/integraldisplay∞ 01 1 + e−tf(t)e−tdt=/integraldisplay∞ 0f(t)1 et+ 1dt. We can again choose between three quadrature methods: Gauss –Laguerre quadrature of the function f(t)/(1+e−t) with d λ(t) = e−tdt, rational/polynomial Gauss–Laguerre of the same function, and Gauss-Fermi quadr ature of f(t) with d λ(t) = dt/(et+ 1) involving the Fermi function 1 /(et+ 1) (also used in solid state physics). Example 17. The series S=∞/summationdisplay k=1(−1)k−1 ke−1/k. One can show that the function fin question here is f(t) =J0(2√ t), with J0the Bessel function of order zero. Errors obtained by the thr ee quadrature methods are displayed in the table below, showing the clear s uperiority of Gauss–Fermi quadrature. 86 nGauss-Laguerre rational Gauss-Laguerre Gauss-Fermi 1 1.6961(–01) 1.0310(–01) 5.6994(–01) 4 4.4754(–03) 4.6605(–05) 9.6454(–07) 7 1.7468(–04) 1.8274(–09) 9.1529(–15) 10 3.7891(–06) 1.5729(–13) 2.8163(–16) 15 2.6569(–07) 1.5490(–15) 20 8.6155(–09) 40 1.8066(–13) timing: 12.7 timing: 19.5 timing: 4.95 20.3 Series generated by the derivative of a Laplace transform These are series (20.1) in which (20.6) ak=−d ds(Lf)(s)/vextendsingle/vextendsingle/vextendsingle/vextendsingle s=k, k= 1,2,3, . . . . In this case one ﬁnds (20.7) S=/integraldisplay∞ 0t 1−e−tf(t)e−tdt=/integraldisplay∞ 0f(t)t et−1dt, and Gauss–Laguerre, rational/polynomial Gauss–Laguerre , and Gauss–Einstein quadrature are again options as in §20.1. Example 18. The series S=∞/summationdisplay k=1(3 2+ 1)k−2(k+ 1)−3/2. The relevant function fis calculated to be f(t) =erf√ t√ t·t1/2, where erf is the error function erf x= (2/√π)/integraltextx 0e−t2dt. Numerical results analogous to those in the two previous tables are shown below . 87 nGauss-Laguerre rational Gauss-Laguerre Gauss-Einstein 1 4.0125(–03) 5.1071(–02) 8.1715(–02) 4 1.5108(–05) 4.5309(–08) 1.6872(–04) 7 4.6576(–08) 1.3226(–13) 3.1571(–07) 10 3.0433(–09) 1.2087(–15) 5.4661(–10) 15 4.3126(–11) 1.2605(–14) 20 7.6664(–14) 30 3.4533(–16) timing: 6.50 timing: 10.8 timing: 1.58 The run time is best for Gauss–Einstein quadrature and the er ror only slightly worse than for the closest competitor, rational Gauss–Lagu erre. 20.4 Slowly convergent series occurring in plate con- tact problems The series of interest here is (20.8) Rp(z) =∞/summationdisplay k=0z2k+1 (2k+ 1)p, z∈C,|z| ≤1, p= 2 or 3 . Rather than expressing the whole general term of the series a s a Laplace transform, we do this only for the coeﬃcient, (20.9)1 (k+1 2)p= (Lf)(k), f(t) =1 (p−1)!tp−1e−t/2. Then Rp(z) =z 2p∞/summationdisplay k=0z2k (k+1 2)p =z 2p∞/summationdisplay k=0z2k/integraldisplay∞ 0e−kt·tp−1e−t/2 (p−1)!dt =z 2p(p−1)!/integraldisplay∞ 0∞/summationdisplay k=0(z2e−t)k·tp−1e−t/2dt =z 2p(p−1)!/integraldisplay∞ 01 1−z2e−ttp−1e−t/2dt, that is, (20.10) Rp(z) =z 2p(p−1)!/integraldisplay∞ 0tp−1et/2 et−z2dt, z−2∈C\[0,1]. 88 The case z= 1 can be treated directly by using the connection with the ze ta function, Rp(1) = (1 −2−p)ζ(p). Assume therefore z/negationslash= 1. When |z|is close to 1, the integrand in (20.10) is rather ill-behaved near t= 0, exhibiting a steep boundary layer. We try to circumvent this by making the change of variables e−t/mapsto→tto obtain Rp(z) =1 2p(p−1)!z/integraldisplay1 0t−1/2[ln(1/t)]p−1 z−2−tdt. This expresses Rp(z) as a Cauchy integral of the measure dλ[p](t) =t−1/2[ln(1/t)]p−1dt. Since by assumption, z−2lies outside the interval [0 ,1], the integral can be evaluated by the continued fraction algorithm, once suﬃc iently many recurrence coeﬃcients for d λ[p]have been precomputed. For the latter, the modiﬁed Chebyshev algorithm is quite eﬀective. The ﬁrst 100 coeﬃcients are available for p= 2 and p= 3 in the OPQﬁlesabsqm1log1 andabsqm1log2 to 25 resp. 20 decimal digits. Example 19. Rp(x), p= 2 and 3 , x=.8, .9, .95, .99, .999 and 1 .0. Numerical results are shown in he table below. x R 2(x) R3(x) .8 0.87728809392147 0.82248858052014 .9 1.02593895111111 0.93414857586540 .95 1.11409957792905 0.99191543992243 .99 1.20207566477686 1.03957223187364 .999 1.22939819733 1.05056774973 1.000 1.233625 1.051795 Full acuracy cannot be achieved for x≥.999 using only 100 recurrence coeﬃcients of d λ[p]. Example 20. Rp(eiα), p= 2 and 3 , α=ωπ/2, ω=.2, .1, .05, .01, .001 and 0 .0. Numerical results are shown in the table below. 89 p ω Re(Rp(z)) Im( Rp(z)) 2 .2 0.98696044010894 0.44740227008596 3 0.96915102126252 0.34882061265337 2 0.1 1.11033049512255 0.27830297928558 3 1.02685555765937 0.18409976778928 2 0.05 1.17201552262936 0.16639152396897 3 1.04449441539672 0.09447224926029 2 0.01 1.22136354463481 0.04592009281744 3 1.05140829197388 0.01928202831056 2 0.001 1.232466849 0.006400460 3 1.051794454 0.001936923 2 0.000 1.2336 0.0000 3 1.0518 0.0000 Here, too, full accuracy is not attainable for ω≤0.001 with only 100 recur- rence coeﬃcients. 20.5 Series involving ratios of hyperbolic functions More of a challenge are series of the type (20.11) Tp(x;b) =∞/summationdisplay k=01 (2k+ 1)pcosh(2 k+ 1)k cosh(2 k+ 1)b,0≤x≤b, b > 0, p= 2,3, which also occur in plate contact problems. Here, we ﬁrst exp and the ratio of hyperbolic cosines into an inﬁnite series, (20.12)cosh(2 k+ 1)x cosh(2 k+ 1)b =∞/summationdisplay n=0(−1)n/braceleftbig e−(2k+1)[(2 n+1)b−x]+ e−(2k+1)[(2 n+1)b+x]/bracerightbig , insert this in (20.11) and apply the Laplace transform techn ique of §20.4. This yields, after an elementary computation (using an inte rchange of the summations over kandn), (20.13) Tp(x, b) =1 2p(p−1)!∞/summationdisplay n=0(−1)ne(2n+1)b[ϕn(−x) +ϕn(x)], 90 where (20.14) ϕn(s) = es/integraldisplay1 0dλ[p](t) e2[(2n+1)b+s]−t,−b≤s≤b. The integral on the right is again amenable to the continued f raction algo- rithm for d λ[p]. 91 Exercises to Part III (Stars indicate more advanced exercises.) 1. With π0, π1, . . ., π N−1denoting the discrete orthogonal polynomials rel- ative to the measure d λN, and ˆ ci(f) the Fourier coeﬃcients of fwith respect to these orthogonal polynomials, show that n/summationdisplay i=0|ˆci(f)|2/bardblπi/bardbl2≤ /bardblf/bardbl2, n < N, with equality holding for n=N−1. 2. Prove the following alternative form for the Fourier coeﬃ cients, ˆci(f) =1 /bardblπi/bardbl2/parenleftBigg f−i−1/summationdisplay j=0ˆcj(f)πj, πi/parenrightBigg , i= 0,1, . . ., n, and discuss its possible advantages over the original form. 3. Discuss the modiﬁcations required in the constrained lea st squares ap- proximation when ν(0≤ν≤m) of the points sjare equal to one of the support points tk. 4. What are pm(f;·),f∗, and σmin Example 11? 5. Calculate the ﬁrst and second derivative of the complemen tary error function of Example 12. 6∗. Prove the unique solvability of the problem (19.8) under th e conditions stated in (19.9)–(19.10), and, in the aﬃrmative case, deriv e (19.11). 7. Derive the measure d λ[m]for the Maxwell distribution of Example 14. 8. Derive the formula for fin Example 15. 9. Derive the formula for fin Example 16. 10. Derive (20.5). 11. Derive the formula for fin Example 17. 12. Derive (20.7). 13. Derive the formula for fin Example 18. 14. Supply the details for deriving (20.13). 92