Mercurial > octave-nkf

comparison scripts/general/quadl.m @ 13008:85dac13a911b

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
quadl.m: Fix integration with large error tolerances (Bug #33792) quadl.m: Force recursion to occur at least once through global state variable. Miscellaneous spacing changes to improve code appearance.
author Rik <octave@nomad.inbox5.com>
date Thu, 25 Aug 2011 09:51:27 -0700
parents 412882f498b4
children e81ddf9cacd5
comparison
equal deleted inserted replaced
13006:61be447052c3 13008:85dac13a911b
60 ## 2004-02-10 Paul Kienzle 60 ## 2004-02-10 Paul Kienzle
61 ## * renamed to quadl for compatibility 61 ## * renamed to quadl for compatibility
62 ## * replace global variable terminate2 with local function need_warning 62 ## * replace global variable terminate2 with local function need_warning
63 ## * add paper ref to docs 63 ## * add paper ref to docs
64 64
65 function q = quadl (f, a, b, tol, trace, varargin) 65 function q = quadl (f, a, b, tol = [], trace = false, varargin)
66 need_warning (1); 66
67 if (nargin < 4) 67 if (nargin < 3)
68 tol = []; 68 print_usage ();
69 endif 69 endif
70 if (nargin < 5) 70
71 trace = [];
72 endif
73 if (isa (a, "single") || isa (b, "single")) 71 if (isa (a, "single") || isa (b, "single"))
74 myeps = eps ("single"); 72 myeps = eps ("single");
75 else 73 else
76 myeps = eps; 74 myeps = eps;
77 endif 75 endif
78 if (isempty (tol)) 76 if (isempty (tol))
79 tol = myeps; 77 tol = myeps;
80 endif 78 endif
81 if (isempty (trace)) 79 if (isempty (trace))
82 trace = 0; 80 trace = false;
83 endif 81 endif
84 if (tol < myeps) 82 if (tol < myeps)
85 tol = myeps; 83 tol = myeps;
86 endif 84 endif
87 85
86 ## Track whether recursion has occurred
87 global __quadl_recurse_done__;
88 __quadl_recurse_done__ = false;
89 ## Track whether warning about machine precision has been issued
90 global __quadl_need_warning__;
91 __quadl_need_warning__ = true;
92
88 m = (a+b)/2; 93 m = (a+b)/2;
89 h = (b-a)/2; 94 h = (b-a)/2;
90 alpha = sqrt(2/3); 95 alpha = sqrt (2/3);
91 beta = 1/sqrt(5); 96 beta = 1/sqrt (5);
92 97
93 x1 = .942882415695480; 98 x1 = .942882415695480;
94 x2 = .641853342345781; 99 x2 = .641853342345781;
95 x3 = .236383199662150; 100 x3 = .236383199662150;
96 101
102 fa = y(1); 107 fa = y(1);
103 fb = y(13); 108 fb = y(13);
104 109
105 i2 = (h/6)*(y(1) + y(13) + 5*(y(5)+y(9))); 110 i2 = (h/6)*(y(1) + y(13) + 5*(y(5)+y(9)));
106 111
107 i1 = (h/1470)*(77*(y(1)+y(13)) 112 i1 = (h/1470)*( 77*(y(1)+y(13))
108 + 432*(y(3)+y(11)) 113 + 432*(y(3)+y(11))
109 + 625*(y(5)+y(9)) 114 + 625*(y(5)+y(9))
110 + 672*y(7)); 115 + 672*y(7));
111 116
112 is = h*(.0158271919734802*(y(1)+y(13)) 117 is = h*( .0158271919734802*(y(1)+y(13))
113 +.0942738402188500*(y(2)+y(12)) 118 +.0942738402188500*(y(2)+y(12))
114 + .155071987336585*(y(3)+y(11)) 119 + .155071987336585*(y(3)+y(11))
115 + .188821573960182*(y(4)+y(10)) 120 + .188821573960182*(y(4)+y(10))
116 + .199773405226859*(y(5)+y(9)) 121 + .199773405226859*(y(5)+y(9))
117 + .224926465333340*(y(6)+y(8)) 122 + .224926465333340*(y(6)+y(8))
118 + .242611071901408*y(7)); 123 + .242611071901408*y(7));
119 124
120 s = sign(is); 125 s = sign (is);
121
122 if (s == 0) 126 if (s == 0)
123 s = 1; 127 s = 1;
124 endif 128 endif
125 erri1 = abs(i1-is); 129 erri1 = abs (i1-is);
126 erri2 = abs(i2-is); 130 erri2 = abs (i2-is);
127 R = 1;
128 if (erri2 != 0) 131 if (erri2 != 0)
129 R = erri1/erri2; 132 R = erri1/erri2;
133 else
134 R = 1;
130 endif 135 endif
131 if (R > 0 && R < 1) 136 if (R > 0 && R < 1)
132 tol = tol/R; 137 tol = tol/R;
133 endif 138 endif
134 is = s*abs(is)*tol/myeps; 139 is = s * abs(is) * tol/myeps;
135 if (is == 0) 140 if (is == 0)
136 is = b-a; 141 is = b-a;
137 endif 142 endif
143
138 q = adaptlobstp (f, a, b, fa, fb, is, trace, varargin{:}); 144 q = adaptlobstp (f, a, b, fa, fb, is, trace, varargin{:});
145
139 endfunction 146 endfunction
140 147
141 ## ADAPTLOBSTP Recursive function used by QUADL. 148 ## ADAPTLOBSTP Recursive function used by QUADL.
142 ## 149 ##
143 ## Q = ADAPTLOBSTP('F', A, B, FA, FB, IS, TRACE) tries to 150 ## Q = ADAPTLOBSTP('F', A, B, FA, FB, IS, TRACE) tries to
144 ## approximate the integral of F(X) from A to B to 151 ## approximate the integral of F(X) from A to B to
145 ## an appropriate relative error. The argument 'F' is 152 ## an appropriate relative error. The argument 'F' is
146 ## a string containing the name of f. The remaining 153 ## a string containing the name of f. The remaining
147 ## arguments are generated by ADAPTLOB or by recursion. 154 ## arguments are generated by ADAPTLOB or by recursion.
148 ## 155 ##
149 ## Walter Gautschi, 08/03/98 156 ## Walter Gautschi, 08/03/98
150 157
151 function q = adaptlobstp (f, a, b, fa, fb, is, trace, varargin) 158 function q = adaptlobstp (f, a, b, fa, fb, is, trace, varargin)
159 global __quadl_recurse_done__;
160 global __quadl_need_warning__;
161
152 h = (b-a)/2; 162 h = (b-a)/2;
153 m = (a+b)/2; 163 m = (a+b)/2;
154 alpha = sqrt(2/3); 164 alpha = sqrt (2/3);
155 beta = 1/sqrt(5); 165 beta = 1 / sqrt(5);
156 mll = m-alpha*h; 166 mll = m-alpha*h;
157 ml = m-beta*h; 167 ml = m-beta*h;
158 mr = m+beta*h; 168 mr = m+beta*h;
159 mrr = m+alpha*h; 169 mrr = m+alpha*h;
160 x = [mll, ml, m, mr, mrr]; 170 x = [mll, ml, m, mr, mrr];
161 y = feval(f, x, varargin{:}); 171 y = feval (f, x, varargin{:});
162 fmll = y(1); 172 fmll = y(1);
163 fml = y(2); 173 fml = y(2);
164 fm = y(3); 174 fm = y(3);
165 fmr = y(4); 175 fmr = y(4);
166 fmrr = y(5); 176 fmrr = y(5);
167 i2 = (h/6)*(fa + fb + 5*(fml+fmr)); 177 i2 = (h/6)*(fa + fb + 5*(fml+fmr));
168 i1 = (h/1470)*(77*(fa+fb) + 432*(fmll+fmrr) + 625*(fml+fmr) + 672*fm); 178 i1 = (h/1470)*(77*(fa+fb) + 432*(fmll+fmrr) + 625*(fml+fmr) + 672*fm);
169 if (is+(i1-i2) == is || mll <= a || b <= mrr) 179 if ((is+(i1-i2) == is || mll <= a || b <= mrr) && __quadl_recurse_done__)
170 if ((m <= a || b <= m) && need_warning ()) 180 if ((m <= a || b <= m) && __quadl_need_warning__)
171 warning ("quadl: interval contains no more machine number"); 181 warning ("quadl: interval contains no more machine number");
172 warning ("quadl: required tolerance may not be met"); 182 warning ("quadl: required tolerance may not be met");
173 need_warning (0); 183 __quadl_need_warning__ = false;
174 endif 184 endif
175 q = i1; 185 q = i1;
176 if (trace) 186 if (trace)
177 disp ([a, b-a, q]); 187 disp ([a, b-a, q]);
178 endif 188 endif
179 else 189 else
180 q = (adaptlobstp (f, a, mll, fa, fmll, is, trace, varargin{:}) 190 __quadl_recurse_done__ = true;
181 + adaptlobstp (f, mll, ml, fmll, fml, is, trace, varargin{:}) 191 q = ( adaptlobstp (f, a , mll, fa , fmll, is, trace, varargin{:})
182 + adaptlobstp (f, ml, m, fml, fm, is, trace, varargin{:}) 192 + adaptlobstp (f, mll, ml , fmll, fml , is, trace, varargin{:})
183 + adaptlobstp (f, m, mr, fm, fmr, is, trace, varargin{:}) 193 + adaptlobstp (f, ml , m , fml , fm , is, trace, varargin{:})
184 + adaptlobstp (f, mr, mrr, fmr, fmrr, is, trace, varargin{:}) 194 + adaptlobstp (f, m , mr , fm , fmr , is, trace, varargin{:})
185 + adaptlobstp (f, mrr, b, fmrr, fb, is, trace, varargin{:})); 195 + adaptlobstp (f, mr , mrr, fmr , fmrr, is, trace, varargin{:})
196 + adaptlobstp (f, mrr, b , fmrr, fb , is, trace, varargin{:}));
186 endif 197 endif
187 endfunction 198 endfunction
188 199
189 function r = need_warning (v)
190 persistent w = [];
191 if (nargin == 0)
192 r = w;
193 else
194 w = v;
195 endif
196 endfunction
197
198 200
199 ## basic functionality 201 ## basic functionality
200 %!assert( quadl (@(x) sin (x), 0, pi, [], []), 2, -3e-16) 202 %!assert (quadl (@(x) sin (x), 0, pi, [], []), 2, -3e-16)
201 203
202 ## the values here are very high so it may be unavoidable that this fails 204 ## the values here are very high so it may be unavoidable that this fails
203 %!assert ( quadl (@(x) sin (3*x).*cosh (x).*sinh (x),10,15), 205 %!assert (quadl (@(x) sin (3*x).*cosh (x).*sinh (x),10,15),
204 %! 2.588424538641647e+10, -9e-15) 206 %! 2.588424538641647e+10, -9e-15)
205 207
206 ## extra parameters 208 ## extra parameters
207 %!assert (quadl (@(x,a,b) sin (a + b*x), 0, 1, [], [], 2, 3), 209 %!assert (quadl (@(x,a,b) sin (a + b*x), 0, 1, [], [], 2, 3),
208 %! cos(2)/3 - cos(5)/3, - 3e-16) 210 %! cos(2)/3 - cos(5)/3, - 3e-16)
209 211
210 ## test different tolerances. This test currently fails for a very high 212 ## test different tolerances.
211 ## tolerances. 213 %!assert (quadl (@(x) sin (2 + 3*x).^2, 0, 10, 0.3, []),
212 %!assert ( quadl (@(x) sin (2 + 3*x).^2, 0, 10, 0.3, []),
213 %! (60 + sin(4) - sin(64))/12, -0.3) 214 %! (60 + sin(4) - sin(64))/12, -0.3)
214 215 %!assert (quadl (@(x) sin (2 + 3*x).^2, 0, 10, 0.1, []),
215
216 ## for lower tolerances the test passes.
217 %!assert ( quadl (@(x) sin (2 + 3*x).^2, 0, 10, 0.1, []),
218 %! (60 + sin(4) - sin(64))/12, -0.1) 216 %! (60 + sin(4) - sin(64))/12, -0.1)
217