source: project/ode/trunk/euler.scm @ 8001

Last change on this file since 8001 was 8001, checked in by Ivan Raikov, 13 years ago

Bug fixes and changes to ode-hhsm and the solvers.

File size: 5.1 KB
Line 
1
2;;
3;; Euler's integration method.
4;;
5;;
6;; Copyright 2007 Ivan Raikov and the Okinawa Institute of Science and Technology
7;;
8;; This program is free software: you can redistribute it and/or
9;; modify it under the terms of the GNU General Public License as
10;; published by the Free Software Foundation, either version 3 of the
11;; License, or (at your option) any later version.
12;;
13;; This program is distributed in the hope that it will be useful, but
14;; WITHOUT ANY WARRANTY; without even the implied warranty of
15;; MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
16;; General Public License for more details.
17;;
18;; A full copy of the GPL license can be found at
19;; <http://www.gnu.org/licenses/>.
20;;
21
22(include "lists.scm")
23
24
25(define (t<tstop tstep t tstop)
26  (if (> tstep 0.0) (< t tstop) (> t tstop)))
27
28(define (make-ode:euler numerror make-vector vector-ref vector-set! vector?)
29
30 (lambda (ode-runtime indep tstart tstop dt initial-tstep  . rest)
31
32  (define guards (if (pair? rest) (car rest) (list)))
33  (define debug  (and (pair? rest) (and (pair? (cdr rest)) (second rest))))
34
35  (define maxerror   (ode-runtime 'maxerror))
36  (define hierror?   (ode-runtime 'hierror?))
37  (define lowerror?  (ode-runtime 'lowerror?))
38
39  (define printq  (ode-runtime 'print))
40  (define eval-rhs  (ode-runtime 'eval-rhs))
41  (define is-state? (ode-runtime 'is-state?))
42  (define is-asgn?  (ode-runtime 'is-asgn?))
43  (define solve-env-ref  (ode-runtime 'solve-env-ref))
44  (define solve-env-set! (ode-runtime 'solve-env-set!))
45  (define eval-env-ref  (ode-runtime 'eval-env-ref))
46  (define eval-env-set! (ode-runtime 'eval-env-set!))
47
48  (define svec-deriv-idx  (ode-runtime 'svec-deriv-idx))
49  (define svec-value-idx  (ode-runtime 'svec-value-idx))
50  (define svec-relerr-idx  (ode-runtime 'svec-relerr-idx))
51  (define svec-abserr-idx  (ode-runtime 'svec-abserr-idx))
52 
53  (define (svec-deriv svec) (vector-ref svec svec-deriv-idx))
54  (define (svec-value svec) (vector-ref svec svec-value-idx))
55  (define (svec-relerr svec) (vector-ref svec svec-relerr-idx))
56  (define (svec-abserr svec) (vector-ref svec svec-abserr-idx))
57
58  (define (svec-deriv-set! svec v)  (vector-set! svec svec-deriv-idx v))
59  (define (svec-value-set! svec v)  (vector-set! svec svec-value-idx v))
60  (define (svec-relerr-set! svec v) (vector-set! svec svec-relerr-idx v))
61  (define (svec-abserr-set! svec v) (vector-set! svec svec-abserr-idx v))
62
63  ;; For each state, this vector contains the value of hy' at the
64  ;; current step
65  (define hy1vect  (make-vector (ode-runtime 'order) 0.0))
66
67  ;; Compensated-summation vector: this is used to reduce accumulated
68  ;; rounding error
69  (define z (make-vector (ode-runtime 'order) 0.0))
70
71  ;; Reduce accumulated rounding error on the independent variable
72  (define tz 0)
73
74  ;; compute y' and hy' for each state
75  (define (eval-derivs! tstep)
76    ((ode-runtime 'eval-poset-fold)
77     (lambda (order i lst)
78       (fold
79        (lambda (x i)
80          (let ((xrec (if (> (length x) 2) x (cons #f x))))
81            (let ((df   (first xrec))
82                  (sym  (second xrec))
83                  (rhs  (third xrec)))
84            (let ((v     (eval-rhs sym rhs))
85                  (svec  (solve-env-ref sym)))
86              (cond ((is-state? sym)  (begin
87                                        (case df
88                                          ((q) (begin
89                                                 (vector-set! hy1vect i 0)
90                                                 (svec-value-set! svec v)))
91                                         
92                                          ((d) (let ((hy1 (+ (* tstep v) (vector-ref z i))))
93                                                 (svec-deriv-set! svec v)
94                                                 (vector-set! hy1vect i hy1))))
95                                        (fx+ 1 i)))
96                  ((is-asgn? sym)   (begin
97                                      (eval-env-set! sym v)
98                                      (if (vector? svec)
99                                          (svec-value-set! svec v)
100                                          (set-car! svec v))
101                                      i))
102                  (else  i)))))) i lst)) 0))
103
104 
105  ;; Standard forward Euler method
106  (define (euler startit t tstep)
107
108    (eval-env-set! indep t)
109    (eval-env-set! dt tstep)
110    (eval-derivs! tstep)
111
112    (let loop ((it startit) (t t))
113
114      (for-each (lambda (g)
115                  (if (not (eval-rhs 'guard g))
116                      (numerror "guard failed: " g))) guards)
117
118      ((ode-runtime 'solve-env-state-foreach) 
119       (lambda (sym svec)
120         (eval-env-set! sym (svec-value svec))))
121
122      ;; update independent variable
123      (solve-env-set! indep svec-value-idx t)
124
125      ;; print quantities
126      (printq)
127
128      (if (t<tstop tstep t tstop)
129          (begin
130            (eval-derivs! tstep)
131            ((ode-runtime 'eval-poset-fold) 
132             (lambda (order i lst)
133               (fold 
134                (lambda (x i)
135                  (let ((xrec (if (> (length x) 2) x (cons #f x))))
136                    (let ((df   (first xrec))
137                          (sym  (second xrec))
138                          (rhs  (third xrec)))
139                      (if df (let ((svec  (solve-env-ref sym))
140                                   (hy1  (vector-ref hy1vect i)))
141                               (if (not (zero? hy1))
142                                   (let* ((oldv (svec-value svec))
143                                          (v    (+ oldv hy1)))
144                                     (svec-value-set! svec v)
145                                     (vector-set! z i (- hy1 (- v  oldv)))))
146                               (fx+ 1 i))
147                          i)))) i lst)) 0)
148
149            (let* ((t1  (+ t tstep tz))
150                   (t2  (- t1 t)))
151              (set! tz (- t2 tstep))
152              (eval-env-set! indep t1)
153              (loop (fx+ 1 it) t1)))
154          (begin
155            (if debug (begin
156                        (print "euler: ")
157                        (for-each (lambda (x) (print x)) (ode-runtime 'debug))))
158            (list it t tstep (list))))))
159
160  (if debug (begin
161              (print "euler: ")
162              (for-each (lambda (x) (print x)) (ode-runtime 'debug))))
163  (euler 0 0.0 initial-tstep)))
164
Note: See TracBrowser for help on using the repository browser.