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

Last change on this file since 25509 was 25509, checked in by Ivan Raikov, 9 years ago

sundials: added M-L example with IDA solver

File size: 4.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(define (ida-main)
69
70  (define (ressc t yy yp)
71    (let ((dd (rhs t yy)))
72
73      (let ((v (- (f64vector-ref dd 0) (f64vector-ref yp 0)))
74            (w (- (f64vector-ref dd 1) (f64vector-ref yp 1))))
75        (f64vector v w))))
76 
77  (let* ((yy (f64vector -60.899 0.0149)) ;; v w
78         (yp (rhs 0.0 yy))
79         
80         ;; Integration limits
81         (t0  0.0)
82         (tf  TEND)
83         (dt  1e-2))
84   
85
86    ;; IDA initialization
87    (let ((solver (ida-create-solver t0 tf yy yp ressc 
88                                     abstol: 1e-14
89                                     reltol: 1e-14)))
90
91      ;; In loop, call IDASolve, print results, and test for error.
92     
93      (let recur ((tnext (+ t0 dt)) (iout 1))
94
95        (let ((flag  (ida-solve solver tnext)))
96          (if (negative? flag) (error 'main "IDA solver error" flag))
97
98         
99          (let ((yy (ida-yy solver)))
100            (printf "~A ~A ~A~%" 
101                    tnext
102                    (f64vector-ref yy 0)
103                    (f64vector-ref yy 1)
104                    ))
105
106
107          (if (< tnext tf)
108              (recur (+ tnext dt) (+ 1 iout)))
109          ))
110     
111      (ida-destroy-solver solver)
112     
113      )))
114
115
116
117(define (cvode-main)
118 
119  (let ((yy (f64vector -60.899 0.0149));; v w
120
121        ;; Integration limits
122        (t0  0.0)
123        (tf  TEND)
124        (dt  1e-2))
125   
126    ;; CVODE initialization
127    (let ((solver (cvode-create-solver
128                   t0 tf yy rhs
129                   abstol: 1e-4
130                   reltol: 1e-4)))
131
132      ;; In loop, call CVodeSolve, print results, and test for error.
133
134      (let recur ((tnext (+ t0 dt)) (iout 1))
135
136        (let ((flag  (cvode-solve solver tnext)))
137          (if (negative? flag) (error 'main "CVODE solver error" flag))
138
139
140          (if (< tnext tf)
141              (recur (+ tnext dt) (+ 1 iout)))
142          ))
143     
144      (let ((yy (cvode-yy solver)))
145        (let ((v (f64vector-ref yy 0))
146              (w (f64vector-ref yy 1)))
147          (assert (< (abs (- v -30.0536210523504)) 1e-12))
148          (assert (< (abs (- w 0.0402640532124851 )) 1e-12))))
149
150      (cvode-destroy-solver solver)
151     
152      )))
153
154(define (cvode-main/unsafe)
155 
156  (let ((yy (f64vector -60.899 0.0149)) ;; v w
157
158        ;; Integration limits
159        (t0  0.0)
160        (tf  TEND)
161        (dt  1e-2))
162   
163    ;; CVODE initialization
164    (let ((solver (cvode-create-solver/unsafe
165                   t0 tf yy rhs/unsafe
166                   abstol: 1e-4
167                   reltol: 1e-4)))
168
169      ;; In loop, call CVodeSolve, print results, and test for error.
170
171      (let recur ((tnext (+ t0 dt)) (iout 1))
172
173        (let ((flag  (cvode-solve solver tnext)))
174          (if (negative? flag) (error 'main "CVODE solver error" flag))
175
176          (if (< tnext tf)
177              (recur (+ tnext dt) (+ 1 iout)))
178          ))
179     
180      (let ((yy (cvode-yy solver)))
181        (let ((v (f64vector-ref yy 0))
182              (w (f64vector-ref yy 1)))
183          (assert (< (abs (- v -30.0536210523504)) 1e-12))
184          (assert (< (abs (- w 0.0402640532124851 )) 1e-12))))
185
186      (cvode-destroy-solver solver)
187     
188      )))
189
190
191(define (print-results solver t)
192  (let ((yy (cvode-yy solver)))
193    (printf "~A ~A ~A ~A ~A ~A~%" 
194            t 
195            (f64vector-ref yy 0)
196            (f64vector-ref yy 1)
197            (cvode-get-last-order solver)
198            (cvode-get-num-steps solver)
199            (cvode-get-last-step solver)
200            )))
201     
202;(ida-main)     
203(cvode-main)
204(cvode-main/unsafe)
Note: See TracBrowser for help on using the repository browser.