Nick Trefethen, September 2010
(Chebfun example approx/ChebCoeffs.m)
Every function defined on [-1,1], so long as it is at least Lipschitz continuous, has an absolutely and uniformly convergent Chebyshev series:
f(x) = a_0 + a_1 T_1(x) + a_2 T_2(x) + ....
The same holds on an interval [a,b] with appropriately scaled and shifted Chebyshev polynomials.
For many functions you can compute these coefficients with the command CHEBPOLY. For example, here we compute the Chebyshev coefficients of a cubic polynomial:
x = chebfun('x'); format long disp('Cheb coeffs of 99x^2 + x^3:') p = 99*x.^2 + x.^3; a = chebpoly(p)'
Cheb coeffs of 99x^2 + x^3: a = 0.250000000000000 49.500000000000000 0.750000000000000 49.500000000000000
Notice that following the usual Matlab convention, the coefficients appear in order from highest degree to lowest. Thus it is often more useful to write
disp('Cheb coeffs of 99x^2 + x^3:') a = chebpoly(p)'; a = a(end:-1:1)
Cheb coeffs of 99x^2 + x^3: a = 49.500000000000000 0.750000000000000 49.500000000000000 0.250000000000000
Similarly, here are the Chebyshev coefficients down to level 1e-15 of exp(x):
disp('Cheb coeffs of exp(x):') a = chebpoly(exp(x))'; a = a(end:-1:1)
Cheb coeffs of exp(x): a = 1.266065877752008 1.130318207984970 0.271495339534077 0.044336849848664 0.005474240442094 0.000542926311914 0.000044977322954 0.000003198436462 0.000000199212481 0.000000011036772 0.000000000550590 0.000000000024980 0.000000000001039 0.000000000000040 0.000000000000001
You can plot the absolute values of these numbers on a log scale with CHEBPOLYPLOT:
FS = 'fontsize'; MS = 'markersize'; LW = 'linewidth'; chebpolyplot(exp(x),'.-',LW,1,MS,16), grid on xlabel('degree n',FS,12) ylabel('|a_n|',FS,12) title('Chebyshev coefficients of exp(x)',FS,16)
Here's a similar plot for a function that needs thousands of terms to be represented to 15 digits. (Can you explain why it looks like a wide stripe?)
chebpolyplot(exp(x)./(1+10000*x.^2)), grid on xlabel('degree n',FS,12) ylabel('|a_n|',FS,12) title('Chebyshev coefficients of exp(x)/(1+10000x^2)',FS,16)
These methods will work for any function f that's represented by a global polynomial, i.e., a chebfun consisting of one fun. (Normally this means that the Chebyshev series needs fewer than 65536 terms, Chebfun's default value of MAXDEGREE. To change this default, type HELP CHEBFUNPREF.) What about Chebyshev coefficients for functions that are not smooth enough for such a representation? Here one can use the TRUNC option in the Chebfun constructor. For example, suppose we are interested in the function
f = sign(x); figure, plot(f,'k',LW,2), ylim([-1.5 1.5]) title('sign(x)',FS,16)
If we try to compute all the Chebyshev coefficients, we'll get an error. On the other hand we can compute the first ten of them like this:
p = chebfun(f,'trunc',10); a = chebpoly(p)'; a = a(end:-1:1)
a = 0.000000000000000 1.273239544735164 -0.000000000000000 -0.424413181578388 -0.000000000000000 0.254647908947030 0.000000000000000 -0.181891363533601 0.000000000000000 0.141471060526123
Here's the degree 9 polynomial obtained by adding up these first terms of the Chebyshev expansion:
hold on plot(p,'m',LW,2) title('sign(x) and truncated Chebyshev series',FS,16)
This is not the same as the degree 9 polynomial interpolant through 10 Chebyshev points:
pinterp = chebfun(f,10); plot(pinterp,'--','color',[0 .8 0],LW,2) title('Same, also with Chebyshev interpolant',FS,16)
 L. N. Trefethen, Approximation Theory and Approximation Practice, draft book available at http://www.maths.ox.ac.uk/chebfun/ATAP/.