Skip to content
Snippets Groups Projects
Commit 159d983c authored by Liam Healy's avatar Liam Healy
Browse files

ODE stepper, evolution as mobjects

Defined mobjects 'ode-step and 'ode-evolution; "reset" is the
reinitialize-instance method.  Test 'ode works.
parent f755beaf
No related branches found
No related tags found
No related merge requests found
...@@ -311,6 +311,7 @@ random-number-generator quasi-random-number-generator discrete-random ...@@ -311,6 +311,7 @@ random-number-generator quasi-random-number-generator discrete-random
polynomial-complex-workspace integration-workspace polynomial-complex-workspace integration-workspace
eigen-symm eigen-symmv eigen-herm eigen-hermv eigen-symm eigen-symmv eigen-herm eigen-hermv
monte-carlo-plain monte-carlo-miser monte-carlo-vegas monte-carlo-plain monte-carlo-miser monte-carlo-vegas
ode-stepper ode-evolution
</pre> </pre>
<p> <p>
An instance may be created with a function whose name is "make-" An instance may be created with a function whose name is "make-"
...@@ -357,7 +358,7 @@ and arrays used internally or for function return. ...@@ -357,7 +358,7 @@ and arrays used internally or for function return.
<!-- Created: Feb 25 2005 --> <!-- Created: Feb 25 2005 -->
<!-- hhmts start --> <!-- hhmts start -->
<small> <small>
Time-stamp: <2009-01-16 09:22:35EST index.html> Time-stamp: <2009-01-20 22:50:29EST index.html>
</small> </small>
<!-- hhmts end --> <!-- hhmts end -->
</div> </div>
......
;; Adaptive step-size control ;; Adaptive step-size control
;; Liam Healy 2008-02-17 17:30:04EST control.lisp ;; Liam Healy 2008-02-17 17:30:04EST control.lisp
;; Time-stamp: <2008-08-21 22:08:55EDT control.lisp> ;; Time-stamp: <2009-01-20 19:39:08EST control.lisp>
;; $Id$ ;; $Id$
(in-package :gsl) (in-package :gsl)
(defmfun new-standard-control (absolute-error relative-error y dydt) (defmfun new-standard-control (absolute-error relative-error y-scaling dydt-scaling)
"gsl_odeiv_control_standard_new" "gsl_odeiv_control_standard_new"
((absolute-error :double) (relative-error :double) ((absolute-error :double) (relative-error :double)
(ay :double) (adydt :double)) (y-scaling :double) (dydt-scaling :double))
:c-return (ptr :pointer) :c-return (ptr :pointer)
:return (ptr) :return (ptr)
:documentation ; FDL :documentation ; FDL
"The standard control object is a four parameter heuristic based on "The standard control object is a four parameter heuristic based on
absolute and relative errors absolute-error and relative-error, and absolute and relative errors absolute-error and relative-error, and
scaling factors ay and adydt for the system state y(t) and derivatives scaling factors y-scaling and dydt-scaling for the system state y(t) and derivatives
y'(t) respectively. y'(t) respectively.
The step-size adjustment procedure for this method begins by computing The step-size adjustment procedure for this method begins by computing
...@@ -47,8 +47,8 @@ ...@@ -47,8 +47,8 @@
"Create a new control object which will keep the local "Create a new control object which will keep the local
error on each step within an absolute error of absolute-error and error on each step within an absolute error of absolute-error and
relative error of relative-error with respect to the solution y_i(t). relative error of relative-error with respect to the solution y_i(t).
This is equivalent to the standard control object with ay=1 and This is equivalent to the standard control object with y-scaling=1 and
adydt=0.") dydt-scaling=0.")
(defmfun new-yp-control (absolute-error relative-error) (defmfun new-yp-control (absolute-error relative-error)
"gsl_odeiv_control_yp_new" "gsl_odeiv_control_yp_new"
...@@ -60,12 +60,13 @@ ...@@ -60,12 +60,13 @@
error on each step within an absolute error of absolute-error and error on each step within an absolute error of absolute-error and
relative error of relative-error with respect to the derivatives of the relative error of relative-error with respect to the derivatives of the
solution y'_i(t). This is equivalent to the standard control solution y'_i(t). This is equivalent to the standard control
object with ay=0 and adydt=1.") object with y-scaling=0 and dydt-scaling=1.")
(defmfun new-scaled-control (defmfun new-scaled-control
(absolute-error relative-error y dydt absolute-scale dimension) (absolute-error relative-error y-scaling dydt-scaling absolute-scale dimension)
"gsl_odeiv_control_scaled_new" "gsl_odeiv_control_scaled_new"
((absolute-error :double) (relative-error :double) (y :double) (dydt :double) ((absolute-error :double) (relative-error :double)
(y-scaling :double) (dydt-scaling :double)
(absolute-scale :pointer) (dimension sizet)) (absolute-scale :pointer) (dimension sizet))
:c-return (ptr :pointer) :c-return (ptr :pointer)
:return (ptr) :return (ptr)
......
;; Evolution functions for ODE integration. ;; Evolution functions for ODE integration.
;; Liam Healy, Sun Sep 30 2007 - 14:31 ;; Liam Healy, Sun Sep 30 2007 - 14:31
;; Time-stamp: <2008-08-21 22:07:27EDT evolution.lisp> ;; Time-stamp: <2009-01-20 22:42:27EST evolution.lisp>
;; $Id$ ;; $Id$
(in-package :gsl) (in-package :gsl)
(defmfun allocate-evolution (dimension) (defmobject ode-evolution "gsl_odeiv_evolve"
"gsl_odeiv_evolve_alloc" ((dimensions sizet))
((dimension sizet)) "evolution for ordinary differential equations"
:c-return :pointer "Make an object to advance the ODE solution."
:documentation ; FDL "reset" nil)
"Allocate a new instance of an evolution function
for a system of dimension dimensions and return the pointer.")
(defmfun apply-evolution (defmfun apply-evolution
(evolve control step dydt time max-time step-size y) (evolve control step dydt time max-time step-size y)
"gsl_odeiv_evolve_apply" "gsl_odeiv_evolve_apply"
((evolve :pointer) (control :pointer) (step :pointer) (((mpointer evolve) :pointer) (control :pointer) ((mpointer step) :pointer)
(dydt :pointer) (time :pointer) (max-time :double) (dydt :pointer) (time :pointer) (max-time :double)
(step-size :pointer) (y :pointer)) (step-size :pointer) (y :pointer))
:documentation ; FDL :documentation ; FDL
...@@ -31,17 +29,3 @@ ...@@ -31,17 +29,3 @@
time max-time is guaranteed not to be exceeded by the time-step. On the time max-time is guaranteed not to be exceeded by the time-step. On the
final time-step the value of time will be set to t1 exactly.") final time-step the value of time will be set to t1 exactly.")
(defmfun reset-evolution (evolve)
"gsl_odeiv_evolve_reset"
((evolve :pointer))
:documentation
"Reset the evolution function evolve. It should be used
whenever the next use of evolve will not be a continuation of a
previous step.")
(defmfun free-evolution (evolve)
"gsl_odeiv_evolve_free"
((evolve :pointer))
:c-return :void
:documentation
"Frees all the memory associated with the evolution function.")
;; ODE system setup ;; ODE system setup
;; Liam Healy, Sun Apr 15 2007 - 14:19 ;; Liam Healy, Sun Apr 15 2007 - 14:19
;; Time-stamp: <2009-01-18 18:23:07EST ode-system.lisp> ;; Time-stamp: <2009-01-20 22:47:16EST ode-system.lisp>
;; $Id$ ;; $Id$
(in-package :gsl) (in-package :gsl)
...@@ -51,23 +51,20 @@ ...@@ -51,23 +51,20 @@
with-c-double." with-c-double."
(let ((ctime (make-symbol "CTIME")) (let ((ctime (make-symbol "CTIME"))
(cstep (make-symbol "CSTEP"))) (cstep (make-symbol "CSTEP")))
`(let ((stepper (step-allocate ,stepper ,dimensions)) `(let ((stepper (make-ode-stepper ,stepper ,dimensions))
(control (new-y-control ,absolute-error ,relative-error)) (control (new-y-control ,absolute-error ,relative-error))
(evolve (allocate-evolution ,dimensions))) (evolve (make-ode-evolution ,dimensions)))
(unwind-protect (unwind-protect
(cffi:with-foreign-objects (cffi:with-foreign-objects
((,(if (listp dependent) (first dependent) dependent) ((,(if (listp dependent) (first dependent) dependent)
:double ,dimensions) (,ctime :double) (,cstep :double)) :double ,dimensions) (,ctime :double) (,cstep :double))
(setf (setf
(dcref ,cstep) ,step-size (dcref ,cstep) ,step-size
(dcref ,ctime) ,time) (dcref ,ctime) ,time)
,(if (listp dependent) ,(if (listp dependent)
`(with-c-double ,dependent `(with-c-double ,dependent
(symbol-macrolet ((,time ,ctime) (,step-size ,cstep)) (symbol-macrolet ((,time ,ctime) (,step-size ,cstep))
,@body)) ,@body))
`(symbol-macrolet ((,time ,ctime) (,step-size ,cstep)) `(symbol-macrolet ((,time ,ctime) (,step-size ,cstep))
,@body))) ,@body)))
(free-evolution evolve) (free-control control)))))
(free-control control)
(step-free stepper)))))
;; Stepping functions for ODE systems. ;; Stepping functions for ODE systems.
;; Liam Healy, Mon Sep 24 2007 - 21:33 ;; Liam Healy, Mon Sep 24 2007 - 21:33
;; Time-stamp: <2008-08-21 22:05:20EDT stepping.lisp> ;; Time-stamp: <2009-01-20 19:09:55EST stepping.lisp>
;; $Id$ ;; $Id$
(in-package :gsl) (in-package :gsl)
(defmfun step-allocate (step-type dim) (defmobject ode-stepper "gsl_odeiv_step"
"gsl_odeiv_step_alloc" ((type :pointer) (dimensions sizet))
((step-type :pointer) (dim sizet)) "stepper for ordinary differential equations"
:c-return :pointer "Make a stepper for ordinary differential equations. The type is
:documentation ; FDL one of the GSL-supplied types defined in stepping.lisp, and
"Allocate a new instance of a stepping function of dimensions is the number of dependent variables. This instance
type step-type for a system of dim dimensions, should be reinitialized whenever the next use of it will not be a
returning the pointer.") continuation of a previous step."
"reset" nil)
(defmfun step-reset (stepper)
"gsl_odeiv_step_reset" (defmfun name ((object ode-stepper))
((stepper :pointer))
:documentation ; FDL
"Reset the stepping function. It should be used whenever
the next use of it will not be a continuation of a previous
step.")
(defmfun step-free (stepper)
"gsl_odeiv_step_free"
((stepper :pointer))
:c-return :void
:documentation ; FDL
"Free all the memory associated with the stepping function.")
(defmfun step-name (stepper)
"gsl_odeiv_step_name" "gsl_odeiv_step_name"
((stepper :pointer)) (((mpointer object) :pointer))
:definition :method
:c-return :string :c-return :string
:documentation ; FDL :documentation ; FDL
"The name of the stepping function.") "The name of the stepping function.")
(defmfun step-order (stepper) (defmfun step-order (stepper)
"gsl_odeiv_step_order" "gsl_odeiv_step_order"
((stepper :pointer)) (((mpointer stepper) :pointer))
:c-return :uint :c-return :uint
:documentation ; FDL :documentation ; FDL
"The order of the stepping function on the previous "The order of the stepping function on the previous
...@@ -47,7 +34,7 @@ ...@@ -47,7 +34,7 @@
(defmfun step-apply (defmfun step-apply
(stepper time step-size y yerr dydt-in dydt-out dydt) (stepper time step-size y yerr dydt-in dydt-out dydt)
"gsl_odeiv_step_apply" "gsl_odeiv_step_apply"
((stepper :pointer) (((mpointer stepper) :pointer)
(time :double) (time :double)
(step-size :double) (step-size :double)
(y :pointer) (y :pointer)
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment