chebyshev.lisp 5.3 KB
Newer Older
liam's avatar
liam committed
1 2
;; Chebyshev Approximations
;; Liam Healy Sat Nov 17 2007 - 20:36
3
;; Time-stamp: <2008-02-03 13:22:58EST chebyshev.lisp>
liam's avatar
liam committed
4
;; $Id: $
liam's avatar
liam committed
5 6 7 8 9 10 11

(in-package :gsl)

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

12
(defgo-s (chebyshev order function lower-limit upper-limit)
13
	 allocate-chebyshev free-chebyshev initialize-chebyshev)
14

liam's avatar
liam committed
15 16 17 18
(defun-gsl allocate-chebyshev (order)
  "gsl_cheb_alloc"
  ((order :size))
  :c-return :pointer
19
  :export nil
20
  :index (letm chebyshev)
liam's avatar
liam committed
21 22 23 24 25 26 27 28
  :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
29
  :export nil
30
  :index (letm chebyshev)
liam's avatar
liam committed
31 32 33 34 35 36 37
  :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))
38
  :export nil
39
  :index (letm chebyshev)
liam's avatar
liam committed
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
  :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
liam's avatar
liam committed
123
;;; to keep things sane.
liam's avatar
liam committed
124

liam's avatar
liam committed
125
(defun-single chebyshev-step (x) (if (< x 0.5d0) 0.25d0 0.75d0))
liam's avatar
liam committed
126

127
(defun chebyshev-table-example ()
128
  (let ((steps 100))
129
    (letm ((cheb (chebyshev 40 chebyshev-step 0.0d0 1.0d0)))
130 131 132 133 134 135 136
      (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)))))))
137 138 139

(defun chebyshev-point-example (x)
  (check-type x double-float)
140 141 142
  (letm ((cheb (chebyshev 40 chebyshev-step 0.0d0 1.0d0))
	 (deriv (chebyshev 40))
	 (integ (chebyshev 40)))
143 144
    (derivative-chebyshev deriv cheb)
    (integral-chebyshev integ cheb)
145 146 147 148
    (list
     (evaluate-chebyshev cheb x)
     (evaluate-chebyshev deriv x)
     (evaluate-chebyshev integ x))))
149 150 151 152 153

(lisp-unit:define-test chebyshev
  (lisp-unit:assert-equal
   '("0.715920990169d+00" "-0.150199666581d+01" "0.172397194040d+00")
   (lisp-unit:fp-sequence (chebyshev-point-example 0.55d0))))