chebyshev.lisp 4.58 KB
Newer Older
liam's avatar
liam committed
1 2
;; Chebyshev Approximations
;; Liam Healy Sat Nov 17 2007 - 20:36
Liam Healy's avatar
Liam Healy committed
3
;; Time-stamp: <2008-10-25 11:33:06EDT chebyshev.lisp>
lhealy's avatar
lhealy 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
(defmfun allocate-chebyshev (order)
liam's avatar
liam committed
16
  "gsl_cheb_alloc"
Liam Healy's avatar
Liam Healy committed
17
  ((order sizet))
liam's avatar
liam committed
18
  :c-return :pointer
19
  :export nil
20
  :index (letm chebyshev)
liam's avatar
liam committed
21
  :documentation			; FDL
liam's avatar
liam committed
22 23 24
  "Allocate a Chebyshev series of specified order
   and return a pointer to it.")

liam's avatar
liam committed
25
(defmfun free-chebyshev (chebyshev)
liam's avatar
liam committed
26 27 28
  "gsl_cheb_free"
  ((chebyshev :pointer))
  :c-return :void
29
  :export nil
30
  :index (letm chebyshev)
liam's avatar
liam committed
31
  :documentation			; FDL
liam's avatar
liam committed
32 33
  "Free a previously allocated Chebyshev series.")

liam's avatar
liam committed
34
(defmfun initialize-chebyshev (chebyshev function lower-limit upper-limit)
liam's avatar
liam committed
35 36 37
  "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
  :documentation			; FDL
liam's avatar
liam committed
41 42
  "Compute the Chebyshev approximation for the function over the range
   (lower-limit, upper-limit) to the previously specified order.  The
liam's avatar
liam committed
43 44
   computation of the Chebyshev approximation is an O(n^2)
   process, and requires n function evaluations.")
liam's avatar
liam committed
45 46 47 48 49 50 51 52 53

;;;;****************************************************************************
;;;; 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.

Liam Healy's avatar
Liam Healy committed
54 55 56 57
(defmfun evaluate-chebyshev (chebyshev x &optional order)
  ("gsl_cheb_eval" "gsl_cheb_eval_n")
  (((chebyshev :pointer) (x :double))
  ((chebyshev :pointer) (order sizet) (x :double)))
liam's avatar
liam committed
58
  :c-return :double
liam's avatar
liam committed
59
  :documentation			; FDL
Liam Healy's avatar
Liam Healy committed
60 61
  "Evaluate the Chebyshev series at a point x.  If order is supplied,
  evaluate to at most the given order.")
liam's avatar
liam committed
62

Liam Healy's avatar
Liam Healy committed
63 64 65 66
(defmfun evaluate-chebyshev-error (chebyshev x &optional order)
  ("gsl_cheb_eval_err" "gsl_cheb_eval_n_err")
  (((chebyshev :pointer) (x :double) (result :double) (abserr :double))
   ((chebyshev :pointer) (order sizet) (x :double) (result :double) (abserr :double)))
liam's avatar
liam committed
67
  :documentation			; FDL
liam's avatar
liam committed
68
  "Evaluate the Chebyshev series at a point x, returning result and
Liam Healy's avatar
Liam Healy committed
69 70
   an estimate of its absolute error.  If order is supplied,
   evaluate to at most the given order.")
liam's avatar
liam committed
71 72 73 74 75

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

liam's avatar
liam committed
76
(defmfun derivative-chebyshev (derivative chebyshev)
liam's avatar
liam committed
77 78
  "gsl_cheb_calc_deriv"
  ((derivative :pointer) (chebyshev :pointer))
liam's avatar
liam committed
79
  :documentation			; FDL
liam's avatar
liam committed
80 81 82 83
  "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
84
(defmfun integral-chebyshev (integral chebyshev)
liam's avatar
liam committed
85 86
  "gsl_cheb_calc_integ"
  ((integral :pointer) (chebyshev :pointer))
liam's avatar
liam committed
87
  :documentation			; FDL
liam's avatar
liam committed
88 89 90 91 92 93 94 95 96 97 98
  "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
99
;;; to keep things sane.
liam's avatar
liam committed
100

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

103
(defun chebyshev-table-example ()
104
  (let ((steps 100))
105
    (letm ((cheb (chebyshev 40 chebyshev-step 0.0d0 1.0d0)))
106 107 108 109 110 111 112
      (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)))))))
113 114 115

(defun chebyshev-point-example (x)
  (check-type x double-float)
116 117 118
  (letm ((cheb (chebyshev 40 chebyshev-step 0.0d0 1.0d0))
	 (deriv (chebyshev 40))
	 (integ (chebyshev 40)))
119 120
    (derivative-chebyshev deriv cheb)
    (integral-chebyshev integ cheb)
121 122 123 124
    (list
     (evaluate-chebyshev cheb x)
     (evaluate-chebyshev deriv x)
     (evaluate-chebyshev integ x))))
125

Liam Healy's avatar
Liam Healy committed
126 127
;;; Unit test
(save-test chebyshev
liam's avatar
liam committed
128
  (chebyshev-point-example 0.55d0))