spline-example.lisp 1.61 KB
Newer Older
1
2
;; Example spline
;; Liam Healy, Sat Nov 10 2007 - 21:18
Liam Healy's avatar
Liam Healy committed
3
;; Time-stamp: <2009-02-16 09:46:57EST spline-example.lisp>
lhealy's avatar
lhealy committed
4
;; $Id$
5
6
7

(in-package :gsl)

Liam Healy's avatar
Liam Healy committed
8
(defun spline-example (&optional (step 0.01d0))
9
  "The first example in Sec. 26.7 of the GSL manual."
Liam Healy's avatar
Liam Healy committed
10
11
12
13
14
15
16
17
18
19
20
21
22
  (let* ((acc (make-acceleration))
	 (xarr
	  (make-marray
	   'double-float
	   :initial-contents
	   (loop for i from 0.0d0 below 10.0d0
	      collect (+ i (* 0.5d0 (sin i))))))
	 (yarr
	  (make-marray
	   'double-float
	   :initial-contents
	   (loop for i from 0.0d0 below 10.0d0
	      collect (+ i (cos (expt i 2))))))
Liam Healy's avatar
Liam Healy committed
23
	 (spline (make-spline +cubic-spline-interpolation+ xarr yarr)))
Liam Healy's avatar
Liam Healy committed
24
    (loop for xi from (maref xarr 0) below (maref xarr 9) by step
25
       collect (list xi (evaluate spline xi :acceleration acc)))))
Liam Healy's avatar
Liam Healy committed
26
27

(save-test interpolation (spline-example 0.1d0))
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48

(defun evaluate-integral-example (&optional (intervals 4))
  "Evaluate integral of sin(x) in interval 0-pi.  sin(x) is tabulated
over a 0-2pi interval and interpolated with
+periodic-cubic-spline-interpolation+"
  (let* ((nodes (1+ intervals))
	 (max-node (1- nodes))
	 (xarr 
	  (loop 
	     with step = (/ (* 2.0 pi) intervals)
	     for i from 0 upto max-node
	     collect (* i step)))
	 (xmarr (make-marray 'double-float :initial-contents xarr))
	 ;; cannot use (loop for x on (cl-array xmarr)...) -- c function gives error
	 (ymarr 
	  (make-marray 'double-float :initial-contents
		       (loop for x in xarr
			  collect (sin x))))
	 (acc (make-acceleration))
	 (spline (make-spline +periodic-cubic-spline-interpolation+ xmarr ymarr)))
  (evaluate-integral spline 0d0 (coerce pi 'double-float) :acceleration acc)))