((section 2 "Outdated egg!" (p "This is an egg for CHICKEN 4, the unsupported old release.  You're almost certainly looking for " (int-link "/eggref/5/seulex" "the CHICKEN 5 version of this egg") ", if it exists.") (p "If it does not exist, there may be equivalent functionality provided by another egg; have a look at the " (link "https://wiki.call-cc.org/chicken-projects/egg-index-5.html" "egg index") ". Otherwise, please consider porting this egg to the current version of CHICKEN.") (tags "eggs") (toc)) (section 2 "seulex" (section 3 "Description" (p "The Chicken " (tt "seulex") " library provides bindings to the Fortran " (tt "SEULEX") " solver for systems of stiff differential and differential algebraic equations of the form " (tt "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.")) (section 3 "Installation notes" (p "The Chicken " (tt "seulex") " library must be compiled along with the " (tt "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:") (dl (dt (tt "FORTRAN_COMPILER")) (dd "path to Fortran compiler (default is \"gfortran\")") (dt (tt "FORTRAN_FLAGS")) (dd "flags to be passed to the Fortran compiler (default is \"-fno-automatic -fPIC -g\")") (dt (tt "FORTRAN_LIBS")) (dd "Fortran libraries to link to  (default is \"-lgfortran -lblas -llapack\")"))) (section 3 "Library procedures" (def (sig (procedure "(seulex-create-solver XINIT YINIT FCN [JACOBIAN] [AUTONOMOUS] [USER-DATA] [DENSITY-COMPONENTS] [RELTOL] [ABSTOL]) => SEULEX-SOLVER" (id seulex-create-solver))) (p "Creates and initializes an object representing a problem to be solved with the SEULEX solver.") (p "Arguments " (tt "XINIT") " and " (tt "YINIT") " represent the initial conditions: " (tt "XINIT") " is a scalar value for the independent variable, and " (tt "YINIT") " is an SRFI-4 " (tt "f64vector") " value containing the initial dependent variable values.") (p "Argument " (tt "FCN") " is used to compute the right-hand side function " (tt "F") " and must be a procedure of the following form:") (pre "(LAMBDA T YY DATA)") (p "or") (pre "(LAMBDA T YY)") (p "depending on whether the " (tt "USER-DATA") " optional argument is set, where") (dl (dt (tt "T")) (dd "real-valued independent variable") (dt (tt "YY")) (dd "SRFI-4 " (tt "f64vector") " with current variable values") (dt (tt "DATA")) (dd "is a user data object (if set)")) (p "This procedure must return a SRFI-4 " (tt "f64vector") " containing the derivative values.") (p "Optional keyword argument " (tt "JACOBIAN") " is a procedure which will be used to compute the partial derivatives of " (tt "F(X,Y)") " with respect to " (tt "Y") ". If not given, it is computed internally by finite differences.") (p "Optional keyword argument " (tt "AUTONOMOUS") " is a boolean value that indicates whether " (tt "F(X,Y)") " is independent of " (tt "X") " (autonomous) or not (non-autonomous). The default is " (tt "#t") " (autonomous).") (p "Optional keyword argument " (tt "USER-DATA") " is an object that will be passed as an additional argument to the right-hand side function.") (p "Optional keyword argument " (tt "DENSITY-COMPONENTS") " must be an SRFI-4 " (tt "s32vector") " which indicates the components for which dense output is required.") (p "Optional keyword arguments " (tt "RELTOL") " and " (tt "ABSTOL") " specify relative and absolute error tolerance, respectively. These both default to 1e-4.")) (def (sig (procedure "(seulex-destroy-solver SEULEX-SOLVER)" (id seulex-destroy-solver))) (p "Deallocates the memory associated with the given solver.")) (def (sig (procedure "(seulex-solve SEULEX-SOLVER T)" (id seulex-solve))) (p "Integrates the system over an interval in the independent variable. This procedure returns either when the given " (tt "T") " is reached, or when a root is found.")) (def (sig (procedure "(seulex-yy SEULEX-SOLVER)" (id seulex-yy))) (p "Returns the SRFI-4 " (tt "f64vector") " value of current state values of the system.")) (def (sig (procedure "(seulex-nfcn SEULEX-SOLVER)" (id seulex-nfcn))) (p "Returns the number of function evaluations done so far.")) (def (sig (procedure "(seulex-njac SEULEX-SOLVER)" (id seulex-njac))) (p "Returns the number of Jacobian function evaluations done so far.")) (def (sig (procedure "(seulex-nstep SEULEX-SOLVER)" (id seulex-nstep))) (p "Returns the number of computed steps.")) (def (sig (procedure "(seulex-ndec SEULEX-SOLVER)" (id seulex-ndec))) (p "Returns the number of LU decompositions.")) (def (sig (procedure "(seulex-nsol SEULEX-SOLVER)" (id seulex-nsol))) (p "Returns the number of backward-forward substitutions."))) (section 3 "Example" (pre ";;\n;; Hodgkin-Huxley model\n;;\n\n(use mathh seulex srfi-4)\n\n(define neg -)\n(define pow expt)\n\n(define TEND  500.0)\n\n  \t                   \n;; Model parameters\n\n(define (I_stim t) 10)\n(define C_m       1)\n(define E_Na      50)\n(define E_K       -77)\n(define E_L       -54.4)\n (define gbar_Na   120)\n(define gbar_K    36)\n(define g_L       0.3)\n\n;; Rate functions\n\n(define (amf v)   (* 0.1    (/ (+ v 40)  (- 1.0 (exp (/ (neg (+ v 40)) 10))))))\n(define (bmf v)   (* 4.0    (exp (/ (neg (+ v 65)) 18))))\n(define (ahf v)   (* 0.07   (exp (/ (neg (+ v 65)) 20))))\n(define (bhf v)   (/ 1.0    (+ 1.0 (exp (/ (neg (+ v 35)) 10)))))\n(define (anf v)   (* 0.01   (/ (+ v 55) (- 1 (exp (/ (neg (+ v 55)) 10))))))\n(define (bnf v)   (* 0.125  (exp (/ (neg (+ v 65)) 80))))\n\n;; State functions\n\n(define (minf v) (* 0.5 (+ 1 (tanh (/ (- v v1) v2)))))\n(define (winf v) (* 0.5 (+ 1 (tanh (/ (- v v3) v4)))))\n(define (lamw v) (* phi (cosh (/ (- v v3) (* 2 v4)))))\n  \t                   \n;; Model equations\n\n(define (rhs t yy)\n\n  (let ((v (f64vector-ref yy 0))\n\t(m (f64vector-ref yy 1))\n\t(h (f64vector-ref yy 2))\n\t(n (f64vector-ref yy 3)))\n\n    ;; transition rates at current step\n    (let ((am  (amf v))\n\t  (an  (anf v))\n\t  (ah  (ahf v))\n\t  (bm  (bmf v))\n\t  (bn  (bnf v))\n\t  (bh  (bhf v))\n\n\t  (g_Na (* gbar_Na  (* h (pow m 3))))\n\t  (g_K  (* gbar_K   (pow n 4))))\n      \n      (let (\n\n\t    ;; currents\n\t    (I_Na   (* (- v E_Na) g_Na))\n\t    (I_K    (* (- v E_K)  g_K))\n\t    (I_L    (* g_L  (- v E_L))))\n\t\t  \n\t(let (\n\t      ;; state equations\n\t      (dm (- (* am (- 1 m))  (* bm m)))\n\t      (dh (- (* ah (- 1 h))  (* bh h)))\n\t      (dn (- (* an (- 1 n))  (* bn n)))\n\t      (dv (/ (- (I_stim t) I_L I_Na I_K) C_m))\n\t      )\n   \n\t  (f64vector dv dm dh dn)\n\t  \n\t  )))\n    ))\n  \n (let ((yy (f64vector -65  0.052 0.596 0.317)) ;; v m h n\n\n\t;; Integration limits \n\t(t0  0.0)\n\t(tf  TEND)\n\t(dt  1e-2))\n   \n    ;; solver initialization \n    (let ((solver (seulex-create-solver\n\t\t   t0 yy rhs  \n\t\t   abstol: 1e-4\n\t\t   reltol: 1e-4)))\n\n      ;; In loop, call solver, print results, and test for error. \n      \n      (let recur ((tnext (+ t0 dt)) (iout 1))\n\n\t(let ((flag  (seulex-solve solver tnext)))\n\t  (if (negative? flag) (error 'main \"SEULEX solver error\" flag))\n \n         (print-results solver tnext)\n\n\t  (if (< tnext tf)\n\t      (recur (+ tnext dt) (+ 1 iout)))\n\t  ))\n      \n\n     (seulex-destroy-solver solver)\n      \n(define (print-results solver t)\n  (let ((yy (seulex-yy solver)))\n    (printf \"~A ~A ~A ~A ~A~%\" \n\t    t \n\t     (f64vector-ref yy 0)\n     (f64vector-ref yy 1)\n     (f64vector-ref yy 2)\n     (f64vector-ref yy 3)\n     )))")) (section 3 "Version History" (ul (li "1.0 Initial release"))) (section 3 "License" (pre " The SEULEX solver was written by E. HAIRER AND G. WANNER, \n UNIVERSITE DE GENEVE, DEPT. DE MATHEMATIQUES\n\n Chicken Scheme bindings for SEULEX are copyright 2011-2012 Ivan Raikov.\n\nThis program is free software: you can redistribute it and/or modify\nit under the terms of the GNU General Public License as published by\nthe Free Software Foundation, either version 3 of the License, or (at\nyour option) any later version.\n\nThis program is distributed in the hope that it will be useful, but\nWITHOUT ANY WARRANTY; without even the implied warranty of\nMERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU\nGeneral Public License for more details.\n\nA full copy of the GPL license can be found at\n<http://www.gnu.org/licenses/>.") (pre " \n "))))