seulex

Outdated egg!

This is an egg for CHICKEN 4, the unsupported old release. You're almost certainly looking for the CHICKEN 5 version of this egg, if it exists.

If it does not exist, there may be equivalent functionality provided by another egg; have a look at the egg index. Otherwise, please consider porting this egg to the current version of CHICKEN.

The Chicken seulex library provides bindings to the Fortran SEULEX solver for systems of stiff differential and differential algebraic equations of the form MY'=F(X,Y). SEULEX was written by E. Hairer AND G. Wanner and stands for Semi-implicit EULer EXtrapolation and computes its error estimates by means of polynomial extrapolation. Extrapolation solvers can achieve good performance for certain types of problems where high precision is required.

Installation notes

The Chicken seulex library must be compiled along with the SEULEX Fortran source code. Therefore, the user must have a Fortran compiler such as gfortran or g95 installed on their system. The following environment variables can be used to control the Fortran compilation process:

FORTRAN_COMPILER: path to Fortran compiler (default is "gfortran")
FORTRAN_FLAGS: flags to be passed to the Fortran compiler (default is "-fno-automatic -fPIC -g")
FORTRAN_LIBS: Fortran libraries to link to (default is "-lgfortran -lblas -llapack")

Library procedures

(seulex-create-solver XINIT YINIT FCN [JACOBIAN] [AUTONOMOUS] [USER-DATA] [DENSITY-COMPONENTS] [RELTOL] [ABSTOL]) => SEULEX-SOLVERprocedure

Creates and initializes an object representing a problem to be solved with the SEULEX solver.

Arguments XINIT and YINIT represent the initial conditions: XINIT is a scalar value for the independent variable, and YINIT is an SRFI-4 f64vector value containing the initial dependent variable values.

Argument FCN is used to compute the right-hand side function F and must be a procedure of the following form:

(LAMBDA T YY DATA)

(LAMBDA T YY)

depending on whether the USER-DATA optional argument is set, where

T: real-valued independent variable
YY: SRFI-4 f64vector with current variable values
DATA: is a user data object (if set)

This procedure must return a SRFI-4 f64vector containing the derivative values.

Optional keyword argument JACOBIAN is a procedure which will be used to compute the partial derivatives of F(X,Y) with respect to Y. If not given, it is computed internally by finite differences.

Optional keyword argument AUTONOMOUS is a boolean value that indicates whether F(X,Y) is independent of X (autonomous) or not (non-autonomous). The default is #t (autonomous).

Optional keyword argument USER-DATA is an object that will be passed as an additional argument to the right-hand side function.

Optional keyword argument DENSITY-COMPONENTS must be an SRFI-4 s32vector which indicates the components for which dense output is required.

Optional keyword arguments RELTOL and ABSTOL specify relative and absolute error tolerance, respectively. These both default to 1e-4.

seulex-destroy-solver SEULEX-SOLVERprocedure: Deallocates the memory associated with the given solver.

seulex-solve SEULEX-SOLVER Tprocedure: Integrates the system over an interval in the independent variable. This procedure returns either when the given T is reached, or when a root is found.

seulex-yy SEULEX-SOLVERprocedure: Returns the SRFI-4 f64vector value of current state values of the system.

seulex-nfcn SEULEX-SOLVERprocedure: Returns the number of function evaluations done so far.

seulex-njac SEULEX-SOLVERprocedure: Returns the number of Jacobian function evaluations done so far.

seulex-nstep SEULEX-SOLVERprocedure: Returns the number of computed steps.

seulex-ndec SEULEX-SOLVERprocedure: Returns the number of LU decompositions.

seulex-nsol SEULEX-SOLVERprocedure: Returns the number of backward-forward substitutions.

Example

;;
;; Hodgkin-Huxley model
;;

(use mathh seulex srfi-4)

(define neg -)
(define pow expt)

(define TEND  500.0)

  	                   
;; Model parameters

(define (I_stim t) 10)
(define C_m       1)
(define E_Na      50)
(define E_K       -77)
(define E_L       -54.4)
 (define gbar_Na   120)
