chebyshev.lisp 5.11 KB
Newer Older
liam's avatar
liam committed
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139
;********************************************************
; file:        chebyshev.lisp                            
; description: Chebyshev Approximations                  
; date:        Sat Nov 17 2007 - 20:36                   
; author:      Liam Healy                                
; modified:    Sat Nov 17 2007 - 22:35
;********************************************************
;;; $Id: $

(in-package :gsl)

;;; Needs with- macro to allocate, initialize, and free object.

;;;;****************************************************************************
;;;; Creation and calculation of Chebyshev series
;;;;****************************************************************************

(defun-gsl allocate-chebyshev (order)
  "gsl_cheb_alloc"
  ((order :size))
  :c-return :pointer
  :documentation
  "Allocate a Chebyshev series of specified order
   and return a pointer to it.")

(defun-gsl free-chebyshev (chebyshev)
  "gsl_cheb_free"
  ((chebyshev :pointer))
  :c-return :void
  :documentation
  "Free a previously allocated Chebyshev series.")

(defun-gsl initialize-chebyshev (chebyshev function lower-limit upper-limit)
  "gsl_cheb_init"
  ((chebyshev :pointer) (function :pointer)
   (lower-limit :double) (upper-limit :double))
  :documentation
  "Compute the Chebyshev approximation for the function over the range
   (lower-limit, upper-limit) to the previously specified order.  The
   computation of the Chebyshev approximation is an @math{O(n^2)}
   process, and requires @math{n} function evaluations.")

;;;;****************************************************************************
;;;; Chebyshev series evaluation
;;;;****************************************************************************

;;; The functions that don't return are defined, but it is recommended
;;; to use the functions that do return error (and ignore it if
;;; desired) in the form of #'evaluate-chebyshev.

(defun-gsl evaluate-chebyshev-noerror (chebyshev x)
  "gsl_cheb_eval"
  ((chebyshev :pointer) (x :double))
  :c-return :double
  :index evaluate-chebyshev
  :export nil
  :documentation
  "Evaluate the Chebyshev series at a point x.")

(defun-gsl evaluate-chebyshev-noerror-order (chebyshev x order)
  "gsl_cheb_eval_n"
  ((chebyshev :pointer) (order :size) (x :double))
  :c-return :double
  :index evaluate-chebyshev
  :export nil
  :documentation
  "Evaluate the Chebyshev series at a point x to at most the given order.")

(defun-gsl evaluate-chebyshev-full (chebyshev x)
  "gsl_cheb_eval_err"
  ((chebyshev :pointer) (x :double) (result :double) (abserr :double))
  :index evaluate-chebyshev
  :export nil
  :documentation
  "Evaluate the Chebyshev series at a point x, returning result and
   an estimate of its absolute error.")

(defun-gsl evaluate-chebyshev-order (chebyshev x order)
  "gsl_cheb_eval_n_err"
  ((chebyshev :pointer) (order :size) (x :double) (result :double) (abserr :double))
  :index evaluate-chebyshev
  :export nil
  :documentation
  "Evaluate the Chebyshev series at a point x to at most the given order,
   returning result and an estimate of its absolute error.")

(export 'evaluate-chebyshev)
(defun-optionals evaluate-chebyshev (chebyshev x &optional order)
  -full -order
  "Evaluate the Chebyshev series at a point x to at most the given order,
   returning result and an estimate of its absolute error.")

;;;;****************************************************************************
;;;; Derivatives and integrals
;;;;****************************************************************************

(defun-gsl derivative-chebyshev (derivative chebyshev)
  "gsl_cheb_calc_deriv"
  ((derivative :pointer) (chebyshev :pointer))
  :documentation
  "Compute the derivative of the Chebyshev series, storing
   the derivative coefficients in the previously allocated series.
   The two series must have been allocated with the same order.")

(defun-gsl integral-chebyshev (integral chebyshev)
  "gsl_cheb_calc_integ"
  ((integral :pointer) (chebyshev :pointer))
  :documentation
  "Compute the integral of the Chebyshev series, storing
   the integral coefficients in the previously allocated series.
   The two series must have been allocated with the same order.
   The lower limit of the integration is taken to be the left hand
   end of the range lower-limit.")

;;;;****************************************************************************
;;;; Example
;;;;****************************************************************************

;;; From Chap. 28.5, except I have set steps = 100 instead of 10000
;;; to things sane.

;;; Calling a callback from CL is not possible, so we define a
;;; separate CL function to show what the answer really is.
(defun chebyshev-step (x) (if (< x 0.5) 0.25 0.75))
(def-gsl-function chebyshev-step-c x (chebyshev-step x))

(defun chebyshev-example ()
  (let ((cheb (allocate-chebyshev 40))
	(steps 100))
    (with-integration-function (step-fn 'chebyshev-step-c)
      (initialize-chebyshev cheb step-fn 0.0d0 1.0d0)
      (dotimes (i steps)
	(let ((x (coerce (/ i steps) 'double-float)))
	  (format t "~&~a ~a ~a ~a"
		  x
		  (chebyshev-step x)
		  (evaluate-chebyshev cheb x 10)
		  (evaluate-chebyshev cheb x))))
      (free-chebyshev cheb))))