;;; -*- Mode: LISP; Syntax: Common-Lisp; -*- ;;; Barycentric interpolation representation of real-valued functions. ;;; This is a fast approximate way of computing with a real-valued ;;; function F of a single real value defined on the real interval ;;; [-1,1] by representing F at interpolated points. This is related ;;; to Chebyshev approximation of F, but is only a small piece of the ;;; total suite of programs that might be built up on this notion. ;;; Various generalizations abound, as can be seen by ;;; Trefethen's Chebfun Project, whose programs and documentation ;;; are far more elaborate, and whose papers explain how such notions can ;;; be used in many circumstances. Chebfun uses Matlab, rather than Lisp. ;;; The main function of interest is bcif, which produces an ;;; approximation of a function, f defined in Lisp R->R, of a ;;; single variable. The approximation is really only set up to be ;;; good in the interval [-1,1]. That is, f: [-1,1]->R ;;; Two input items are specified in bcif: the callable function f, ;;; and an integer k, probably less than 1024. That is ;;; (bcif f n x) computes the value of f at location x by ;;; interpolation at n points. If you have a rapidly varying function, ;;; or one that has "edges", you should use more points. ;;; Example.. (setf s (bcif(#'sin 8)), followed by (funcall s 0.5) should give ;;; approx value for sin(0.5). ;;; Here we set up a framework of the interpolation points, ;;; which we store in a global hash table *glpa*, which is filled up as needed. ;;; These key interpolation data points are known as gauss-lobatto points. ;;; Weights for interpolation are also stored. ;;; Oh, Barycentric interpolation is equivalent to Lagrange interpolation, but ;;; the computation is nicer for reasons you can read about elsewhere. (defvar *glpa* (make-hash-table)) (defvar *wta* (make-hash-table)) (defun glp (n) ;; return a vector of Gauss Lobatto points of length n. (or (gethash n *glpa*) ;; if we have already computed them, just return that entry. (setf (gethash n *glpa*) ;; otherwise compute GL points, store them, return them. (let((h (/ pi n))) (declare (double-float h)) (make-array n :element-type 'double-float :initial-contents (loop for i from 0 to (1- n) collect (cos(* h (+ i 0.5d0))))))))) ;; Note that glp[0]=-glp[n]. ;; Indeed, if n is even, glp[k]=-glp[n-k] for k=0..n/2-1 ;; If n is odd, glp[(n+1)/2]=0, and glp[k]=-glp[n-k] for k=0..(n+1)/2-1 (defun wta (n);; weight array (or (gethash n *wta*);; if we have already computed them, just return that entry. (setf (gethash n *wta*);; otherwise compute weights (let((h (glp n))) (make-array n :element-type 'double-float :initial-contents (loop for j fixnum from 0 to (1- n) collect (loop for i fixnum from 0 to (1- n) ;; very un-lisp-like code here... with wj double-float = (aref h j) with res double-float = 1.0d0 finally (return (/ 1.0d0 res)) unless (= i j) do (setf res (* res (- wj (aref h i))))))))))) ;; Note that wta[0]=-wta[n]. ;; Indeed, if n is even, wta[k]=-wta[n-k] for k=0..n/2-1 ;; if n is odd, wta[k]=-wta[n-k] for k=0..(n+1)/2-1 ;; (The midpoint weight is NOT zero.) ;; Make an interpolating function out of f: (defun bcif(f n) (let ((fa (make-array n :element-type 'double-float :initial-contents (map 'simple-array f (glp n)))) (g (glp n)) (w (wta n))) #'(lambda(x)(bcirun fa x g w)))) (defun bcirun(fa xin xa w) ;; bcirun uses an array of function-values of f, called fa. ;; and uses an array of weights. most likely already pre-computed, wta. ;; it evaluates f(xin) using barycentric interpolation from values in fa. ;; as a special trick, when xin=nil, bcirun returns the array fa, rather than a value. ;; this trick allows us to do arithmetic on interpolation arrays. See bciplus. (if (null xin) fa;; returns fa. Useful info from this function when x=nil (let*((sn 0d0)(sd 0d0)(term 0d0)(n (length fa)) (x (float xin 1.0d0)) ; make sure it is a double float (nm1 (1- n))) (declare (double-float sd sn term x)(fixnum n nm1) (type (simple-array double-float (*)) fa xa w)) (loop for i fixnum from 0 to nm1 do (if (= x (aref xa i)) (return-from bcirun (aref fa i))) (setf term (/ (the double-float (aref w i)) (- x (the double-float (aref xa i))))) (incf sn (* term (the double-float (aref fa i)))) (incf sd term)) (if (= sd 0.0d0) #+allegro #.excl::*nan-double* #-allegro (error "Division by zero in interpolation") (/ sn sd))))) ;; that's all that is needed. ;; a totally optional example function that shows we can do functional operations ;; on bci functions. Like adding them. We use a nil argument to extract the data needed ;; from inside the function. We redo this function in a more general setting, as bci+. ;; (see below) ;; We can apply an ordinarily (numerical) lisp function of a single real to real to ;; a bci function this way. We create a new bci fun. (defun bcifuncall(f s) ;; create h=f.s so that (funcall h x) is like (funcall f(funcall s x)). (let ((s-fa (funcall s nil))) (let ((fa (map 'simple-array f s-fa)) ;; compute a new set of points. (g (glp (length s-fa))) (w (wta (length s-fa)))) #'(lambda(x)(bcirun fa x g w ))))) ;; here is the generalization to a function of 2 arguments, and then ;; we use it to add, multiply, etc (defun bci2argfun(%%f r s) ;; try computing %%f(r,s) (let ((r-fa (funcall r nil)) (s-fa (funcall s nil))) (cond ((= (length r-fa)(length s-fa)) ;; compatible arrays (let ((fa (map 'simple-array %%f r-fa s-fa)) (g (glp (length r-fa))) (w (wta (length r-fa)))) #'(lambda(x)(bcirun fa x g w)))) (t(error "lengths do not match so we cannot use ~s on ~s and ~s as BCI forms" %%f r s))))) ;;then (defun bci+(r s)(bci2argfun #'+ r s)) (defun bci*(r s)(bci2argfun #'* r s)) (defun bci-(r s)(bci2argfun #'- r s)) (defun bci/(r s)(bci2argfun #'/ r s)) (defun bcimax(r s)(bci2argfun #'max r s)) ;; etc etc ;; we could do integration, differentiation, finding zeros. See chebfun project web page for ;; ideas. ;; an example ;; (setf s (bcif #'sin 8)) ;; make s a function that evaluates sin(x) for x in [-1,1] using ;; 8 point chebyshev / barycentric interpolation. It can be ;; used this way: (funcall s 0.5) computes (sin 0.5), approximately. ;; it gives 0.47942555438789375d0 whereas (sin 0.5d0) ;; gives 0.479425538604203d0. ;; Using (bcif #'sin 16) is a more accurate function: it gives exactly the same for sin(0.5). ;; if you don't like using funcall, you can do this [also] ;; (setf (symbol-function 's) (bcif #'sin 16)) . This ;; means you can write (s 0.5) instead of (funcall s 0.5) ;; (defun h(x) (+ (cos(exp(* 3.0 x)))(sin(exp(* 2.0d0 x))) 12.00d0 (log(+ x 4.0d0)))) ;; (setf hfun (bcif #'h 64)) ;; (defun g(x)(funcall hfun x)) ;; (defun testgh(x)(/ (abs(- (h x)(g x))) (max (abs(h x))(abs (g x))))) #| ;; a heuristic proof that sin^2 + cos^2 =1 (setf s (bcif #'sin 16)) (setf c (bcif #'cos 16)) (setf g (bci+ (bci* c c ) (bci* s s))) ;; g is the function sin^2 + cos^2 (pprint (loop for i from -1 to 1 by 1/10 collect(funcall g i))) ;; the approximation gets bad outside [-1,1]: (loop for i from 1 to 5 by 1/4 do (print (1- (funcall g i)))) ;; to find the coefficients in the Chebyshev series needed to encode ;; some f, do something like this: (setf r (bcif #'sin 128)) ;; best to use a power of 2 so we can use FAST discrete cosine transform (setf rcopy (copy-seq(funcall r nil))) ;; do this because the next step will over-write its argument. (setf r_transform (fct rcopy 128)) ;; fct is fast cosine transform and is defined in file algo749.lisp. ;; The basis for the success of the Chebfun project is that one can ;; judge how many numerical points are needed to represent some ;; function sufficiently accurately by converting the sequence of ;; interpolation points to a Chebyshev series and looking at the ;; coefficients in this transformed space. If, in the Chebyshev space, ;; all coefficients of order > N are consistently at the level of ;; "noise" compared to the larger coefficients, then they may ;; plausibly be dropped entirely. Heuristics based on examining some ;; number of terms can be used to approximate this decision. ;; If the terms do not drop fast enough, we blame the function being ;; approximated: it is ill-behaved in some way. One cure is to split ;; the function into pieces and approximate each piece separately. |# ;; if maxima.. #+GCL (defun maxima::\$bcifun(f n) (let ((gs (intern (symbol-name (gensym "\$BF")) :maxima)) (ff `(lambda(r)(float (mfuncall (quote ,f) r) 1.0d0)))) (setf (symbol-function gs)(bcif ff n)) gs)) #+GCL (defun maxima::\$bcitimes(f g) (let ((gs (intern (symbol-name (gensym "\$BF")) :maxima))) (setf (symbol-function gs)(bci2argfun #'* f g)) gs)) ;; etc etc #| maxima example.. load("bci"); h(x):=sin(10*x)*exp(x); rr:bcifun(h,30); plot2d(['rr(x),h(x)],[x,-1.5,1.5]); plot2d('rr(x)-h(x),[x,-1,1]); s3:bcifun(lambda([x],x^3),10); s6: bcitimes(s3,s3); plot2d(['s6(x),x^6],[x,-1,1]); |#