chebyshev.lisp 5.55 KB
Newer Older
1 2
;; Chebyshev Approximations
;; Liam Healy Sat Nov 17 2007 - 20:36
3
;; Time-stamp: <2012-01-13 12:01:09EST chebyshev.lisp>
4
;;
5
;; Copyright 2007, 2008, 2009, 2010, 2011 Liam M. Healy
6 7 8 9 10 11 12 13 14 15 16 17 18 19
;; Distributed under the terms of the GNU General Public License
;;
;; This program is free software: you can redistribute it and/or modify
;; it under the terms of the GNU General Public License as published by
;; the Free Software Foundation, either version 3 of the License, or
;; (at your option) any later version.
;;
;; This program is distributed in the hope that it will be useful,
;; but WITHOUT ANY WARRANTY; without even the implied warranty of
;; MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
;; GNU General Public License for more details.
;;
;; You should have received a copy of the GNU General Public License
;; along with this program.  If not, see <http://www.gnu.org/licenses/>.
liam committed
20 21 22

(in-package :gsl)

23 24
;;; /usr/include/gsl/gsl_chebyshev.h

25 26 27 28
;;;;****************************************************************************
;;;; Creation and calculation of Chebyshev series
;;;;****************************************************************************

29
(defmobject chebyshev "gsl_cheb"
30
  ((order :sizet))
31 32
  "Chebyshev series"
  :documentation			; FDL
33
  "Make a Chebyshev series of specified order."
34
  :callbacks
35
  (callback (:struct fnstruct) nil (function :double (:input :double) :slug))
36 37
  :initialize-suffix "init"
  :initialize-args
38 39 40
  ((callback :pointer) (lower-limit :double) (upper-limit :double))
  :singular (function))

41 42 43 44
(defmfun order ((object chebyshev))
  "gsl_cheb_order"
  (((mpointer object) :pointer))
  :definition :method
45
  :documentation "The order of Chebyshev series."
46
  :c-return :sizet
47 48
  :gsl-version (1 12))

Liam Healy committed
49
(defmfun size ((chebyshev chebyshev))
50 51
  "gsl_cheb_size"
  (((mpointer chebyshev) :pointer))
Liam Healy committed
52
  :definition :method
53
  :documentation "The length of the Chebyshev coefficient array."
54
  :c-return :sizet
55 56 57 58 59
  :gsl-version (1 12))

(defmfun coefficients (chebyshev)
  "gsl_cheb_coeffs"
  (((mpointer chebyshev) :pointer))
60
  :c-return (crtn :pointer)
61 62
  :documentation
  "The Chebyshev coefficient array as a CL array (foreign-friendly)."
63 64
  :return ((grid:make-foreign-array-from-pointer
	    crtn (list (size chebyshev)) 'double-float t))
65 66
  :gsl-version (1 12))

liam committed
67 68 69 70 71 72
;;;;****************************************************************************
;;;; 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
73
;;; desired) in the form of #'evaluate.
liam committed
74

75
(defmfun evaluate ((object chebyshev) x &key order)
Liam Healy committed
76
  ("gsl_cheb_eval" "gsl_cheb_eval_n")
77
  ((((mpointer object) :pointer) (x :double))
78
   (((mpointer object) :pointer) (order :sizet) (x :double)))
79
  :definition :method
80
  :callback-object object
liam committed
81
  :c-return :double
liam committed
82
  :documentation			; FDL
Liam Healy committed
83 84
  "Evaluate the Chebyshev series at a point x.  If order is supplied,
  evaluate to at most the given order.")
liam committed
85

Liam Healy committed
86 87
(defmfun evaluate-chebyshev-error (chebyshev x &optional order)
  ("gsl_cheb_eval_err" "gsl_cheb_eval_n_err")
88 89
  ((((mpointer chebyshev) :pointer) (x :double)
    (result (:pointer :double)) (abserr (:pointer :double)))
90
   (((mpointer chebyshev) :pointer) (order :sizet) (x :double)
91
    (result (:pointer :double)) (abserr (:pointer :double))))
liam committed
92
  :documentation			; FDL
liam committed
93
  "Evaluate the Chebyshev series at a point x, returning result and
Liam Healy committed
94 95
   an estimate of its absolute error.  If order is supplied,
   evaluate to at most the given order.")
liam committed
96 97 98 99 100

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

liam committed
101
(defmfun derivative-chebyshev (derivative chebyshev)
liam committed
102
  "gsl_cheb_calc_deriv"
103
  (((mpointer derivative) :pointer) ((mpointer chebyshev) :pointer))
liam committed
104
  :documentation			; FDL
liam committed
105 106 107 108
  "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 committed
109
(defmfun integral-chebyshev (integral chebyshev)
liam committed
110
  "gsl_cheb_calc_integ"
111
  (((mpointer integral) :pointer) ((mpointer chebyshev) :pointer))
liam committed
112
  :documentation			; FDL
liam committed
113 114 115 116 117 118 119 120 121 122 123
  "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 committed
124
;;; to keep things sane.
liam committed
125

126
(defun chebyshev-step (x) (if (< x 0.5d0) 0.25d0 0.75d0))
liam committed
127

128
(defun chebyshev-table-example ()
129
  (let ((steps 100))
130
    (let ((cheb (make-chebyshev 40 'chebyshev-step 0.0d0 1.0d0)))
131 132 133 134 135
      (dotimes (i steps)
	(let ((x (coerce (/ i steps) 'double-float)))
	  (format t "~&~a ~a ~a ~a"
		  x
		  (chebyshev-step x)
136 137
		  (evaluate cheb x :order 10)
		  (evaluate cheb x)))))))
138 139 140

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

Liam Healy committed
151 152
;;; Unit test
(save-test chebyshev
liam committed
153
  (chebyshev-point-example 0.55d0))