/* Some Maxima functions for playing with representation of functions of a single "real" variable (double-precision floats, actually) by their Chebyshev series coefficients. This particular version is set up for generation and manipulation of functions defined on [-1,1]. We refer elsewhere (insert good reference) for a discussion of the nature of functions that can be appropriately modeled by a Chebyshev expansion. The functions could be translated or scaled, or pasted together piecewise. as done by the later work of the Chebfun group (see their web page). It has been said that "the Chebyshev polynomials are everywhere dense in numerical analysis." |# To Do: 1. I haven't debugged how to use FFT for this, so indeed this is probably much slower than it could be. The discrete cosine transform of size n is equivalent to an FFT of size 4 n where half the entries are zero and half of the remaining entries are the negatives of the original. There are so many implementations of the FFT that I could not be sure of using the best one, but if I pick one, (including my own), for sure that would not be best. The Chebfun people presumably finess this by using whatever Matlab uses, which is probably some version from FFTW. (which includes versions of the DCT as well). And then I would feel some obligation to deal with the "give special consideration to powers of two". Possibilities for fastest include DJ Bernstein's code, djbfft, FFTPACK, FFTW. And I don't like to irritate friends here. So all-in-all, I don't feel so bad about using one or two line programs, even if they are n^2 rather than nlogn. But if I find one "fast" implementation that works without effort and can be called from Lisp, and it falls into my lap, I might use it. 2. I am probably naive in the calculation heuristics for how many terms are needed to represent a function. 3. I have not dealt at all with NaN or Infinities. 4. There are many extensions (cf chebfun group software in Matlab) that are possible. For example, using splitting of functions to handle discontinuities. Using arbitrary or even infinite domains. 5. Many functions that are supplied could be done in a cleverer fashion. For example, conversion of rational functions can be done by partial fraction expansion, rootfinding, and many special cases. We can multiply directly in chebyshev series, and (in the works) do function composition. 6. Many functions that are described in chebfun literature are not implemented at all. Rootfinding, quadrature for a start. 7. These functions are not using the full generic arithmetic facilities. Some assume (or insist on) double-float coefficients. The package could be make more generic to allow quad-double, or arbitrary precision (software) floats. There is less argument for using (say) rationals or intervals, and it is probably worthwhile to be fast. On the other hand, some functions are written to produce exact chebyshev coefficients when given polynomial inputs. Even polynomials with symbolic entries. */ /*going back to Einwohner (& Fateman) paper in 1989. Take an expression and convert to Cheby form via Taylor Series. */ kill(tc,FromCheby, ToCheby, ToCheby2)$ tc[0](x):=1$ tc[1](x):=rat(x)$ tc[n](x):=2*x*tc[n-1](x)-tc[n-2](x)$ /*here's another way, if you are in a hurry */ tc2[0](x):=1$ tc2[1](x):=x$ tc2[n](x):= if evenp(n) then 2*tc2[n/2](x)^2-1 else 2*tc2[(n+1)/2](x)*tc2[(n-1)/2](x)-x$ ToCheby(e,x,n):= /* convert analytic expression, depending on x, to Chebyshev series. The expression e may be a polynomial or function in other variables, too. The result will be a list of the Chebyshev coefficients. Here we try to make them all floats. */ block([ tay: taylor(e,x,0,n),q,r,ans:[], keepfloat:true], tay:ratdisrep(tay), for i: n step -1 thru 0 do ([q,r]:divide(tay,tc[i](x),x), ans:cons(float(q),ans),tay:r), return(apply('chebseries,ans)))$ ToChebyt(e,x,n):= /* convert analytic expression, depending on x, to Chebyshev series, via taylor series */ block([ tay: taylor(e,x,0,n),q,r,ans:[]], power2cheby(makelist(ratcoef(tay,x,i),i,0,n))) /* another version, which has been written in lisp for maxima */ /* convert analytic function lambda([x]...expression in x) to Chebyshev series The expression f may be a polynomial or function in other variables, too. The result will be a list of the Chebyshev coefficients. Here we try to make them all floats. */ /*BROKEN*/ ToCheby2(%%f,n):= /* this one takes a function, e,g, lambda([x],..) */ apply ('chebseries, DCT(makelist(apply(%%f,[(cos(?pi*(k+0.5d0)/n))]),k,0,n-1)))$ /* gauss-lobatto points*/ glp(n):= makelist(cos(?pi*(i+0.5d0)/n),i,0,n-1)$ glpf(f,n):= makelist(apply(f, [cos(?pi*(i+0.5d0)/n)]),i,0,n-1)$ FromCheby(a,x):= /* If x is a symbol, this will convert Chebyshev series to an expression in the power basis, i.e. polynomial in x. If x is a number, this will evaluate the series at that point and provide a number. (if a is just numbers). This algorithm is based on recurrence in Clenshaw alg. */ block([N:length(a)-1,keepfloat:true, bj2:0,bj1:0,bj0:0, z:reverse(a)], for j: 1 thru N do (bj0:2*x*bj1-bj2+first(z), bj2:bj1, bj1:bj0,z:rest(z)), x*bj1+first(z)/2 -bj2)$ /* a test */ FromCheby(ToCheby(a*x^5+b*x^2+45,x,5),rat(x)); /* Given a Chebyshev series, return another one in which we may have removed terms that don't matter much; that is, trailing coeffs were small relative to large coeffs. */ TrimCS(a,relerr):= block([max: apply(max,args(a)), min: apply(min,args(a)), r: reverse(a), bigmaginv], bigmaginv: 1/max(abs(max),abs(min)), while (abs(first(r))*bigmaginv < relerr) do r:rest(r), return(reverse(r)))$ small_to_zero(expr,eps):= /* Replaces all numbers smaller in absolute value than EPS by zero . */ fullmap(lambda([x], if numberp(x) and abs(x)