Skip to content
Snippets Groups Projects
chebyshev.lisp 5.55 KiB
Newer Older
;; Chebyshev Approximations
;; Liam Healy Sat Nov 17 2007 - 20:36
;; Time-stamp: <2012-01-13 12:01:09EST chebyshev.lisp>
;; Copyright 2007, 2008, 2009, 2010, 2011 Liam M. Healy
;; 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's avatar
liam committed

(in-package :gsl)

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

(defmobject chebyshev "gsl_cheb"
  "Chebyshev series"
  :documentation			; FDL
  "Make a Chebyshev series of specified order."
  (callback (:struct fnstruct) nil (function :double (:input :double) :slug))
  :initialize-suffix "init"
  :initialize-args
  ((callback :pointer) (lower-limit :double) (upper-limit :double))
  :singular (function))

(defmfun order ((object chebyshev))
  "gsl_cheb_order"
  (((mpointer object) :pointer))
  :definition :method
  :documentation "The order of Chebyshev series."
  :gsl-version (1 12))

Liam Healy's avatar
Liam Healy committed
(defmfun size ((chebyshev chebyshev))
  "gsl_cheb_size"
  (((mpointer chebyshev) :pointer))
Liam Healy's avatar
Liam Healy committed
  :definition :method
  :documentation "The length of the Chebyshev coefficient array."
  :gsl-version (1 12))

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

liam's avatar
liam committed
;;;;****************************************************************************
;;;; 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.
liam's avatar
liam committed

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

Liam Healy's avatar
Liam Healy committed
(defmfun evaluate-chebyshev-error (chebyshev x &optional order)
  ("gsl_cheb_eval_err" "gsl_cheb_eval_n_err")
  ((((mpointer chebyshev) :pointer) (x :double)
    (result (:pointer :double)) (abserr (:pointer :double)))
   (((mpointer chebyshev) :pointer) (order :sizet) (x :double)
    (result (:pointer :double)) (abserr (:pointer :double))))
liam's avatar
liam committed
  :documentation			; FDL
liam's avatar
liam committed
  "Evaluate the Chebyshev series at a point x, returning result and
Liam Healy's avatar
Liam Healy committed
   an estimate of its absolute error.  If order is supplied,
   evaluate to at most the given order.")
liam's avatar
liam committed

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

liam's avatar
liam committed
(defmfun derivative-chebyshev (derivative chebyshev)
liam's avatar
liam committed
  "gsl_cheb_calc_deriv"
  (((mpointer derivative) :pointer) ((mpointer chebyshev) :pointer))
liam's avatar
liam committed
  :documentation			; FDL
liam's avatar
liam committed
  "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
(defmfun integral-chebyshev (integral chebyshev)
liam's avatar
liam committed
  "gsl_cheb_calc_integ"
  (((mpointer integral) :pointer) ((mpointer chebyshev) :pointer))
liam's avatar
liam committed
  :documentation			; FDL
liam's avatar
liam committed
  "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
;;; to keep things sane.
(defun chebyshev-step (x) (if (< x 0.5d0) 0.25d0 0.75d0))
(defun chebyshev-table-example ()
    (let ((cheb (make-chebyshev 40 'chebyshev-step 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 cheb x :order 10)
		  (evaluate cheb x)))))))

(defun chebyshev-point-example (x)
  (check-type x double-float)
  (let ((cheb (make-chebyshev 40 'chebyshev-step 0.0d0 1.0d0))
	(deriv (make-chebyshev 40))
	(integ (make-chebyshev 40)))
    (derivative-chebyshev deriv cheb)
    (integral-chebyshev integ cheb)
     (evaluate cheb x)
     (evaluate deriv x)
     (evaluate integ x))))
Liam Healy's avatar
Liam Healy committed
;;; Unit test
(save-test chebyshev
liam's avatar
liam committed
  (chebyshev-point-example 0.55d0))