(define gbar_K    36)
(define g_L       0.3)

;; Rate functions

(define (amf v)   (* 0.1    (/ (+ v 40)  (- 1.0 (exp (/ (neg (+ v 40)) 10))))))
(define (bmf v)   (* 4.0    (exp (/ (neg (+ v 65)) 18))))
(define (ahf v)   (* 0.07   (exp (/ (neg (+ v 65)) 20))))
(define (bhf v)   (/ 1.0    (+ 1.0 (exp (/ (neg (+ v 35)) 10)))))
(define (anf v)   (* 0.01   (/ (+ v 55) (- 1 (exp (/ (neg (+ v 55)) 10))))))
(define (bnf v)   (* 0.125  (exp (/ (neg (+ v 65)) 80))))

;; State functions

(define (minf v) (* 0.5 (+ 1 (tanh (/ (- v v1) v2)))))
(define (winf v) (* 0.5 (+ 1 (tanh (/ (- v v3) v4)))))
(define (lamw v) (* phi (cosh (/ (- v v3) (* 2 v4)))))
  	                   
;; Model equations

(define (rhs t yy)

  (let ((v (f64vector-ref yy 0))
	(m (f64vector-ref yy 1))
	(h (f64vector-ref yy 2))
	(n (f64vector-ref yy 3)))

    ;; transition rates at current step
    (let ((am  (amf v))
	  (an  (anf v))
	  (ah  (ahf v))
	  (bm  (bmf v))
	  (bn  (bnf v))
	  (bh  (bhf v))

	  (g_Na (* gbar_Na  (* h (pow m 3))))
	  (g_K  (* gbar_K   (pow n 4))))
      
      (let (

	    ;; currents
	    (I_Na   (* (- v E_Na) g_Na))
	    (I_K    (* (- v E_K)  g_K))
	    (I_L    (* g_L  (- v E_L))))
		  
	(let (
	      ;; state equations
	      (dm (- (* am (- 1 m))  (* bm m)))
	      (dh (- (* ah (- 1 h))  (* bh h)))
	      (dn (- (* an (- 1 n))  (* bn n)))
	      (dv (/ (- (I_stim t) I_L I_Na I_K) C_m))
	      )
   
	  (f64vector dv dm dh dn)
	  
	  )))
    ))
  
 (let ((yy (f64vector -65  0.052 0.596 0.317)) ;; v m h n

	;; Integration limits 
	(t0  0.0)
	(tf  TEND)
	(dt  1e-2))
   
    ;; solver initialization 
    (let ((solver (seulex-create-solver
		   t0 yy rhs  
		   abstol: 1e-4
		   reltol: 1e-4)))

      ;; In loop, call solver, print results, and test for error. 
      
      (let recur ((tnext (+ t0 dt)) (iout 1))

	(let ((flag  (seulex-solve solver tnext)))
	  (if (negative? flag) (error 'main "SEULEX solver error" flag))
 
         (print-results solver tnext)

	  (if (< tnext tf)
	      (recur (+ tnext dt) (+ 1 iout)))
	  ))
      

     (seulex-destroy-solver solver)
      
(define (print-results solver t)
  (let ((yy (seulex-yy solver)))
    (printf "~A ~A ~A ~A ~A~%" 
	    t 
	     (f64vector-ref yy 0)
     (f64vector-ref yy 1)
     (f64vector-ref yy 2)
     (f64vector-ref yy 3)
     )))

Version History

1.0 Initial release

License

 The SEULEX solver was written by E. HAIRER AND G. WANNER, 
 UNIVERSITE DE GENEVE, DEPT. DE MATHEMATIQUES

 Chicken Scheme bindings for SEULEX are copyright 2011-2012 Ivan Raikov.

This program is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or (at
your option) any later version.

This program is distributed in the hope that it will be useful, but
WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
General Public License for more details.

A full copy of the GPL license can be found at
<http://www.gnu.org/licenses/>.

chickadee » seulex

Identifier search