source: project/ode/trunk/extensions/stalin-numerr.scm @ 8365

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

Added stalin-related extension.

File size: 1.9 KB
Line 
1
2;; (define-record errinfo absmax absname relmax relname)
3
4(define make-errinfo     (primitive-procedure make-structure errinfo 4))
5(define errinfo-absmax   (primitive-procedure structure-ref errinfo 0))
6(define errinfo-absname  (primitive-procedure structure-ref errinfo 1))
7(define errinfo-relmax   (primitive-procedure structure-ref errinfo 2))
8(define errinfo-relname  (primitive-procedure structure-ref errinfo 3))
9
10(define (maxerror)
11  (let ((erelmax 0.0) (eabsmax 0.0)
12        (relname #f)  (absname #f))
13    (let ((proc
14           (lambda (q)
15             (let ((svec (vector-ref solve-env (qsym-i q))))
16               (let ((relerr (vector-ref svec svec-relerr-idx))
17                     (abserr (vector-ref svec svec-abserr-idx)))
18                 (cond ((< erelmax relerr)
19                        (set! erelmax relerr)
20                        (set! relname q))
21                       ((< eabsmax abserr)
22                        (set! eabsmax abserr)
23                        (set! absname q))))))))
24      (for-each proc states)
25      (for-each proc asgns)
26      (make-errinfo eabsmax absname erelmax relname ))))
27
28(define (hierror? relmax absmax hmin)
29  (lambda (t tstep err)
30    (let ((relname (errinfo-relname err))
31          (absname (errinfo-absname err))
32          (erelmax (errinfo-relmax err))
33          (eabsmax (errinfo-absmax err)))
34      (if (= t (+ tstep t))
35          (error "step size is smaller than machine precision: " "t = " t " tstep = " tstep)
36          (cond  ((and (<= erelmax relmax) (<= eabsmax absmax)) #f)
37                 ((<= (abs hmin) (abs tstep)) (max erelmax eabsmax))
38                 ((< relmax erelmax) 
39                  (error "relative error limit exceeded while calculating " relname))
40                 ((< absmax eabsmax)
41                  (error "absolute error limit exceeded while calculating " absname
42                         "; absmax = " absmax))
43                 (else #f))))))
44
45(define (lowerror? relmin absmin hmax)
46  (lambda (t tstep err)
47    (let ((relname (errinfo-relname err))
48          (absname (errinfo-absname err))
49          (erelmax (errinfo-relmax err))
50          (eabsmax (errinfo-absmax err)))
51      (if (or (< erelmax (* 0.25 relmin)) (< eabsmax (* 0.25 absmin)))
52          (and (<= (abs tstep) (abs hmax)) (max erelmax eabsmax))
53          #f))))
Note: See TracBrowser for help on using the repository browser.