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

Last change on this file since 29469 was 29469, checked in by Ivan Raikov, 7 years ago

sundials: removing assertions from examples with no closed-form solution

File size: 4.1 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  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-4
125                   reltol: 1e-4)))
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      (cvode-destroy-solver solver)
146     
147      )))
148
149(define (cvode-main/unsafe)
150 
151  (let ((yy (f64vector -60.899 0.0149)) ;; v w
152
153        ;; Integration limits
154        (t0  0.0)
155        (tf  TEND)
156        (dt  0.1))
157   
158    ;; CVODE initialization
159    (let ((solver (cvode-create-solver/unsafe
160                   t0 yy rhs/unsafe
161                   tstop: tf
162                   abstol: 1e-4
163                   reltol: 1e-4)))
164
165      ;; In loop, call CVodeSolve, print results, and test for error.
166
167      (let recur ((tnext (+ t0 dt)) (iout 1))
168
169        (let ((flag  (cvode-solve solver tnext)))
170          (if (negative? flag) (error 'main "CVODE solver error" flag))
171
172;;        (print-results solver tnext)
173
174          (if (< tnext tf)
175              (recur (+ tnext dt) (+ 1 iout)))
176          ))
177     
178      (let ((yy (cvode-yy solver)))
179        (let ((v (f64vector-ref yy 0))
180              (w (f64vector-ref yy 1)))
181          (printf "v = ~A w = ~A~%" v w)
182          ))
183
184      (cvode-destroy-solver solver)
185     
186      )))
187
188
189(define (print-results solver t)
190  (let ((yy (cvode-yy solver)))
191    (printf "~A ~A ~A ~A ~A ~A~%" 
192            t 
193            (f64vector-ref yy 0)
194            (f64vector-ref yy 1)
195            (cvode-get-last-order solver)
196            (cvode-get-num-steps solver)
197            (cvode-get-last-step solver)
198            )))
199     
200(ida-main)     
201(cvode-main)
202(cvode-main/unsafe)
Note: See TracBrowser for help on using the repository browser.