diff --git a/documentation/index.html b/documentation/index.html index 2215b1e1369885e32d7d4f0dc7f6a30a45a1e758..13d0225bdcbbdeb16662ba3ef6f1d679fa1ed51a 100644 --- a/documentation/index.html +++ b/documentation/index.html @@ -340,7 +340,8 @@ random-number-generator quasi-random-number-generator discrete-random polynomial-complex-workspace integration-workspace eigen-symm eigen-symmv eigen-herm eigen-hermv monte-carlo-plain monte-carlo-miser monte-carlo-vegas -ode-stepper ode-evolution +ode-stepper ode-evolution standard-control y-control +yp-control scaled-control </pre> <p> An instance may be created with a function whose name is "make-" @@ -387,7 +388,7 @@ and arrays used internally or for function return. <!-- Created: Feb 25 2005 --> <!-- hhmts start --> <small> -Time-stamp: <2009-01-24 21:04:19EST index.html> +Time-stamp: <2009-01-25 21:38:05EST index.html> </small> <!-- hhmts end --> </div> diff --git a/init/callback.lisp b/init/callback.lisp index b90d6a8cf7455882cb6d0cc23dd7c7a98c3df64e..ad40f45f55e0d62b0e79cbc7402838337f5a514a 100644 --- a/init/callback.lisp +++ b/init/callback.lisp @@ -1,6 +1,6 @@ ;; Foreign callback functions. ;; Liam Healy -;; Time-stamp: <2009-01-24 19:31:51EST callback.lisp> +;; Time-stamp: <2009-01-25 19:21:02EST callback.lisp> ;; $Id$ (in-package :gsl) @@ -31,8 +31,7 @@ ;;; desirable to read and set marrays, in which case :pointer is the ;;; right specification. -(export - '(make-single-function undef-cbstruct with-c-double with-c-doubles)) +(export '(make-single-function undef-cbstruct)) ;;;;**************************************************************************** ;;;; Setting slots @@ -270,33 +269,3 @@ structure `((defcbstruct ,cbsymb ,structure ,(if dim `((dimensions ,dim)))))))))) - -;;;;**************************************************************************** -;;;; Vector of doubles -;;;;**************************************************************************** - -;;; Unfortunately, because vectors must be copied between languages -;;; (even with vector-sap in callbacks, unless vector-sap can be -;;; setfed), there is no way to provide a function of a vector and -;;; have it work in both languages. As a consolation the macro -;;; #'with-c-double is provided to make things easier. - -(defmacro with-c-double - ((c-vector &rest element-names) &body body) - "Provide named access to each element of a C array of doubles, for either - reading or setting." - `(symbol-macrolet - ,(loop for i from 0 for a in element-names - collect `(,a (dcref ,c-vector ,i))) - ,@body)) - -(defmacro with-c-doubles - ((&rest cvector-names) &body body) - "Provide named access to each element of a C array of doubles, for either - reading or setting." - `(symbol-macrolet - ,(loop for (c-vector . element-names) in cvector-names - append - (loop for i from 0 for a in element-names - collect `(,a (dcref ,c-vector ,i)))) - ,@body)) diff --git a/ordinary-differential-equations/evolution.lisp b/ordinary-differential-equations/evolution.lisp index e764bd7fac7ff507631bbfefa5119a08f997dd62..2d8832b5acb413ec434071bd0cf022fdf6b98824 100644 --- a/ordinary-differential-equations/evolution.lisp +++ b/ordinary-differential-equations/evolution.lisp @@ -1,6 +1,6 @@ ;; Evolution functions for ODE integration. ;; Liam Healy, Sun Sep 30 2007 - 14:31 -;; Time-stamp: <2009-01-25 17:11:41EST evolution.lisp> +;; Time-stamp: <2009-01-25 17:52:45EST evolution.lisp> ;; $Id$ (in-package :gsl) @@ -18,8 +18,9 @@ "gsl_odeiv_evolve_apply" (((mpointer evolve) :pointer) ((mpointer control) :pointer) ((mpointer step) :pointer) - (dydt :pointer) (time :pointer) (max-time :double) - (step-size :pointer) (y :pointer)) + (dydt :pointer) ((c-pointer time) :pointer) + (max-time :double) + ((c-pointer step-size) :pointer) ((c-pointer y) :pointer)) :documentation ; FDL "Advance the system (e, dydt) from time and position y using the stepping function step. diff --git a/ordinary-differential-equations/ode-example.lisp b/ordinary-differential-equations/ode-example.lisp index a4fa2871c9f36e6f27905b55b399101b081c0204..e9724c10d9394fc028a174e95bf39927242ee2a5 100644 --- a/ordinary-differential-equations/ode-example.lisp +++ b/ordinary-differential-equations/ode-example.lisp @@ -1,6 +1,6 @@ ;; Example ODE ;; Liam Healy Sat Sep 29 2007 - 17:49 -;; Time-stamp: <2009-01-24 14:34:20EST ode-example.lisp> +;; Time-stamp: <2009-01-25 19:15:39EST ode-example.lisp> ;; $Id$ ;;; van der Pol as given in Section 25.5 of the GSL manual. To @@ -26,22 +26,34 @@ (defparameter *max-iter* 2000) -(defun integrate-vanderpol (max-time &optional (step-size 1.0d-6) (print-steps t)) +(defun integrate-vanderpol + (max-time &optional (step-size 1.0d-6) (stepper *step-rk8pd*) (print-steps t)) "Integrate the van der Pol oscillator as given in Section 25.5 of the GSL manual. To reproduce that example, (integrate-vanderpol 100.0d0)." - (let ((mu 10.0d0) (time 0.0d0) (iter 0)) + (let ((mu 10.0d0) (initial-time 0.0d0) (iter 0)) (declare (special mu)) - (with-ode-integration (time step-size (dependent dep0 dep1) 2) - (setf dep0 1.0d0 dep1 0.0d0) - (loop (when (or (>= (dcref time) max-time) (> iter *max-iter*)) - (return (values (dcref time) dep0 dep1))) - (apply-evolution - evolve control stepper - (make-ode-functions vanderpol vanderpol-jacobian 2) - time max-time step-size dependent) + (with-ode-integration + ((make-ode-functions vanderpol vanderpol-jacobian 2) + time step max-time (dep0 dep1) 2 stepper) + (setf dep0 1.0d0 dep1 0.0d0 step step-size time initial-time) + (loop + (when (or (>= time max-time) (> iter *max-iter*)) + (return (values iter time dep0 dep1))) + make-next-step (incf iter) (when print-steps - (format t "~12,6f~10t~12,6f~24t~12,6f~&" - (dcref time) dep0 dep1)))))) - -(save-test ode (integrate-vanderpol 1.0d0 1.d-4 nil)) + (format t "~12,6f~10t~12,6f~24t~12,6f~&" time dep0 dep1)))))) + +(save-test + ode + (integrate-vanderpol 1.0d0 1.d-4 *step-rk2* nil) + (integrate-vanderpol 1.0d0 1.d-4 *step-rk4* nil) + (integrate-vanderpol 1.0d0 1.d-4 *step-rkf45* nil) + (integrate-vanderpol 1.0d0 1.d-4 *step-rkck* nil) + (integrate-vanderpol 1.0d0 1.d-4 *step-rk8pd* nil) + (integrate-vanderpol 1.0d0 1.d-4 *step-rk2imp* nil) + (integrate-vanderpol 1.0d0 1.d-4 *step-rk4imp* nil) + (integrate-vanderpol 1.0d0 1.d-4 *step-bsimp* nil) + (let ((*max-iter* 12000)) + (integrate-vanderpol 1.0d0 1.d-4 *step-gear1* nil)) + (integrate-vanderpol 1.0d0 1.d-4 *step-gear2* nil)) diff --git a/ordinary-differential-equations/ode-system.lisp b/ordinary-differential-equations/ode-system.lisp index ba5615e7d1ba583dfe1b5a92ae056c91ef168f64..38b422a3f15cb48d9659319d5fcba51e11507304 100644 --- a/ordinary-differential-equations/ode-system.lisp +++ b/ordinary-differential-equations/ode-system.lisp @@ -1,6 +1,6 @@ ;; ODE system setup ;; Liam Healy, Sun Apr 15 2007 - 14:19 -;; Time-stamp: <2009-01-25 17:10:32EST ode-system.lisp> +;; Time-stamp: <2009-01-25 19:01:57EST ode-system.lisp> ;; $Id$ (in-package :gsl) @@ -42,28 +42,30 @@ ((dimension ,dimension)))))) (defmacro with-ode-integration - ((time step-size dependent dimensions &optional (stepper '*step-rk8pd*) - (absolute-error 1.0d-6) (relative-error 0.0d0)) + ((function time step-size max-time + dependent dimensions &optional (stepper '*step-rk8pd*) + (absolute-error 1.0d-6) (relative-error 0.0d0)) &body body) - "Environment for integration of ordinary differential equations. - The variables time and step-size will become C doubles in the body; - to convert back, use #'dcref The dependent variable may - be specified as a list being the same as the first argument to - with-c-double." - (let ((ctime (make-symbol "CTIME")) + "Environment for integration of ordinary differential equations." + (let ((dep (make-symbol "DEP")) + (ctime (make-symbol "CTIME")) (cstep (make-symbol "CSTEP"))) - `(let ((stepper (make-ode-stepper ,stepper ,dimensions)) + `(let ((stepperobj (make-ode-stepper ,stepper ,dimensions)) (control (make-y-control ,absolute-error ,relative-error)) - (evolve (make-ode-evolution ,dimensions))) - (cffi:with-foreign-objects - ((,(if (listp dependent) (first dependent) dependent) - :double ,dimensions) (,ctime :double) (,cstep :double)) - (setf - (dcref ,cstep) ,step-size - (dcref ,ctime) ,time) - ,(if (listp dependent) - `(with-c-double ,dependent - (symbol-macrolet ((,time ,ctime) (,step-size ,cstep)) - ,@body)) - `(symbol-macrolet ((,time ,ctime) (,step-size ,cstep)) - ,@body)))))) + (evolve (make-ode-evolution ,dimensions)) + (,dep + (make-marray 'double-float :dimensions ,dimensions)) + (,ctime (make-marray 'double-float :dimensions 1)) + (,cstep (make-marray 'double-float :dimensions 1))) + (symbol-macrolet + ((,time (maref ,ctime 0)) + (,step-size (maref ,cstep 0)) + ,@(loop for symb in dependent + for i from 0 + collect `(,symb (maref ,dep ,i))) + (make-next-step + (apply-evolution + evolve control stepperobj + ,function ,ctime + ,max-time ,cstep ,dep))) + ,@body)))) diff --git a/tests/ode.lisp b/tests/ode.lisp index 2a760d4f706db69afdd493b685e5d1bd591d1f87..971eeaf81c7c043e1de034fe6254e7c6e6852014 100644 --- a/tests/ode.lisp +++ b/tests/ode.lisp @@ -4,8 +4,55 @@ (LISP-UNIT:DEFINE-TEST ODE (LISP-UNIT::ASSERT-NUMERICAL-EQUAL - (LIST 1.0d0 -1.4568622636249005d0 + (LIST 495 1.0d0 -1.4568657425802234d0 + -11.547345633897558d0) + (MULTIPLE-VALUE-LIST + (INTEGRATE-VANDERPOL 1.0d0 1.d-4 *STEP-RK2* NIL))) + (LISP-UNIT::ASSERT-NUMERICAL-EQUAL + (LIST 40 1.0d0 -1.4568569264026898d0 + -11.547449151779395d0) + (MULTIPLE-VALUE-LIST + (INTEGRATE-VANDERPOL 1.0d0 1.d-4 *STEP-RK4* NIL))) + (LISP-UNIT::ASSERT-NUMERICAL-EQUAL + (LIST 35 1.0d0 -1.456874342553472d0 + -11.547250693698407d0) + (MULTIPLE-VALUE-LIST + (INTEGRATE-VANDERPOL 1.0d0 1.d-4 *STEP-RKF45* NIL))) + (LISP-UNIT::ASSERT-NUMERICAL-EQUAL + (LIST 27 1.0d0 -1.4568588825970334d0 + -11.547432342449643d0) + (MULTIPLE-VALUE-LIST + (INTEGRATE-VANDERPOL 1.0d0 1.d-4 *STEP-RKCK* NIL))) + (LISP-UNIT::ASSERT-NUMERICAL-EQUAL + (LIST 16 1.0d0 -1.4568622636249005d0 -11.547385179410822d0) (MULTIPLE-VALUE-LIST - (INTEGRATE-VANDERPOL 1.0d0 1.d-4 NIL)))) + (INTEGRATE-VANDERPOL 1.0d0 1.d-4 *STEP-RK8PD* NIL))) + (LISP-UNIT::ASSERT-NUMERICAL-EQUAL + (LIST 359 1.0d0 -1.4569507371916566d0 + -11.546406519788615d0) + (MULTIPLE-VALUE-LIST + (INTEGRATE-VANDERPOL 1.0d0 1.d-4 *STEP-RK2IMP* NIL))) + (LISP-UNIT::ASSERT-NUMERICAL-EQUAL + (LIST 40 1.0d0 -1.4568644436344684d0 + -11.547361181933427d0) + (MULTIPLE-VALUE-LIST + (INTEGRATE-VANDERPOL 1.0d0 1.d-4 *STEP-RK4IMP* NIL))) + (LISP-UNIT::ASSERT-NUMERICAL-EQUAL + (LIST 9 1.0d0 -1.4568620806209944d0 + -11.547387321938151d0) + (MULTIPLE-VALUE-LIST + (INTEGRATE-VANDERPOL 1.0d0 1.d-4 *STEP-BSIMP* NIL))) + (LISP-UNIT::ASSERT-NUMERICAL-EQUAL + (LIST 11176 1.0d0 -1.459959741392502d0 + -11.510866371379606d0) + (MULTIPLE-VALUE-LIST + (LET ((*MAX-ITER* 12000)) + (INTEGRATE-VANDERPOL 1.0d0 1.d-4 *STEP-GEAR1* + NIL)))) + (LISP-UNIT::ASSERT-NUMERICAL-EQUAL + (LIST 52 1.0d0 -1.4568645170220367d0 + -11.54736019525012d0) + (MULTIPLE-VALUE-LIST + (INTEGRATE-VANDERPOL 1.0d0 1.d-4 *STEP-GEAR2* NIL))))