1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
|
;;; -*- Mode: Lisp; Package: Maxima; Syntax: Common-Lisp; Base: 10 -*- ;;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;; The data in this file contains enhancements. ;;;;;
;;; ;;;;;
;;; Copyright (c) 1984,1987 by William Schelter,University of Texas ;;;;;
;;; All rights reserved ;;;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;; (c) Copyright 1982 Massachusetts Institute of Technology ;;;
;;; Original code by CFFK. Modified to interface correctly with TRANSL ;;;
;;; and the rest of macsyma by GJC ;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
(in-package :maxima)
(macsyma-module rombrg)
(load-macsyma-macros transm)
(defmvar $rombergit 11 "the maximum number of iterations")
(defmvar $rombergmin 0 "the minimum number of iterations")
(defmvar $rombergtol 1e-4 "the relative tolerance of error")
(defmvar $rombergabs 0.0 "the absolute tolerance of error")
(defmvar $rombergit_used 0 "the number of iterations actually used.")
(defun romberg-subr (f left right)
(let (a b
(x 0.0)
(tt (make-array $rombergit :element-type 'flonum))
(rr (make-array $rombergit :element-type 'flonum)))
(let (($numer t) ($%enumer t))
(setq
a ($float left)
b ($float right)))
(if (not (and (numberp a) (numberp b)))
(return-from romberg-subr (values nil a b)))
(setq x (- b a))
(let
((fa (funcall f a))
(fb (funcall f b))
(fab (funcall f (* (+ b a) 0.5))))
(if
(or
(not (numberp fa))
(not (numberp fb))
(not (numberp fab)))
(return-from romberg-subr (values nil a b))
(progn
(setf (aref tt 0) (* x (+ fb fa) 0.5))
(setf (aref rr 0) (* x fab)))))
(do ((l 1 (1+ l))
(m 4 (* m 2))
(y 0.0)
(z 0.0)
(cerr 0.0))
((= l $rombergit)
(merror "`romberg' failed to converge"))
(declare (fixnum l m))
(setq y (float m) z (/ x y))
(setf (aref tt l) (* (+ (aref tt (1- l))
(aref rr (1- l))) 0.5))
(setf (aref rr l) 0.0)
(do ((i 1 (+ i 2)))
((> i m))
(let
((fzi+a (funcall f (+ (* z (float i)) a))))
(if (not (numberp fzi+a))
(return-from romberg-subr (values nil a b))
(setf (aref rr l) (+ fzi+a (aref rr l))))))
(setf (aref rr l) (* z (aref rr l) 2))
(setq y 0.0)
(do ((k l (1- k))) ((= k 0))
(declare (fixnum k))
(setq y (+ (* y 4) 3))
(setf (aref tt (1- k))
(+ (/ (- (aref tt k)
(aref tt (1- k))) y)
(aref tt k)))
(setf (aref rr (1- k))
(+ (/ (- (aref rr k) (aref rr (1- k))) y)
(aref rr k))))
(setq y (* (+ (aref tt 0) (aref rr 0)) 0.5))
;;; this is the WIN condition test.
(cond ((and (or (not (< $rombergabs
(setq cerr (abs (- (aref tt 0) (aref rr 0))))))
(not (< $rombergtol
;; cerr = "calculated error"; used for ^]
(setq cerr (/ cerr (if (zerop y) 1.0 (abs y)))))))
(> l $rombergmin))
(setq $rombergit_used l)
(return-from romberg-subr (values y nil nil)))))))
(defun $romberg (&rest args)
(cond
((= (length args) 3)
(multiple-value-bind
(result left right)
;; BIND EVIL SPECIAL VARIABLE *PLOT-REALPART* HERE ...
(let (($numer t) ($%enumer t) (*plot-realpart* nil))
(romberg-subr
(coerce-float-fun (first args))
(second args)
(third args)))
(if (numberp result)
result
`(($romberg) ,(first args) ,left ,right))))
((= (length args) 4)
(when ($constantp (second args))
(merror
"romberg: variable of integration cannot be a constant; found ~M"
(second args)))
(multiple-value-bind
(result left right)
;; BIND EVIL SPECIAL VARIABLE *PLOT-REALPART* HERE ...
(let (($numer t) ($%enumer t) (*plot-realpart* nil))
(romberg-subr
(coerce-float-fun (first args) `((mlist) ,(second args)))
(third args)
(fourth args)))
(if (numberp result)
result
`(($romberg) ,(first args) ,(second args) ,left ,right))))
(t
(wna-err '$romberg))))
|
