;;; -*- Mode: LISP; Package: Maxima; Syntax: Common-Lisp; ;;; a package for chebyshev approximations. Comments at end (in-package :maxima) ;; tocheb computes all the a_j, j=0 to n-1 for the Maxima function f. (defun $tocheb (f n) (cons '($chebseries) (ajs f n))) ;; compute a[j]s for Maxima this way (defun ajs(f n)(loop for j from 0 to (1- n) collect (aj f j n))) ;; it turns out we will be re-using certain key values repeatedly, ;; so we save them in a hash table. (defvar *cosc* (make-hash-table)) ;; store useful values of cosine here. (defun cachecos(r)(or (gethash r *cosc*)(setf (gethash r *cosc*) (cos r)))) (defun aj (f j n) ;; compute just one coefficient, a_j out of n for function f (let ((in (/ 1.0d0 n))) ; inverse of n (* in 2 (- 2 (if (= j n) 0 1)) (loop for k from 0 to (1- n) sum (let ((h(* pi (+ k 0.5d0) in))) (* (cachecos (* h j)) (mapply1 f (list (cachecos h)) f f))))))) (defun $cheb_call(a x) ;; clenshaw algorithm ;; evaluate p=sum'(a_nT_n(x), i=0 to N) a is chebseries ;; sum' multiplies a_0 by 1/2 ;; works. (declare (optimize (speed 3)(safety 0)) (double-float x)) (cond ((eq (caar a) '$chebseries) (let* ((bj2 0.0d0) (bj1 0.0d0) (bj0 0.0d0)) (declare(double-float bj2 bj1 bj0)) (setf a (reverse(cdr a))) ;; could hack this away. maybe later (loop while (cdr a) do (setf bj0 (+ (* 2.0d0 x bj1) (- bj2) (the double-float (pop a)))) (setf bj2 bj1 bj1 bj0)) (+ (* x bj1) (* 0.5d0(pop a))(- bj2)))) (t (error "expecting chebseries,not ~s" a)))) (defun $chebadd(a b) ;; a slower version, chebplus, is below (let* ((ans nil)) (cond ((> (length a)(length b)) (setf a (append (cdr a) nil)) ;copy (setf ans a) (setf b (cdr b)) (while b (incf (car a)(pop b)) (pop a))) (t (setf b (append (cdr b) nil)) ;copy (setf ans b) (setf a (cdr a)) (while a (incf (car b)(pop a)) (pop b)))) (cons '($chebseries) ans))) (defun $chebmul(a b) ;; assumes each is a chebseries. no error checking here. ;; a slower version is below. (let* ((aa (coerce (rest a) 'simple-array)) (bb (coerce (rest b) 'simple-array)) (la (length aa)) (lb (length bb)) (term 0.0d0) (ans (make-array (+ la lb -1) :initial-element 0.0d0))) (setf (aref aa 0)(* 0.5d0 (aref aa 0))) (setf (aref bb 0)(* 0.5d0 (aref bb 0))) (loop for i from 0 to (1- la) do (loop for j from 0 to (1- lb) do (setf term (* (aref aa i)(aref bb j))) (incf (aref ans (+ i j)) term) (incf (aref ans (abs(- i j))) term) )) (loop for i from 1 to (+ la lb -2) do;; not ans[0] (setf (aref ans i) (* (aref ans i)0.5d0)) ) (cons '($chebseries)(coerce ans 'list)))) ;; how about composition? (defun $chebcomp(f g n) ;; f, g are chebseries ;; simple hack is to evaluate f(g(x)) as needed, e.g. ($tocheb #'(lambda(r)($cheb_call f ($cheb_call g r))) ;; how many points are justified though?? n)) (defun $chebplus(f g) ;; f, g are chebseries ;; simple hack is to evaluate f and g as needed, e.g. ($tocheb #'(lambda(r)(+($cheb_call f r) ($cheb_call g r))) ;; how many points are justified though?? (max (length f)(length g)))) ;; this works (defun $chebtimes(f g) ;; f, g are chebseries ;; simple hack is to evaluate f and g as needed, e.g. ($tocheb #'(lambda(r)(*($cheb_call f r) ($cheb_call g r))) ;; how many points are justified though?? (+ (length f)(length g) -1))) (defun $chebinv(f n) ;; f is a chebseries. Compute 1/f ;; simple hack is to evaluate f and g as needed, e.g. ($tocheb #'(lambda(r)(/ 1.0d0 ($cheb_call f r))) ;; how many points are justified though?? n)) (defun $chebapply(fun cs n) ;; apply any maxima function to a chebseries ;; e.g. chebapply(lambda([r](exp(r^2+1)),chebs); ;; constraint: fun must be a function of one float argument that ;; returns a double-float. ($tocheb #'(lambda(r) (mapply1 fun (list r) fun fun)) cs ;; how many points are justified though?? n)) (defun $chebint(a) ;; indefinite integration of the function given by a. ;;integrate(t[n](x),x) = ;; if n=0 then t[1]. ;; if n=1 then 1/2*t[0] +1/4*t[2] else ;; 1/2*(t[n+1]/(n+1) - t[n-1]/(n-1)) ;; ;; (let* ((aa (coerce (rest a) 'simple-array)) (la (length aa)) (term 0.0d0) (ans (make-array (max 2(1+ la)) :initial-element 0.0d0))) (setf (aref aa 0) (* 0.5d0 (aref aa 0))) (setf (aref ans 1) (aref aa 0)) ;int(const)= const*x (setf (aref ans 0) (* 0.25d0 (aref aa 1))) ;int(x) = x^2/2 =1/4*T[0]+ (setf (aref ans 2) (* 0.25d0 (aref aa 1))) ;...+ 1/4*T[2] (loop for i from 2 to (1- la) do (setf term (* 0.5 (aref aa i))) (incf (aref ans (1+ i)) (/ term (1+ i))) (decf (aref ans (1- i)) (/ term (1- i)))) (setf (aref ans 0)(* 0.5d0 (aref ans 0))) ;so sum' works (cons '($chebseries)(coerce ans 'list)))) ;; elaboration: fix up the above programs to combine non-chebseries with ;; chebseries. Compute the appropriate number of terms. #| This file defines 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 (http://www2.maths.ox.ac.uk/chebfun) for a discussion of the kinds of computations that can be considered, and the nature of functions that can be appropriately modeled by a Chebyshev expansion. The chebfun web site discusses a Matlab project which includes facilities for functions being translated or scaled, or pasted together piecewise. A key piece (not included here), is a method for truncation of these chebyshev series, based on the relative size of the trailing coefficients. Drop em if them are too small to matter. Also, these series can be directly multiplied, added, composed. There are neat ways of computing them for polynomials and rational functions, though not as briefly as here. And there is a faster way of computing the coefficients from the evaluation points using an FFT (or a Discrete Cosine Transform), which we do not include because it would require other programs, and more debugging. These lisp programs should be compiled for speed. It has been said that "the Chebyshev polynomials are everywhere dense in numerical analysis." examples load("thisfile") s: tocheb(lambda([x],sin(x)),20) ; cheb_call(s,0.1234) ; sin(0.1234) ; gg(x):=sin(20*x^2+exp(x^5+sqrt(x^4+1))) ; g:tocheb(gg,60); gg(0.4) ; cheb_call(g,0.4) ; plot2d('[gg(x),cheb_call(g,x)],[x,0,1.025]) ; /* note they are the same until x>1 */ oh if you don't like s: tocheb(lambda([x],sin(x)),20) ; then we can do this: afun(ex,x):=buildq([expr:ex,var:x],lambda([var],expr)) $ cheb3(expr,var,count):= tocheb(afun(expr,var),count)$ cheb3(sin(z),z,20) then is the same as tocheb(lambda([z],sin(z)),20) which happens also to be the same as tocheb(sin,20). Here are a few additional ideas. Elaborate on chebseries so that they are callable. e.g. (($chebseries) a b c) becomes #'(lambda(r)(cheb_call '((chebseries) a b c) r)) in effect. applying a chebseries to an argument becomes simpler. include various limits (now default -1,1) package up a piecewise function as the chebfun project does. extract a coefficient from the series so that all coefficients are between -1 and 1. e.g. [1,2,3] becomes 3* [1/3,2/3,1]. The max coefficient would then always be 1. In fact, the coefficients could all be fixed point numbers between -.111111... and +.11111. Low precision would look like .0000001.. . Given that float arithmetic is ubiquitous and fast, this may not be a winner. other advantage: a method for truncating the product of chebyshev series. a= A* csa, where csa has max coefficient of +-1. [apparently not new,] b= B* csb, similarly. the product is then A*B* (csa*csb). without loss of generality, consider only csa*csb. We agree that if we know that all coefficients in the answer beyond some N have magnitude less than epsilon (say, 1e-13) then we will drop all those terms beyond N. We will not necessarily know this to be the case when we generate terms ab initio from some expression we are evaluating at Gauss-Lobatto points, but if we assert the two cs to be multiplied are in fact correct in all given terms, and furthermore, all terms beyond those explicitly represented are exactly zero, then we are doing about as well as can be expected with the information at our disposal. Typically the input cs will have rapidly declining coefficients, especially just before they are cut off. For purposes of being explicit, let us say that the coefficients start with magnitude 1 and decrease by a factor of 10 for each coefficient. that would suggest coefficients of magnitude 1,0.1, 0.01, ... 1e-12, 1e-13, 0.... in each input, and csa and csb are both of length 13. We know that any terms in the answer will be of magnitude 1e-13 or greater, so there is not much point in accumulating answers that are 1e-7 times terms in the other series that are 1e-7 or less. At worst, there will be 13 terms all of the same sign added together in a 13x13 multiply.. Consider csa = csa_big + sqrt(eps/n)*csa_small =A1+sqrt(eps/n)A2 n is length of csa csb = csb_big + sqrt(eps/m)*csb_small =B1+sqrt(eps/m)B2 m is length of csb The product contains A1*B1+sqrt(eps/m)*A1*B2 +sqrt(eps/n)*A2*B1+ eps/sqrt(nm)*A2*B2. This last component promises to be mostly negligible. Why multiply these? So instead of computing 13^2 =169 terms to produce 26 coefficients, we compute about 3/4 of them, or about 127 terms. We will also produce about 13 coefficients. maybe this cuts off too much? what if we use eps^2? a term of order eps^2 would not affect even the low order bits of a term of order eps, which is the smallest retained term, anyway. how would terms of order eps^2 be generated -- hardly at all. ugh. more thoughts, if we are going to make this into a generic system for dealing with functions. integral (just done) derivative zeros (eh, maybe) positive (eh, heuristic) test for equality majorization? f>g, f