source: project/release/4/sundials/trunk/tests/ml.scm @ 25082

Last change on this file since 25082 was 25082, checked in by Ivan Raikov, 10 years ago

sundials: introduced safe and unsafe variants of the solver creation procedures

File size: 3.2 KB
Line 
1;;
2;; Morris-Lecar model
3;;
4
5(use mathh sundials srfi-4)
6
7
8(define TEND  500.0)
9
10                           
11;; Model parameters
12
13(define Istim  20)
14(define vk     -84)
15(define vl     -60)
16(define vca    120)
17(define gk     2.5)
18(define gl     0.5)
19(define gca    1)
20(define c      2)
21(define v1     -1.2)
22(define v2     18)
23(define v3     12)
24(define v4     17)
25(define phi    0.0667)
26
27
28;; State functions
29
30(define (minf v) (* 0.5 (+ 1 (tanh (/ (- v v1) v2)))))
31(define (winf v) (* 0.5 (+ 1 (tanh (/ (- v v3) v4)))))
32(define (lamw v) (* phi (cosh (/ (- v v3) (* 2 v4)))))
33                           
34;; Model equations
35
36(define (rhs t yy)
37  (let ((v (f64vector-ref yy 0))
38        (w (f64vector-ref yy 1)))
39
40  (let ((ica (* gca (* (minf v)  (- vca v))))
41        (ik  (* gk  (* w (- vk v )))))
42   
43    (let ((dv (/ (+ Istim (* gl (- vl v)) ica ik) c))
44          (dw (* (lamw v) (- (winf v) w))))
45   
46      (f64vector dv dw)
47
48      ))
49  ))
50
51(define (rhs/unsafe t yy yp _)
52  (let ((v (pointer-f64-ref yy))
53        (w (pointer-f64-ref (pointer+-f64 yy 1))))
54
55  (let ((ica (* gca (* (minf v)  (- vca v))))
56        (ik  (* gk  (* w (- vk v )))))
57   
58    (let ((dv (/ (+ Istim (* gl (- vl v)) ica ik) c))
59          (dw (* (lamw v) (- (winf v) w))))
60   
61      (pointer-f64-set! yp dv)
62      (pointer-f64-set! (pointer+-f64 yp 1) dw)
63
64      ))
65  ))
66
67
68
69(define (main)
70 
71  (let ((yy (f64vector -60.899 0.0149)) ;; v w
72
73        ;; Integration limits
74        (t0  0.0)
75        (tf  TEND)
76        (dt  1e-2))
77   
78    ;; CVODE initialization
79    (let ((solver (cvode-create-solver
80                   t0 tf yy rhs
81                   abstol: 1e-4
82                   reltol: 1e-4)))
83
84      ;; In loop, call CVodeSolve, print results, and test for error.
85
86      (let recur ((tnext (+ t0 dt)) (iout 1))
87
88        (let ((flag  (cvode-solve solver tnext)))
89          (if (negative? flag) (error 'main "CVODE solver error" flag))
90
91          (if (< tnext tf)
92              (recur (+ tnext dt) (+ 1 iout)))
93          ))
94     
95      (let ((yy (cvode-yy solver)))
96        (let ((v (f64vector-ref yy 0))
97              (w (f64vector-ref yy 1)))
98          (assert (< (abs (- v -30.0536210523504)) 1e-12))
99          (assert (< (abs (- w 0.0402640532124851 )) 1e-12))))
100
101      (cvode-destroy-solver solver)
102     
103      )))
104
105(define (main/unsafe)
106 
107  (let ((yy (f64vector -60.899 0.0149)) ;; v w
108
109        ;; Integration limits
110        (t0  0.0)
111        (tf  TEND)
112        (dt  1e-2))
113   
114    ;; CVODE initialization
115    (let ((solver (cvode-create-solver/unsafe
116                   t0 tf yy rhs/unsafe
117                   abstol: 1e-4
118                   reltol: 1e-4)))
119
120      ;; In loop, call CVodeSolve, print results, and test for error.
121
122      (let recur ((tnext (+ t0 dt)) (iout 1))
123
124        (let ((flag  (cvode-solve solver tnext)))
125          (if (negative? flag) (error 'main "CVODE solver error" flag))
126
127          (if (< tnext tf)
128              (recur (+ tnext dt) (+ 1 iout)))
129          ))
130     
131      (let ((yy (cvode-yy solver)))
132        (let ((v (f64vector-ref yy 0))
133              (w (f64vector-ref yy 1)))
134          (assert (< (abs (- v -30.0536210523504)) 1e-12))
135          (assert (< (abs (- w 0.0402640532124851 )) 1e-12))))
136
137      (cvode-destroy-solver solver)
138     
139      )))
140
141
142(define (print-results solver t)
143  (let ((yy (cvode-yy solver)))
144    (printf "~A ~A ~A ~A ~A ~A~%" 
145            t 
146            (f64vector-ref yy 0)
147            (f64vector-ref yy 1)
148            (cvode-get-last-order solver)
149            (cvode-get-num-steps solver)
150            (cvode-get-last-step solver)
151            )))
152     
153     
154(main)
155(main/unsafe)
Note: See TracBrowser for help on using the repository browser.