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

Last change on this file since 30965 was 30965, checked in by Ivan Raikov, 6 years ago

sundials: added regular spiking izhikevich model

File size: 4.1 KB
Line 
1;;
2;; Morris-Lecar model
3;;
4
5(use mathh sundials srfi-4)
6
7
8(define TEND  1000.0)
9
10                           
11;; Model parameters
12
13(define Istim  50.0)
14(define vl     -50)
15(define vk     -70)
16(define vca    100)
17(define gl     2.0)
18(define gk     8.0)
19(define gca    4.0)
20(define c      20.0)
21(define v1     -1.0)
22(define v2     15)
23(define v3     10)
24(define v4     14.5)
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 yy yp ressc 
88                                     tstop: tf
89                                     abstol: 1e-14
90                                     reltol: 1e-14)))
91
92      ;; In loop, call IDASolve, print results, and test for error.
93     
94      (let recur ((tnext (+ t0 dt)) (iout 1))
95
96        (let ((flag  (ida-solve solver tnext)))
97          (if (negative? flag) (error 'main "IDA solver error" flag))
98
99;;        (print-results solver tnext)
100
101          (if (< tnext tf)
102              (recur (+ tnext dt) (+ 1 iout)))
103          ))
104     
105      (ida-destroy-solver solver)
106     
107      )))
108
109
110
111(define (cvode-main)
112 
113  (let ((yy (f64vector -60.899 0.0149));; v w
114
115        ;; Integration limits
116        (t0  0.0)
117        (tf  TEND)
118        (dt  0.1))
119   
120    ;; CVODE initialization
121    (let ((solver (cvode-create-solver
122                   t0 yy rhs
123                   tstop: tf
124                   abstol: 1e-6
125                   reltol: 1e-6)))
126
127      ;; In loop, call CVodeSolve, print results, and test for error.
128
129      (let recur ((tnext (+ t0 dt)) (iout 1))
130
131        (let ((flag  (cvode-solve solver tnext)))
132          (if (negative? flag) (error 'main "CVODE solver error" flag))
133
134          (print-results solver tnext)
135
136          (if (< tnext tf)
137              (recur (+ tnext dt) (+ 1 iout)))
138          ))
139     
140#;      (let ((yy (cvode-yy solver)))
141        (let ((v (f64vector-ref yy 0))
142              (w (f64vector-ref yy 1)))
143          (printf "v = ~A w = ~A~%" v w)
144          ))
145
146      (cvode-destroy-solver solver)
147     
148      )))
149
150(define (cvode-main/unsafe)
151 
152  (let ((yy (f64vector -60.899 0.0149)) ;; v w
153
154        ;; Integration limits
155        (t0  0.0)
156        (tf  TEND)
157        (dt  0.1))
158   
159    ;; CVODE initialization
160    (let ((solver (cvode-create-solver/unsafe
161                   t0 yy rhs/unsafe
162                   tstop: tf
163                   abstol: 1e-4
164                   reltol: 1e-4)))
165
166      ;; In loop, call CVodeSolve, print results, and test for error.
167
168      (let recur ((tnext (+ t0 dt)) (iout 1))
169
170        (let ((flag  (cvode-solve solver tnext)))
171          (if (negative? flag) (error 'main "CVODE solver error" flag))
172
173;         (print-results solver tnext)
174
175          (if (< tnext tf)
176              (recur (+ tnext dt) (+ 1 iout)))
177          ))
178     
179      (let ((yy (cvode-yy solver)))
180        (let ((v (f64vector-ref yy 0))
181              (w (f64vector-ref yy 1)))
182          (printf "v = ~A w = ~A~%" v w)
183          ))
184
185      (cvode-destroy-solver solver)
186     
187      )))
188
189
190(define (print-results solver t)
191  (let ((yy (cvode-yy solver)))
192    (printf "~A ~A ~A ~A ~A ~A~%" 
193            t 
194            (f64vector-ref yy 0)
195            (f64vector-ref yy 1)
196            (cvode-get-last-order solver)
197            (cvode-get-num-steps solver)
198            (cvode-get-last-step solver)
199            )))
200     
201(cvode-main)
202;(cvode-main/unsafe)
203;(ida-main)     
Note: See TracBrowser for help on using the repository browser.