diff --git a/init/interface.lisp b/init/interface.lisp index 0fac265f1c2e71fa776bde5dfcdfda1728235099..13a088df4b5c594d0c050154ad9a1e26f03aa66d 100644 --- a/init/interface.lisp +++ b/init/interface.lisp @@ -3,7 +3,7 @@ ; description: Macros to interface GSL functions. ; date: Mon Mar 6 2006 - 22:35 ; author: Liam M. Healy -; modified: Mon Jan 1 2007 - 17:26 +; modified: Sat Sep 29 2007 - 22:43 ;******************************************************** (in-package :gsl) @@ -102,24 +102,35 @@ ;;;;**************************************************************************** (defvar *special-c-return* - '(:error-code :number-of-answers :success-failure :true-false)) + '(:error-code :number-of-answers :success-failure :true-false :enumerate)) -;;; c-arguments List of (symbol c-type). Anything not in arglist will be allocated. ;;; arglist List of CL arguments. +;;; gsl-name Name of the GSL C function, as a quoted string. +;;; c-arguments List of (symbol c-type). Anything not in arglist will be allocated. ;;; c-return, a symbol naming a type, (e.g. :int, :double, :void), ;;; or a list of (symbol type) to name the value, ;;; or :error-code, :number-of-answers, :success-failure, -;;; :true-false +;;; :true-false, :enumerate. If :enumeration is given, +;;; the :enumeration keyword argument will supply the name +;;; of the enumeration. ;;; return, a list of quantities to return. ;;; May be or include :c-return to include the c-return value ;;; or its derivatives. ;;; Default are allocated quantities in c-arguments, or :c-return if none. +;;; type :function or :method +;;; index Name under which this function should be cross-referenced +;;; export Whether to export the symbol. +;;; null-pointer-info Return value if C function returns a null pointer. +;;; documentation +;;; invalidate Invalidate the CL array/matrix cache. +;;; after After method. +;;; enumeration The name of the enumeration return. (defmacro defun-gsl (name arglist gsl-name c-arguments &key (c-return :error-code) (return nil return-supplied-p) (type :function) index (export (not (eq type :method))) - null-pointer-info documentation invalidate after) + null-pointer-info documentation invalidate after enumeration) (let* ((cargs (substitute '(mode sf-mode) :mode c-arguments)) (carg-symbs (remove-if-not #'symbolp @@ -193,6 +204,8 @@ `(,@clret (success-failure ,cret-name))))) (:true-false `((not (zerop ,cret-name)))) + (:enumerate + `((cffi:foreign-enum-keyword ',enumeration ,cret-name))) (t (unless (or (and (eq c-return :error-code) diff --git a/ordinary-differential-equations/control.lisp b/ordinary-differential-equations/control.lisp new file mode 100644 index 0000000000000000000000000000000000000000..6df8bd6c25df646fa9ed96e0e0a56bd3b629e4cc --- /dev/null +++ b/ordinary-differential-equations/control.lisp @@ -0,0 +1,140 @@ +;******************************************************** +; file: control.lisp +; description: Adaptive step-size control +; date: Sat Sep 29 2007 - 18:51 +; author: Liam Healy +; modified: Sat Sep 29 2007 - 22:51 +;******************************************************** +;;; $Id: $ + +(in-package :gsl) + +(defun-gsl new-standard-control (absolute-error relative-error y dydt) + "gsl_odeiv_control_standard_new" + ((absolute-error :double) (relative-error :double) (y :double) (dydt :double)) + :c-return (ptr :pointer) + :return (ptr) + :documentation + "The standard control object is a four parameter heuristic based on + absolute and relative errors @var{eps_absolute} and @var{eps_relative}, and + scaling factors @var{a_y} and @var{a_dydt} for the system state + @math{y(t)} and derivatives @math{y'(t)} respectively. + + The step-size adjustment procedure for this method begins by computing + the desired error level @math{D_i} for each component, + D_i = \epsilon_{absolute} + \epsilon_{relative} * (a_{y} |y_i| + a_{dydt} h |y'_i|) + and comparing it with the observed error @math{E_i = |yerr_i|}. If the + observed error @var{E} exceeds the desired error level @var{D} by more + than 10% for any component then the method reduces the step-size by an + appropriate factor, + h_{new} = h_{old} * S * (E/D)^{-1/q} + where @math{q} is the consistency order of the method (e.g. @math{q=4} for + 4(5) embedded RK), and @math{S} is a safety factor of 0.9. The ratio + @math{E/D} is taken to be the maximum of the ratios @math{E_i/D_i}. + + If the observed error @math{E} is less than 50% of the desired error + level @var{D} for the maximum ratio @math{E_i/D_i} then the algorithm + takes the opportunity to increase the step-size to bring the error in + line with the desired level, + h_{new} = h_{old} * S * (E/D)^{-1/(q+1)} + This encompasses all the standard error scaling methods. To avoid + uncontrolled changes in the stepsize, the overall scaling factor is + limited to the range @math{1/5} to 5.") + +(defun-gsl new-y-control (absolute-error relative-error) + "gsl_odeiv_control_y_new" + ((absolute-error :double) (relative-error :double)) + :c-return (ptr :pointer) + :return (ptr) + :documentation + "Create a new control object which will keep the local + error on each step within an absolute error of @var{eps_absolute} and + relative error of @var{eps_relative} with respect to the solution @math{y_i(t)}. + This is equivalent to the standard control object with @var{a_y}=1 and + @var{a_dydt}=0.") + +(defun-gsl new-yp-control (absolute-error relative-error) + "gsl_odeiv_control_yp_new" + ((absolute-error :double) (relative-error :double)) + :c-return (ptr :pointer) + :return (ptr) + :documentation + "Create a new control object which will keep the local + error on each step within an absolute error of @var{eps_absolute} and + relative error of @var{eps_relative} with respect to the derivatives of the + solution @math{y'_i(t)}. This is equivalent to the standard control + object with @var{a_y}=0 and @var{a_dydt}=1.") + +(defun-gsl new-scaled-control + (absolute-error relative-error y dydt absolute-scale dimension) + "gsl_odeiv_control_scaled_new" + ((absolute-error :double) (relative-error :double) (y :double) (dydt :double) + (absolute-scale :pointer) (dimension :size)) + :c-return (ptr :pointer) + :return (ptr) + :documentation + "Create a new control object which uses the same algorithm + as @code{gsl_odeiv_control_standard_new} but with an absolute error + which is scaled for each component by the array @var{scale_abs}. + The formula for @math{D_i} for this control object is + D_i = \epsilon_{abs} s_i + \epsilon_{rel} * (a_{y} |y_i| + a_{dydt} h |y'_i|) + where @math{s_i} is the @math{i}-th component of the array @var{scale_abs}. + The same error control heuristic is used by the Matlab @sc{ode} suite.") + +(defun-gsl control-alloc (control-type) + "gsl_odeiv_control_alloc" + ((control-type :pointer)) + :c-return (ptr :pointer) + :return (ptr) + :documentation + "Return a pointer to a newly allocated instance of a + control function of type @var{T}. This function is only needed for + defining new types of control functions. For most purposes the standard + control functions described above should be sufficient.") + +(defun-gsl initialize-control + (control absolute-error relative-error y-scaling dydt-scaling) + "gsl_odeiv_control_init" + ((control :pointer) (absolute-error :double) (relative-error :double) + (y-scaling :double) (dydt-scaling :double)) + :documentation + "Initialize the control function @var{control} with the + parameters absolute-error, relative-error, + y-scaling (scaling factor for y) and dydt-scaling (scaling + factor for derivatives).") + +(defun-gsl free-control (control) + "gsl_odeiv_control_free" + ((control :pointer)) + :c-return :void + :documentation + "Free all the memory associated with the control function + @var{control}.") + +(cffi:defcenum step-size-adjustment + (:step-size-decreased -1) :step-size-unchanged :step-size-increased) + +(defun-gsl adjust-stepsize (control stepper current-y y-error dydt step-size) + "gsl_odeiv_control_hadjust" + ((control :pointer) (stepper :pointer) (current-y :pointer) + (y-error :pointer) (dydt :pointer) (step-size :pointer)) + :c-return :enumerate + :enumeration step-size-adjustment + :documentation + "Adjust the step-size @var{h} using the control function + and the current values of current-y, y-error and @var{dydt}. + The stepping function stepper is also needed to determine the order + of the method. If the error in the y-values y-error is found to be + too large then the step-size is reduced and the function returns + :step-size-decreased. If the error is sufficiently small then + step-size may be increased and :step-size-increased is returned. The + function returns :step-size-unchanged if the step-size is + unchanged. The goal of the function is to estimate the largest + step-size which satisfies the user-specified accuracy requirements for + the current point.") + +(defun-gsl control-name (control) + "gsl_odeiv_control_name" + ((control :pointer)) + :c-return :string + :documentation "The name of the control function.") diff --git a/ordinary-differential-equations/ode-example.lisp b/ordinary-differential-equations/ode-example.lisp new file mode 100644 index 0000000000000000000000000000000000000000..7ff32c679cceb1b0808f9f446073e779a73d3f6f --- /dev/null +++ b/ordinary-differential-equations/ode-example.lisp @@ -0,0 +1,33 @@ +;******************************************************** +; file: ode-example.lisp +; description: Example ODE +; date: Sat Sep 29 2007 - 17:49 +; author: Liam Healy +; modified: Sat Sep 29 2007 - 18:51 +;******************************************************** +;;; $Id: $ + +(in-package :gsl) + +(def-ode-function vanderpol-function (time y dydt) + (declare (special mu)) + (setf (double-to-cl dydt 0) (double-to-cl y 1) + (double-to-cl dydt 1) + (- (- (double-to-cl y 0)) + (* mu (double-to-cl y 1) + (- (* (double-to-cl y 0) (double-to-cl y 0)) 1))))) + +(def-jacobian-function vanderpol-jacobian (time y dfdy dfdt) + (declare (special mu)) + (setf (double-to-cl dfdt 0) 0.0d0 + (double-to-cl dfdt 1) 0.0d0 + (double-to-cl dfdy 0) 0.0d0 + (double-to-cl dfdy 1) 1.0d0 + (double-to-cl dfdy 2) + (- (* -2 mu (double-to-cl y 0) (double-to-cl y 1)) 1.0d0) + (double-to-cl dfdy 3) + (* -1 mu (- (* (double-to-cl y 0) (double-to-cl y 0)) 1.0d0)))) + +(let ((stepper (step-allocate *step-rk8pd* 2) + + (step-free stepper)))) diff --git a/ordinary-differential-equations/ode-system.lisp b/ordinary-differential-equations/ode-system.lisp index 7d757593227460cdcb4d7364f0766d6a858ce952..82807e1a1682c2bf284f66774f8aa48d2b316536 100644 --- a/ordinary-differential-equations/ode-system.lisp +++ b/ordinary-differential-equations/ode-system.lisp @@ -3,7 +3,7 @@ ; description: ODE system setup ; date: Sun Apr 15 2007 - 14:19 ; author: Liam Healy -; modified: Mon Sep 24 2007 - 21:33 +; modified: Sat Sep 29 2007 - 17:54 ;******************************************************** ;;; $Id: $ @@ -19,7 +19,7 @@ (export 'def-ode-function 'def-jacobian-function) -(defmacro def-ode-function (name time dependent derivatives &body body) +(defmacro def-ode-function (name (time dependent derivatives) &body body) "Define a function that will evaluate the right-hand sides (derivatives) defining a set of ordinary differential equations (ODE). The function should take as input the time (a double-float) and @@ -46,7 +46,7 @@ (double-to-cl dydt 1) (double-to-cl y 0))) |# -(defmacro def-jacobian-function (name time dependent dfdy dfdt &body body) +(defmacro def-jacobian-function (name (time dependent dfdy dfdt) &body body) "Define a function that will evaluate the Jacobian (partial derivative) of the set of ordinary differential equations (ODE). The function should take as input the time (a double-float) and