1 ////////////////////////////////////////////////////////////////////////
2 //
3 // Copyright (C) 1996-2021 The Octave Project Developers
4 //
5 // See the file COPYRIGHT.md in the top-level directory of this
6 // distribution or <https://octave.org/copyright/>.
7 //
8 // This file is part of Octave.
9 //
10 // Octave is free software: you can redistribute it and/or modify it
11 // under the terms of the GNU General Public License as published by
12 // the Free Software Foundation, either version 3 of the License, or
13 // (at your option) any later version.
14 //
15 // Octave is distributed in the hope that it will be useful, but
16 // WITHOUT ANY WARRANTY; without even the implied warranty of
17 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
18 // GNU General Public License for more details.
19 //
20 // You should have received a copy of the GNU General Public License
21 // along with Octave; see the file COPYING. If not, see
22 // <https://www.gnu.org/licenses/>.
23 //
24 ////////////////////////////////////////////////////////////////////////
25
26 #if defined (HAVE_CONFIG_H)
27 # include "config.h"
28 #endif
29
30 #include <string>
31
32 #include "Quad.h"
33 #include "lo-mappers.h"
34
35 #include "defun.h"
36 #include "error.h"
37 #include "errwarn.h"
38 #include "interpreter-private.h"
39 #include "pager.h"
40 #include "parse.h"
41 #include "ov.h"
42 #include "ovl.h"
43 #include "unwind-prot.h"
44 #include "utils.h"
45 #include "variables.h"
46
47 #include "Quad-opts.cc"
48
49 // Global pointer for user defined function required by quadrature functions.
50 static octave_value quad_fcn;
51
52 // Have we warned about imaginary values returned from user function?
53 static bool warned_imaginary = false;
54
55 // Is this a recursive call?
56 static int call_depth = 0;
57
58 double
quad_user_function(double x)59 quad_user_function (double x)
60 {
61 double retval = 0.0;
62
63 octave_value_list args;
64 args(0) = x;
65
66 if (quad_fcn.is_defined ())
67 {
68 octave_value_list tmp;
69
70 try
71 {
72 tmp = octave::feval (quad_fcn, args, 1);
73 }
74 catch (octave::execution_exception& e)
75 {
76 err_user_supplied_eval (e, "quad");
77 }
78
79 if (! tmp.length () || ! tmp(0).is_defined ())
80 err_user_supplied_eval ("quad");
81
82 if (! warned_imaginary && tmp(0).iscomplex ())
83 {
84 warning ("quad: ignoring imaginary part returned from user-supplied function");
85 warned_imaginary = true;
86 }
87
88 retval = tmp(0).xdouble_value ("quad: expecting user supplied function to return numeric value");
89 }
90
91 return retval;
92 }
93
94 float
quad_float_user_function(float x)95 quad_float_user_function (float x)
96 {
97 float retval = 0.0;
98
99 octave_value_list args;
100 args(0) = x;
101
102 if (quad_fcn.is_defined ())
103 {
104 octave_value_list tmp;
105
106 try
107 {
108 tmp = octave::feval (quad_fcn, args, 1);
109 }
110 catch (octave::execution_exception& e)
111 {
112 err_user_supplied_eval (e, "quad");
113 }
114
115 if (! tmp.length () || ! tmp(0).is_defined ())
116 err_user_supplied_eval ("quad");
117
118 if (! warned_imaginary && tmp(0).iscomplex ())
119 {
120 warning ("quad: ignoring imaginary part returned from user-supplied function");
121 warned_imaginary = true;
122 }
123
124 retval = tmp(0).xfloat_value ("quad: expecting user supplied function to return numeric value");
125 }
126
127 return retval;
128 }
129
130 DEFMETHODX ("quad", Fquad, interp, args, ,
131 doc: /* -*- texinfo -*-
132 @deftypefn {} {@var{q} =} quad (@var{f}, @var{a}, @var{b})
133 @deftypefnx {} {@var{q} =} quad (@var{f}, @var{a}, @var{b}, @var{tol})
134 @deftypefnx {} {@var{q} =} quad (@var{f}, @var{a}, @var{b}, @var{tol}, @var{sing})
135 @deftypefnx {} {[@var{q}, @var{ier}, @var{nfun}, @var{err}] =} quad (@dots{})
136 Numerically evaluate the integral of @var{f} from @var{a} to @var{b} using
137 Fortran routines from @w{@sc{quadpack}}.
138
139 @var{f} is a function handle, inline function, or a string containing the
140 name of the function to evaluate. The function must have the form @code{y =
141 f (x)} where @var{y} and @var{x} are scalars.
142
143 @var{a} and @var{b} are the lower and upper limits of integration. Either
144 or both may be infinite.
145
146 The optional argument @var{tol} is a vector that specifies the desired
147 accuracy of the result. The first element of the vector is the desired
148 absolute tolerance, and the second element is the desired relative
149 tolerance. To choose a relative test only, set the absolute
150 tolerance to zero. To choose an absolute test only, set the relative
151 tolerance to zero. Both tolerances default to @code{sqrt (eps)} or
152 approximately 1.5e-8.
153
154 The optional argument @var{sing} is a vector of values at which the
155 integrand is known to be singular.
156
157 The result of the integration is returned in @var{q}.
158
159 @var{ier} contains an integer error code (0 indicates a successful
160 integration).
161
162 @var{nfun} indicates the number of function evaluations that were
163 made.
164
165 @var{err} contains an estimate of the error in the solution.
166
167 The function @code{quad_options} can set other optional parameters for
168 @code{quad}.
169
170 Note: because @code{quad} is written in Fortran it cannot be called
171 recursively. This prevents its use in integrating over more than one
172 variable by routines @code{dblquad} and @code{triplequad}.
173 @seealso{quad_options, quadv, quadl, quadgk, quadcc, trapz, dblquad, triplequad}
174 @end deftypefn */)
175 {
176 int nargin = args.length ();
177
178 if (nargin < 3 || nargin > 5)
179 print_usage ();
180
181 warned_imaginary = false;
182
183 octave::unwind_protect frame;
184
185 frame.protect_var (call_depth);
186 call_depth++;
187
188 if (call_depth > 1)
189 error ("quad: invalid recursive call");
190
191 quad_fcn = octave::get_function_handle (interp, args(0), "x");
192
193 octave_value_list retval;
194
195 if (args(1).is_single_type () || args(2).is_single_type ())
196 {
197 float a = args(1).xfloat_value ("quad: lower limit of integration A must be a scalar");
198 float b = args(2).xfloat_value ("quad: upper limit of integration B must be a scalar");
199
200 int indefinite = 0;
201 FloatIndefQuad::IntegralType indef_type
202 = FloatIndefQuad::doubly_infinite;
203 float bound = 0.0;
204 if (octave::math::isinf (a) && octave::math::isinf (b))
205 {
206 indefinite = 1;
207 indef_type = FloatIndefQuad::doubly_infinite;
208 }
209 else if (octave::math::isinf (a))
210 {
211 indefinite = 1;
212 bound = b;
213 indef_type = FloatIndefQuad::neg_inf_to_bound;
214 }
215 else if (octave::math::isinf (b))
216 {
217 indefinite = 1;
218 bound = a;
219 indef_type = FloatIndefQuad::bound_to_inf;
220 }
221
222 octave_idx_type ier = 0;
223 octave_idx_type nfun = 0;
224 float abserr = 0.0;
225 float val = 0.0;
226 bool have_sing = false;
227 FloatColumnVector sing;
228 FloatColumnVector tol;
229
230 switch (nargin)
231 {
232 case 5:
233 if (indefinite)
234 error ("quad: singularities not allowed on infinite intervals");
235
236 have_sing = true;
237
238 sing = args(4).xfloat_vector_value ("quad: fifth argument SING must be a vector of singularities");
239 OCTAVE_FALLTHROUGH;
240
241 case 4:
242 tol = args(3).xfloat_vector_value ("quad: TOL must be a 1 or 2-element vector");
243
244 switch (tol.numel ())
245 {
246 case 2:
247 quad_opts.set_single_precision_relative_tolerance (tol (1));
248 OCTAVE_FALLTHROUGH;
249
250 case 1:
251 quad_opts.set_single_precision_absolute_tolerance (tol (0));
252 break;
253
254 default:
255 error ("quad: TOL must be a 1 or 2-element vector");
256 }
257 OCTAVE_FALLTHROUGH;
258
259 case 3:
260 if (indefinite)
261 {
262 FloatIndefQuad iq (quad_float_user_function, bound,
263 indef_type);
264 iq.set_options (quad_opts);
265 val = iq.float_integrate (ier, nfun, abserr);
266 }
267 else
268 {
269 if (have_sing)
270 {
271 FloatDefQuad dq (quad_float_user_function, a, b, sing);
272 dq.set_options (quad_opts);
273 val = dq.float_integrate (ier, nfun, abserr);
274 }
275 else
276 {
277 FloatDefQuad dq (quad_float_user_function, a, b);
278 dq.set_options (quad_opts);
279 val = dq.float_integrate (ier, nfun, abserr);
280 }
281 }
282 break;
283
284 default:
285 panic_impossible ();
286 break;
287 }
288
289 retval = ovl (val, ier, nfun, abserr);
290
291 }
292 else
293 {
294 double a = args(1).xdouble_value ("quad: lower limit of integration A must be a scalar");
295 double b = args(2).xdouble_value ("quad: upper limit of integration B must be a scalar");
296
297 int indefinite = 0;
298 IndefQuad::IntegralType indef_type = IndefQuad::doubly_infinite;
299 double bound = 0.0;
300 if (octave::math::isinf (a) && octave::math::isinf (b))
301 {
302 indefinite = 1;
303 indef_type = IndefQuad::doubly_infinite;
304 }
305 else if (octave::math::isinf (a))
306 {
307 indefinite = 1;
308 bound = b;
309 indef_type = IndefQuad::neg_inf_to_bound;
310 }
311 else if (octave::math::isinf (b))
312 {
313 indefinite = 1;
314 bound = a;
315 indef_type = IndefQuad::bound_to_inf;
316 }
317
318 octave_idx_type ier = 0;
319 octave_idx_type nfun = 0;
320 double abserr = 0.0;
321 double val = 0.0;
322 bool have_sing = false;
323 ColumnVector sing;
324 ColumnVector tol;
325
326 switch (nargin)
327 {
328 case 5:
329 if (indefinite)
330 error ("quad: singularities not allowed on infinite intervals");
331
332 have_sing = true;
333
334 sing = args(4).xvector_value ("quad: fifth argument SING must be a vector of singularities");
335 OCTAVE_FALLTHROUGH;
336
337 case 4:
338 tol = args(3).xvector_value ("quad: TOL must be a 1 or 2-element vector");
339
340 switch (tol.numel ())
341 {
342 case 2:
343 quad_opts.set_relative_tolerance (tol (1));
344 OCTAVE_FALLTHROUGH;
345
346 case 1:
347 quad_opts.set_absolute_tolerance (tol (0));
348 break;
349
350 default:
351 error ("quad: TOL must be a 1 or 2-element vector");
352 }
353 OCTAVE_FALLTHROUGH;
354
355 case 3:
356 if (indefinite)
357 {
358 IndefQuad iq (quad_user_function, bound, indef_type);
359 iq.set_options (quad_opts);
360 val = iq.integrate (ier, nfun, abserr);
361 }
362 else
363 {
364 if (have_sing)
365 {
366 DefQuad dq (quad_user_function, a, b, sing);
367 dq.set_options (quad_opts);
368 val = dq.integrate (ier, nfun, abserr);
369 }
370 else
371 {
372 DefQuad dq (quad_user_function, a, b);
373 dq.set_options (quad_opts);
374 val = dq.integrate (ier, nfun, abserr);
375 }
376 }
377 break;
378
379 default:
380 panic_impossible ();
381 break;
382 }
383
384 retval = ovl (val, ier, nfun, abserr);
385 }
386
387 return retval;
388 }
389
390 /*
391 %!function y = __f (x)
392 %! y = x + 1;
393 %!endfunction
394
395 %!test
396 %! [v, ier, nfun, err] = quad ("__f", 0, 5);
397 %! assert (ier, 0);
398 %! assert (v, 17.5, sqrt (eps));
399 %! assert (nfun > 0);
400 %! assert (err < sqrt (eps));
401
402 %!test
403 %! [v, ier, nfun, err] = quad ("__f", single (0), single (5));
404 %! assert (ier, 0);
405 %! assert (v, 17.5, sqrt (eps ("single")));
406 %! assert (nfun > 0);
407 %! assert (err < sqrt (eps ("single")));
408
409 %!function y = __f (x)
410 %! y = x .* sin (1 ./ x) .* sqrt (abs (1 - x));
411 %!endfunction
412
413 %!test
414 %! [v, ier, nfun, err] = quad ("__f", 0.001, 3);
415 %! assert (ier == 0 || ier == 1);
416 %! assert (v, 1.98194120273598, sqrt (eps));
417 %! assert (nfun > 0);
418
419 %!test
420 %! [v, ier, nfun, err] = quad (@__f, 0.001, 3);
421 %! assert (ier == 0 || ier == 1);
422 %! assert (v, 1.98194120273598, sqrt (eps));
423 %! assert (nfun > 0);
424
425 %!test
426 %! fstr = "x .* sin (1 ./ x) .* sqrt (abs (1 - x))";
427 %! [v, ier, nfun, err] = quad (fstr, 0.001, 3);
428 %! assert (ier == 0 || ier == 1);
429 %! assert (v, 1.98194120273598, sqrt (eps));
430 %! assert (nfun > 0);
431
432 %!test
433 %! anon_fcn = @(x) x .* sin (1 ./ x) .* sqrt (abs (1 - x));
434 %! [v, ier, nfun, err] = quad (anon_fcn, 0.001, 3);
435 %! assert (ier == 0 || ier == 1);
436 %! assert (v, 1.98194120273598, sqrt (eps));
437 %! assert (nfun > 0);
438
439 %!test
440 %! inline_fcn = inline ("x .* sin (1 ./ x) .* sqrt (abs (1 - x))", "x");
441 %! [v, ier, nfun, err] = quad (inline_fcn, 0.001, 3);
442 %! assert (ier == 0 || ier == 1);
443 %! assert (v, 1.98194120273598, sqrt (eps));
444 %! assert (nfun > 0);
445
446 %!test
447 %! [v, ier, nfun, err] = quad ("__f", single (0.001), single (3));
448 %! assert (ier == 0 || ier == 1);
449 %! assert (v, 1.98194120273598, sqrt (eps ("single")));
450 %! assert (nfun > 0);
451
452 %!error quad ()
453 %!error quad ("__f", 1, 2, 3, 4, 5)
454
455 %!test
456 %! quad_options ("absolute tolerance", eps);
457 %! assert (quad_options ("absolute tolerance") == eps);
458
459 %!error quad_options (1, 2, 3)
460 */
461