Commit b6d90200 authored by Liam Healy's avatar Liam Healy

New macro with-ode-integration, new ODE tests

The macro with-ode-integration has been significantly altered, using
marrays and symbol macros to make integration forms easier.  There are
now 10 ode tests, one for each of the supplied steppers.
parent 4e83fcc4
......@@ -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>
......
;; 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))
;; 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.
......
;; 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))
;; 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))))
......@@ -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))))
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment