chebyshev.lisp 4.62 KB
Newer Older
liam's avatar
liam committed
1 2
;; Chebyshev Approximations
;; Liam Healy Sat Nov 17 2007 - 20:36
3
;; Time-stamp: <2009-03-10 22:45:03EDT chebyshev.lisp>
lhealy's avatar
lhealy committed
4
;; $Id$
liam's avatar
liam committed
5 6 7

(in-package :gsl)

8 9
;;; /usr/include/gsl/gsl_chebyshev.h

10 11 12 13
;;;;****************************************************************************
;;;; Creation and calculation of Chebyshev series
;;;;****************************************************************************

14 15
(defmobject chebyshev "gsl_cheb"
  ((order sizet))
Liam Healy's avatar
Liam Healy committed
16 17
  "Chebyshev series"
  :documentation			; FDL
18
  "Make a Chebyshev series of specified order."
19
  ;; Old arguments:
20 21
  :superclasses (callback-included)
  :ci-class-slots (gsl-function nil (function))
22 23 24
  ;; :ci-class-slots (gsl-mfunction-fdf marray (function df fdf))

  ;; New arguments:
25
  :callbacks (callback gsl-function (function))
26 27
  ;:callback-dynamic ((function))

Liam Healy's avatar
Liam Healy committed
28 29
  :initialize-suffix "init"
  :initialize-args
30 31 32
  ((callback :pointer) (lower-limit :double) (upper-limit :double))
  :singular (function))

33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51
(defmfun order ((object chebyshev))
  "gsl_cheb_order"
  (((mpointer object) :pointer))
  :definition :method
  :c-return sizet
  :gsl-version (1 12))

(defmfun size (chebyshev)
  "gsl_cheb_size"
  (((mpointer chebyshev) :pointer))
  :c-return sizet
  :gsl-version (1 12))

(defmfun coefficients (chebyshev)
  "gsl_cheb_coeffs"
  (((mpointer chebyshev) :pointer))
  :c-return sizet
  :gsl-version (1 12))

liam's avatar
liam committed
52 53 54 55 56 57
;;;;****************************************************************************
;;;; 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
58
;;; desired) in the form of #'evaluate.
liam's avatar
liam committed
59

60
(defmfun evaluate ((object chebyshev) x &key order)
Liam Healy's avatar
Liam Healy committed
61
  ("gsl_cheb_eval" "gsl_cheb_eval_n")
62 63 64
  ((((mpointer object) :pointer) (x :double))
   (((mpointer object) :pointer) (order sizet) (x :double)))
  :definition :method
liam's avatar
liam committed
65
  :c-return :double
liam's avatar
liam committed
66
  :documentation			; FDL
Liam Healy's avatar
Liam Healy committed
67 68
  "Evaluate the Chebyshev series at a point x.  If order is supplied,
  evaluate to at most the given order.")
liam's avatar
liam committed
69

Liam Healy's avatar
Liam Healy committed
70 71
(defmfun evaluate-chebyshev-error (chebyshev x &optional order)
  ("gsl_cheb_eval_err" "gsl_cheb_eval_n_err")
72 73 74 75
  ((((mpointer chebyshev) :pointer) (x :double) (result :double)
    (abserr :double))
   (((mpointer chebyshev) :pointer) (order sizet) (x :double)
    (result :double) (abserr :double)))
liam's avatar
liam committed
76
  :documentation			; FDL
liam's avatar
liam committed
77
  "Evaluate the Chebyshev series at a point x, returning result and
Liam Healy's avatar
Liam Healy committed
78 79
   an estimate of its absolute error.  If order is supplied,
   evaluate to at most the given order.")
liam's avatar
liam committed
80 81 82 83 84

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

liam's avatar
liam committed
85
(defmfun derivative-chebyshev (derivative chebyshev)
liam's avatar
liam committed
86
  "gsl_cheb_calc_deriv"
87
  (((mpointer derivative) :pointer) ((mpointer chebyshev) :pointer))
liam's avatar
liam committed
88
  :documentation			; FDL
liam's avatar
liam committed
89 90 91 92
  "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.")

liam's avatar
liam committed
93
(defmfun integral-chebyshev (integral chebyshev)
liam's avatar
liam committed
94
  "gsl_cheb_calc_integ"
95
  (((mpointer integral) :pointer) ((mpointer chebyshev) :pointer))
liam's avatar
liam committed
96
  :documentation			; FDL
liam's avatar
liam committed
97 98 99 100 101 102 103 104 105 106 107
  "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
108
;;; to keep things sane.
liam's avatar
liam committed
109

110
(defun chebyshev-step (x) (if (< x 0.5d0) 0.25d0 0.75d0))
liam's avatar
liam committed
111

112
(make-callbacks single-function chebyshev-step)
113

114
(defun chebyshev-table-example ()
115
  (let ((steps 100))
116
    (let ((cheb (make-chebyshev 40 'chebyshev-step 0.0d0 1.0d0)))
117 118 119 120 121
      (dotimes (i steps)
	(let ((x (coerce (/ i steps) 'double-float)))
	  (format t "~&~a ~a ~a ~a"
		  x
		  (chebyshev-step x)
122 123
		  (evaluate cheb x :order 10)
		  (evaluate cheb x)))))))
124 125 126

(defun chebyshev-point-example (x)
  (check-type x double-float)
127
  (let ((cheb (make-chebyshev 40 'chebyshev-step 0.0d0 1.0d0))
128 129
	(deriv (make-chebyshev 40))
	(integ (make-chebyshev 40)))
130 131
    (derivative-chebyshev deriv cheb)
    (integral-chebyshev integ cheb)
132
    (list
133 134 135
     (evaluate cheb x)
     (evaluate deriv x)
     (evaluate integ x))))
136

Liam Healy's avatar
Liam Healy committed
137 138
;;; Unit test
(save-test chebyshev
liam's avatar
liam committed
139
  (chebyshev-point-example 0.55d0))