Algorithms for Polynomial and Rational Approximation

141 Pages·2014·2.57 MB·English

Algorithms for Polynomial and
Rational Approximation
(cid:0)
Ricardo Pacho´n
Exeter College
University of Oxford
A thesis submitted for the degree of
Doctor of Philosophy
Trinity Term 2010
a mi hermana
Acknowledgements
I am deeply grateful to Nick Trefethen, an inspiring supervisor and a mentor
since I arrived to Oxford. It has been a privilege to be his student and I am
indebted to him for his advice and encouragement through this project.
I had the pleasure of working very closely with postdoctoral researchers Rodrigo
Platte, Pedro Gonnet and Joris van Deun. My sincere thanks go to the three of
themforsharingtheirknowledgeandcontributingtomyownresearch. Together
withA´sgeirBirkisson,TobyDriscoll,NickHaleandMarkRichardson,theywere
mycolleaguesintheChebfunTeam,ofwhichIamproudtohavebeenamember.
My gratefulness goes to all of them.
I offer my thanks and enduring gratitude to everyone in the Numerical Analysis
GroupatOxford. InparticularIwouldliketoacknowledgeNickGould,Raphael
Hauser, Endre Su¨li and Andy Wathen, and my friends and colleagues Siobhan
Burke, Leonardo Figueroa, Jari Fowkes, Bernhard Langwallner, Tyrone Rees,
Thomas Schmelzer and Martin Stoll. A special mention goes to Lotti Ekert, for
her invaluable help on a day-to-day basis.
My gratitude goes to Jean Paul Berrut and Ian Sobey, the examiners of this
thesis, who provided me with thorough comments and suggestions to improve
theoriginalmanuscript. Needlesstosay, anyerrorsoroversightsareminealone.
IwasoneofthefoundersoftheUniversityofOxfordSIAMStudentChapterand
a member of the committee for two years. This was a tremendously enriching
experience and I would like to acknowledge David Hewett, David Knezevic,
Sarah McBurnie, Mason Porter and Thomas Woolley.
My sincere gratitude to all my friends at Oxford, many at Exeter College, for
making the last four years a unique and vibrant experience.
I would like to acknowlege the funding by the Engineering and Physical Sci-
ences Research Council, Oxford University Computing Laboratory, Universidad
Cat´olica de Colombia, and Exeter College.
Finally, all my love goes to my father, my mother and my sister, who have been
supporting me at every step since I left Colombia six years ago.
Abstract
Robustalgorithmsfortheapproximationoffunctionsarestudiedanddeveloped
in this thesis. Novel results and algorithms on piecewise polynomial interpola-
tion, rational interpolation and best polynomial and rational approximations
are presented.
Algorithms for the extension of Chebfun, a software system for the numerical
computationwithfunctions,aredescribed. Thesealgorithmsallowtheconstruc-
tion and manipulation of piecewise smooth functions numerically with machine
precision. Breakpoints delimiting subintervals are introduced explicitly, implic-
itly or automatically, the latter method combining recursive subdivision and
edge detection techniques.
For interpolation by rational functions with free poles, a novel method is pre-
sented. When the interpolation nodes are roots of unity or Chebyshev points
the algorithm is particularly simple and relies on discrete Fourier transform ma-
trices, which results in a fast implementation using the Fast Fourier Transform.
The method is generalised for arbitrary grids, which requires the construction
of polynomials orthogonal on the set of interpolation nodes. The new algo-
rithm has connections with other methods, particularly the work of Jacobi and
Kronecker, Berrut and Mittelmann, and E˜gecio˜glu and Ko¸c.
Computedrationalinterpolantsarecomparedwiththebehaviourexpectedfrom
the theory of convergence of these approximants, and the difficulties due to
truncated arithmetic are explained. The appearance of common factors in the
numerator and denominator due to finite precision arithmetic is characterised
by the behaviour of the singular values of the linear system associated with the
rational interpolation problem.
Finally, new Remez algorithms for the computation of best polynomial and
rational approximations are presented. These algorithms rely on interpolation,
for the computation of trial functions, and on Chebfun, for the location of trial
references. For polynomials, the algorithm is particularly robust and efficient,
andwereportexperimentswithdegreesinthethousands. Forrationalfunctions,
we clarify the numerical issues that affect its application.
Contents
1 Introduction 1
2 Approximation by Polynomial Interpolants 8
2.1 Computation of polynomial interpolants . . . . . . . . . . . . . . . . . . . . 10
2.1.1 Monomials and interpolation in the unit circle . . . . . . . . . . . . 11
2.1.2 Chebyshev polynomials and interpolation in a real interval . . . . . 13
2.1.3 Lagrange interpolation and the barycentric formula . . . . . . . . . . 17
2.1.4 Newton interpolation. . . . . . . . . . . . . . . . . . . . . . . . . . . 20
2.2 Convergence and divergence of polynomial interpolants . . . . . . . . . . . . 20
2.2.1 Illusory interpolants . . . . . . . . . . . . . . . . . . . . . . . . . . . 21
2.2.2 Convergence for analytic functions . . . . . . . . . . . . . . . . . . . 23
2.2.3 Convergence for piecewise analytic functions . . . . . . . . . . . . . 28
2.2.4 The Gibbs phenomenon . . . . . . . . . . . . . . . . . . . . . . . . . 30
2.2.5 Number of points per wavelength . . . . . . . . . . . . . . . . . . . . 31
2.3 Lebesgue constants and truncations of series . . . . . . . . . . . . . . . . . . 32
2.4 Continuation from polynomial interpolants . . . . . . . . . . . . . . . . . . 34
2.5 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35
3 Chebfuns and Piecewise-Smooth Chebfuns 37
3.1 Smooth chebfuns . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40
3.2 Piecewise smooth chebfuns . . . . . . . . . . . . . . . . . . . . . . . . . . . 41
3.3 Automatic subdivision . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43
3.4 Other commands . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48
3.5 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48
i
4 A New Method for Rational Interpolation 50
4.1 Rational interpolation in roots of unity . . . . . . . . . . . . . . . . . . . . . 52
4.1.1 Derivation of the method . . . . . . . . . . . . . . . . . . . . . . . . 52
4.1.2 Implementation with FFTs and barycentric formula . . . . . . . . . 54
4.2 Rational interpolation in other grids . . . . . . . . . . . . . . . . . . . . . . 56
4.2.1 Chebyshev points . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57
4.2.2 Equidistant points . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59
4.2.3 Arbitrary Grids . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61
4.3 History and Related Work . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61
4.4 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 64
5 Approximation by Rational Functions 66
5.1 Convergence analysis of rational interpolation . . . . . . . . . . . . . . . . . 68
5.2 Practical considerations in approximation by rational interpolants . . . . . 77
5.3 Pad´e and Chebyshev–Pad´e approximations . . . . . . . . . . . . . . . . . . 83
5.4 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 89
6 Algorithms for Best Polynomial and Rational Approximation 92
6.1 Overview of the Remez algorithm . . . . . . . . . . . . . . . . . . . . . . . . 95
6.2 Interpolation and best approximation: From a reference to a trial function . 98
6.3 Location of extrema: From a trial function to a reference . . . . . . . . . . 102
6.4 Implementation for best polynomial approximation . . . . . . . . . . . . . . 103
6.5 Implementation for best rational approximation . . . . . . . . . . . . . . . . 113
6.6 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 118
7 Conclusion 120
Bibliography 123
ii
List of Figures
1.1 Best polynomial approximants of degrees 5, 7, 9 and 11 to x .. . . . . . . . 2
| |
1.2 Error of the best polynomial approximant of degree 100 to exp( x ). . . . . 3
| |
1.3 Structure of the thesis. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4
2.1 Contour plots of T (z) 1/N = constant, with N = 5,9,17 and as N . . 15
N
| | → ∞
2.2 Conditioning of Vandermonde and Vandermonde-type matrices. . . . . . . . 16
2.3 Partial products of the polynomial barycentric weights. . . . . . . . . . . . 19
2.4 Errors of polynomial approximants to the function (1+25x2)−1. . . . . . . 23
2.5 Contour plots of (cid:96)(z) 1/N = constant, with N = 4,8,16 and as N . . . 28
| | → ∞
2.6 Hermite integral contours in the proof of Theorem 2.8 for convergence of
polynomial interpolants to piecewise analytic functions. . . . . . . . . . . . 30
3.1 Plot of a chebfun length 265 . . . . . . . . . . . . . . . . . . . . . . . . . . 39
3.2 Chebyshev coefficients of intermediate interpolants when computing a chebfun. 41
3.3 Example of a piecewise smooth chebfun. . . . . . . . . . . . . . . . . . . . . 43
3.4 Algorithm for the computation of piecewise smooth chebfuns using recursive
subdivision. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45
3.5 Algorithm for the computation of piecewise smooth chebfuns using recursive
subdivision and edge-detection. . . . . . . . . . . . . . . . . . . . . . . . . . 46
3.6 Algorithm for the detection of edges in Chebfun (proposed by Platte). . . . 47
4.1 Matlab code for rational interpolation in roots of unity. . . . . . . . . . . . 55
4.2 Matlab code for rational interpolation in Chebyshev points (1st kind). . . . 58
4.3 Matlab code for rational interpolation in Chebyshev points (2nd kind). . . . 60
5.1 Number of publications in MathSciNet on rational approximation 1970–2010. 68
5.2 Convergence of rational interpolants in Chebyshev points with fixed degree
of the denominator to a meromorphic function. . . . . . . . . . . . . . . . . 70
iii
5.3 Convergence of rational interpolants in roots of unity with fixed degree of
the denominator to a meromorphic function. . . . . . . . . . . . . . . . . . 72
5.4 Analytic continuation through rational interpolation in Chebyshev points. . 73
5.5 Analytic continuation through rational interpolation in roots of unity. . . . 74
5.6 Rational interpolation error to x on grids of Chebyshev, Newman and ad-
| |
justed Chebyshev points.. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 75
5.7 Gibbs phenomenon for rational interpolation in Chebyshev points. . . . . . 76
5.8 Example of a rational interpolant with poles in the domain of definition. . . 78
5.9 Comparison between computed and theoretical rational interpolants. . . . . 80
5.10 Error of diagonal rational interpolants in Chebyshev points (1st and 2nd kind). 81
5.11 Singular values of Z with monomials and Chebyshev polynomials. . . . . . 82
5.12 Error of rational interpolant with reduced type. . . . . . . . . . . . . . . . . 83
5.13 Error of Chebyshev–Pad´e approximation (Maehly type). . . . . . . . . . . . 86
5.14 Error of Chebyshev–Pad´e approximation (Clenshaw–Lord type). . . . . . . 88
5.15 Matlab code for Chebyshev–Pad´e approximation. . . . . . . . . . . . . . . . 89
6.1 Best polynomial approximants to sin(3πx)exp(x). . . . . . . . . . . . . . . 96
6.2 Best rational approximants to tanh(50x) (diagonal). . . . . . . . . . . . . . 97
6.3 Best rational approximants of types to tanh(50x) (anti-diagonal). . . . . . . 97
6.4 Four steps of the Remez algorithm for the best approximation of type [5/5]
to the gamma function. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 99
6.5 Matlab code of the exchange algorithm in the command remez in Chebfun. 104
6.6 Matlab code of the command remez in Chebfun. . . . . . . . . . . . . . . . 105
6.7 Best polynomial approximations of degree 10 to the functions in Table 6.1. 108
6.8 Errors of best polynomial approximants to f in Table 6.1.. . . . . . . . . . 109
6
6.9 Error of best polynomial approximant of degree 1000 to f in Table 6.1. . . 110
6
6.10 Computed errors of best polynomial approximation for Lipschitz continuous
functions. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 110
6.11 Computed errors of best polynomial approximation for continuous functions. 111
6.12 Computed values for the disproof of Bernstein’s conjecture. . . . . . . . . . 112
6.13 Matlab code of the command rationalremez.. . . . . . . . . . . . . . . . . 114
6.14 Chebyshev–Pad´e approximants to the gamma function. . . . . . . . . . . . 116
6.15 Best rational approximations to the gamma function. . . . . . . . . . . . . . 117
6.16 Rows and columns of the Walsh table for the gamma function. . . . . . . . 118
6.17 Walsh table for the gamma function. . . . . . . . . . . . . . . . . . . . . . . 119
iv
1
CHAPTER
Introduction
Although this may seem a paradox, all exact science is dominated by the idea
of approximation. When a man tells you that he knows the exact truth about
anything, you are safe in inferring that he is an inexact man.
Bertrand Russell
The Scientific Outlook (1931)
In 1934, Evgeny Yakovlevich Remez published in a series of three papers the algorithm
that now bears his name for the computation of the best polynomial approximation, that
is, the polynomial among those of a fixed degree that minimises the deviation in the supre-
mum norm from a given continuous function on a given interval. We learn from Remez
himself, in a footnote in one of these publications, that the first users of his algorithm
were three female students at the University of Kiev—Mesdemoiselles Linnik, Kouniavska,
and Masanovska. This revelation is not surprising, coming from a time when the role of
women in mathematics was largely confined to that of a human computer, workers with the
heavy labour of performing long and tedious numerical calculations for scientific purposes.
The three women obtained the coefficients and errors of the best approximations to x
| |
by polynomials of degrees 5, 7, 9, and 11, and though we do not know how tiresome this
computation was, we know it was carried out meticulously, yielding all the figures accurate
to about 4 decimal places.
Seventy-five years after these numerical experiments were performed, how large can
we make the degree when computing best polynomial approximations? The advent of
computers in the mid-20th century ignited interest in the Remez algorithm and we find
numerous publications about it in the late 1950s and early 1960s, reporting examples with
degrees of a couple of dozens. However, after this burst of investigation, research on the
computation of best approximations diminished, perhaps in part due to the paucity of
applications aside from digital filter design, a subject where it found its niche. The Remez
1
1 1
N = 5 N = 7
0.8 0.8
0.6 0.6
0.4 0.4
0.2 0.2
0 0
−1 −0.5 0 0.5 1 −1 −0.5 0 0.5 1
1 1
N = 9 N = 11
0.8 0.8
0.6 0.6
0.4 0.4
0.2 0.2
0 0
−1 −0.5 0 0.5 1 −1 −0.5 0 0.5 1
Figure 1.1: Best polynomial approximation (thin lines) of degrees 5, 7, 9, and 11 to x on [ 1,1].
| | −
Thecoefficientsanderrorsofthesepolynomials,accurateto4decimalplaces,werecomputedforthe
first time by Linnik, Kouniavska, and Masanovska, and published by Remez in 1934. The methods
presentedinthisthesisallowonetocomputebestapproximantstonumerousfunctionswithdegrees
in the thousands, typically accurate to 12–13 decimal places, in a matter of seconds.
algorithm became a mandatory topic in books of approximation theory but distant from
practical usage, and the exploration of best polynomials with degrees in the hundreds or
thousands never became a priority.
At the end of this thesis, in Chapter 6, we will present a new version of the Remez
algorithm that makes it possible to compute best approximations with a single Chebfun
command. For example, here we find the degree 100 best approximation to exp( x ) in less
| |
than a second:
>> f = chebfun(@(x) exp(abs(x)),’splitting’,’on’);
>> tic, p = remez(f,100); toc
Elapsed time is 0.307301 seconds.
>> norm(f-p,inf)
ans = 0.0028014408940777
Withonefurthercommand,plot(f-p),weobtainacurvewith103pointsofequioscillation.
The simplicity and speed of these computations changes the nature of the exploration of
the Remez algorithm, making it an easy matter to quickly compare best approximants of
many degrees, or for various different functions. See Figures 1.1 and 1.2.
This new way of computing best polynomial approximations is a good example of the
main theme in this thesis: The study and development of robust algorithms for the ap-
proximation of functions. In the rest of this chapter we give a panoramic view of the main
2