GNU Octave  4.0.0
A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
quad.cc
Go to the documentation of this file.
1 /*
2 
3 Copyright (C) 1996-2015 John W. Eaton
4 
5 This file is part of Octave.
6 
7 Octave is free software; you can redistribute it and/or modify it
8 under the terms of the GNU General Public License as published by the
9 Free Software Foundation; either version 3 of the License, or (at your
10 option) any later version.
11 
12 Octave is distributed in the hope that it will be useful, but WITHOUT
13 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
14 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
15 for more details.
16 
17 You should have received a copy of the GNU General Public License
18 along with Octave; see the file COPYING. If not, see
19 <http://www.gnu.org/licenses/>.
20 
21 */
22 
23 #ifdef HAVE_CONFIG_H
24 #include <config.h>
25 #endif
26 
27 #include <string>
28 
29 #include <iomanip>
30 #include <iostream>
31 
32 #include "Quad.h"
33 #include "lo-mappers.h"
34 
35 #include "defun.h"
36 #include "error.h"
37 #include "gripes.h"
38 #include "pager.h"
39 #include "oct-obj.h"
40 #include "ov-fcn.h"
41 #include "unwind-prot.h"
42 #include "utils.h"
43 #include "variables.h"
44 
45 #include "Quad-opts.cc"
46 
47 #if defined (quad)
48 #undef quad
49 #endif
50 
51 // Global pointer for user defined function required by quadrature functions.
53 
54 // Have we warned about imaginary values returned from user function?
55 static bool warned_imaginary = false;
56 
57 // Is this a recursive call?
58 static int call_depth = 0;
59 
60 double
62 {
63  double retval = 0.0;
64 
65  octave_value_list args;
66  args(0) = x;
67 
68  if (quad_fcn)
69  {
70  octave_value_list tmp = quad_fcn->do_multi_index_op (1, args);
71 
72  if (error_state)
73  {
74  quad_integration_error = 1; // FIXME
75  gripe_user_supplied_eval ("quad");
76  return retval;
77  }
78 
79  if (tmp.length () && tmp(0).is_defined ())
80  {
81  if (! warned_imaginary && tmp(0).is_complex_type ())
82  {
83  warning ("quad: ignoring imaginary part returned from user-supplied function");
84  warned_imaginary = true;
85  }
86 
87  retval = tmp(0).double_value ();
88 
89  if (error_state)
90  {
91  quad_integration_error = 1; // FIXME
92  gripe_user_supplied_eval ("quad");
93  }
94  }
95  else
96  {
97  quad_integration_error = 1; // FIXME
98  gripe_user_supplied_eval ("quad");
99  }
100  }
101 
102  return retval;
103 }
104 
105 float
107 {
108  float retval = 0.0;
109 
110  octave_value_list args;
111  args(0) = x;
112 
113  if (quad_fcn)
114  {
115  octave_value_list tmp = quad_fcn->do_multi_index_op (1, args);
116 
117  if (error_state)
118  {
119  quad_integration_error = 1; // FIXME
120  gripe_user_supplied_eval ("quad");
121  return retval;
122  }
123 
124  if (tmp.length () && tmp(0).is_defined ())
125  {
126  if (! warned_imaginary && tmp(0).is_complex_type ())
127  {
128  warning ("quad: ignoring imaginary part returned from user-supplied function");
129  warned_imaginary = true;
130  }
131 
132  retval = tmp(0).float_value ();
133 
134  if (error_state)
135  {
136  quad_integration_error = 1; // FIXME
137  gripe_user_supplied_eval ("quad");
138  }
139  }
140  else
141  {
142  quad_integration_error = 1; // FIXME
143  gripe_user_supplied_eval ("quad");
144  }
145  }
146 
147  return retval;
148 }
149 
150 #define QUAD_ABORT() \
151  do \
152  { \
153  if (fcn_name.length ()) \
154  clear_function (fcn_name); \
155  return retval; \
156  } \
157  while (0)
158 
159 #define QUAD_ABORT1(msg) \
160  do \
161  { \
162  ::error ("quad: " msg); \
163  QUAD_ABORT (); \
164  } \
165  while (0)
166 
167 #define QUAD_ABORT2(fmt, arg) \
168  do \
169  { \
170  ::error ("quad: " fmt, arg); \
171  QUAD_ABORT (); \
172  } \
173  while (0)
174 
175 DEFUN (quad, args, nargout,
176  "-*- texinfo -*-\n\
177 @deftypefn {Built-in Function} {@var{q} =} quad (@var{f}, @var{a}, @var{b})\n\
178 @deftypefnx {Built-in Function} {@var{q} =} quad (@var{f}, @var{a}, @var{b}, @var{tol})\n\
179 @deftypefnx {Built-in Function} {@var{q} =} quad (@var{f}, @var{a}, @var{b}, @var{tol}, @var{sing})\n\
180 @deftypefnx {Built-in Function} {[@var{q}, @var{ier}, @var{nfun}, @var{err}] =} quad (@dots{})\n\
181 Numerically evaluate the integral of @var{f} from @var{a} to @var{b} using\n\
182 Fortran routines from @w{@sc{quadpack}}.\n\
183 \n\
184 @var{f} is a function handle, inline function, or a string containing the\n\
185 name of the function to evaluate. The function must have the form @code{y =\n\
186 f (x)} where @var{y} and @var{x} are scalars.\n\
187 \n\
188 @var{a} and @var{b} are the lower and upper limits of integration. Either\n\
189 or both may be infinite.\n\
190 \n\
191 The optional argument @var{tol} is a vector that specifies the desired\n\
192 accuracy of the result. The first element of the vector is the desired\n\
193 absolute tolerance, and the second element is the desired relative\n\
194 tolerance. To choose a relative test only, set the absolute\n\
195 tolerance to zero. To choose an absolute test only, set the relative\n\
196 tolerance to zero. Both tolerances default to @code{sqrt (eps)} or\n\
197 approximately @math{1.5e^{-8}}.\n\
198 \n\
199 The optional argument @var{sing} is a vector of values at which the\n\
200 integrand is known to be singular.\n\
201 \n\
202 The result of the integration is returned in @var{q}.\n\
203 \n\
204 @var{ier} contains an integer error code (0 indicates a successful\n\
205 integration).\n\
206 \n\
207 @var{nfun} indicates the number of function evaluations that were\n\
208 made.\n\
209 \n\
210 @var{err} contains an estimate of the error in the solution.\n\
211 \n\
212 The function @code{quad_options} can set other optional parameters for\n\
213 @code{quad}.\n\
214 \n\
215 Note: because @code{quad} is written in Fortran it cannot be called\n\
216 recursively. This prevents its use in integrating over more than one\n\
217 variable by routines @code{dblquad} and @code{triplequad}.\n\
218 @seealso{quad_options, quadv, quadl, quadgk, quadcc, trapz, dblquad, triplequad}\n\
219 @end deftypefn")
220 {
221  octave_value_list retval;
222 
223  std::string fcn_name;
224 
225  warned_imaginary = false;
226 
227  unwind_protect frame;
228 
229  frame.protect_var (call_depth);
230  call_depth++;
231 
232  if (call_depth > 1)
233  QUAD_ABORT1 ("invalid recursive call");
234 
235  int nargin = args.length ();
236 
237  if (nargin > 2 && nargin < 6 && nargout < 5)
238  {
239  if (args(0).is_function_handle () || args(0).is_inline_function ())
240  quad_fcn = args(0).function_value ();
241  else
242  {
243  fcn_name = unique_symbol_name ("__quad_fcn__");
244  std::string fname = "function y = ";
245  fname.append (fcn_name);
246  fname.append ("(x) y = ");
247  quad_fcn = extract_function (args(0), "quad", fcn_name, fname,
248  "; endfunction");
249  }
250 
251  if (! quad_fcn)
252  QUAD_ABORT ();
253 
254  if (args(1).is_single_type () || args(2).is_single_type ())
255  {
256  float a = args(1).float_value ();
257 
258  if (error_state)
259  QUAD_ABORT1 ("expecting second argument to be a scalar");
260 
261  float b = args(2).float_value ();
262 
263  if (error_state)
264  QUAD_ABORT1 ("expecting third argument to be a scalar");
265 
266  int indefinite = 0;
269  float bound = 0.0;
270  if (xisinf (a) && xisinf (b))
271  {
272  indefinite = 1;
273  indef_type = FloatIndefQuad::doubly_infinite;
274  }
275  else if (xisinf (a))
276  {
277  indefinite = 1;
278  bound = b;
280  }
281  else if (xisinf (b))
282  {
283  indefinite = 1;
284  bound = a;
285  indef_type = FloatIndefQuad::bound_to_inf;
286  }
287 
288  octave_idx_type ier = 0;
289  octave_idx_type nfun = 0;
290  float abserr = 0.0;
291  float val = 0.0;
292  bool have_sing = false;
293  FloatColumnVector sing;
294  FloatColumnVector tol;
295 
296  switch (nargin)
297  {
298  case 5:
299  if (indefinite)
300  QUAD_ABORT1 ("singularities not allowed on infinite intervals");
301 
302  have_sing = true;
303 
304  sing = FloatColumnVector (args(4).float_vector_value ());
305 
306  if (error_state)
307  QUAD_ABORT1 ("expecting vector of singularities as fourth argument");
308 
309  case 4:
310  tol = FloatColumnVector (args(3).float_vector_value ());
311 
312  if (error_state)
313  QUAD_ABORT1 ("expecting vector of tolerances as fifth argument");
314 
315  switch (tol.capacity ())
316  {
317  case 2:
319 
320  case 1:
322  break;
323 
324  default:
325  QUAD_ABORT1 ("expecting tol to contain no more than two values");
326  }
327 
328  case 3:
329  if (indefinite)
330  {
332  indef_type);
333  iq.set_options (quad_opts);
334  val = iq.float_integrate (ier, nfun, abserr);
335  }
336  else
337  {
338  if (have_sing)
339  {
340  FloatDefQuad dq (quad_float_user_function, a, b, sing);
341  dq.set_options (quad_opts);
342  val = dq.float_integrate (ier, nfun, abserr);
343  }
344  else
345  {
347  dq.set_options (quad_opts);
348  val = dq.float_integrate (ier, nfun, abserr);
349  }
350  }
351  break;
352 
353  default:
354  panic_impossible ();
355  break;
356  }
357 
358  retval(3) = abserr;
359  retval(2) = nfun;
360  retval(1) = ier;
361  retval(0) = val;
362 
363  }
364  else
365  {
366  double a = args(1).double_value ();
367 
368  if (error_state)
369  QUAD_ABORT1 ("expecting second argument to be a scalar");
370 
371  double b = args(2).double_value ();
372 
373  if (error_state)
374  QUAD_ABORT1 ("expecting third argument to be a scalar");
375 
376  int indefinite = 0;
378  double bound = 0.0;
379  if (xisinf (a) && xisinf (b))
380  {
381  indefinite = 1;
382  indef_type = IndefQuad::doubly_infinite;
383  }
384  else if (xisinf (a))
385  {
386  indefinite = 1;
387  bound = b;
388  indef_type = IndefQuad::neg_inf_to_bound;
389  }
390  else if (xisinf (b))
391  {
392  indefinite = 1;
393  bound = a;
394  indef_type = IndefQuad::bound_to_inf;
395  }
396 
397  octave_idx_type ier = 0;
398  octave_idx_type nfun = 0;
399  double abserr = 0.0;
400  double val = 0.0;
401  bool have_sing = false;
402  ColumnVector sing;
403  ColumnVector tol;
404 
405  switch (nargin)
406  {
407  case 5:
408  if (indefinite)
409  QUAD_ABORT1 ("singularities not allowed on infinite intervals");
410 
411  have_sing = true;
412 
413  sing = ColumnVector (args(4).vector_value ());
414 
415  if (error_state)
416  QUAD_ABORT1 ("expecting vector of singularities as fourth argument");
417 
418  case 4:
419  tol = ColumnVector (args(3).vector_value ());
420 
421  if (error_state)
422  QUAD_ABORT1 ("expecting vector of tolerances as fifth argument");
423 
424  switch (tol.capacity ())
425  {
426  case 2:
428 
429  case 1:
431  break;
432 
433  default:
434  QUAD_ABORT1 ("expecting tol to contain no more than two values");
435  }
436 
437  case 3:
438  if (indefinite)
439  {
440  IndefQuad iq (quad_user_function, bound, indef_type);
441  iq.set_options (quad_opts);
442  val = iq.integrate (ier, nfun, abserr);
443  }
444  else
445  {
446  if (have_sing)
447  {
448  DefQuad dq (quad_user_function, a, b, sing);
449  dq.set_options (quad_opts);
450  val = dq.integrate (ier, nfun, abserr);
451  }
452  else
453  {
454  DefQuad dq (quad_user_function, a, b);
455  dq.set_options (quad_opts);
456  val = dq.integrate (ier, nfun, abserr);
457  }
458  }
459  break;
460 
461  default:
462  panic_impossible ();
463  break;
464  }
465 
466  retval(3) = abserr;
467  retval(2) = nfun;
468  retval(1) = ier;
469  retval(0) = val;
470  }
471 
472  if (fcn_name.length ())
473  clear_function (fcn_name);
474  }
475  else
476  print_usage ();
477 
478  return retval;
479 }
480 
481 /*
482 %!function y = __f (x)
483 %! y = x + 1;
484 %!endfunction
485 
486 %!test
487 %! [v, ier, nfun, err] = quad ("__f", 0, 5);
488 %! assert (ier, 0);
489 %! assert (v, 17.5, sqrt (eps));
490 %! assert (nfun > 0);
491 %! assert (err < sqrt (eps));
492 
493 %!test
494 %! [v, ier, nfun, err] = quad ("__f", single (0), single (5));
495 %! assert (ier, 0);
496 %! assert (v, 17.5, sqrt (eps ("single")));
497 %! assert (nfun > 0);
498 %! assert (err < sqrt (eps ("single")));
499 
500 %!function y = __f (x)
501 %! y = x .* sin (1 ./ x) .* sqrt (abs (1 - x));
502 %!endfunction
503 
504 %!test
505 %! [v, ier, nfun, err] = quad ("__f", 0.001, 3);
506 %! assert (ier == 0 || ier == 1);
507 %! assert (v, 1.98194120273598, sqrt (eps));
508 %! assert (nfun > 0);
509 
510 %!test
511 %! [v, ier, nfun, err] = quad ("__f", single (0.001), single (3));
512 %! assert (ier == 0 || ier == 1);
513 %! assert (v, 1.98194120273598, sqrt (eps ("single")));
514 %! assert (nfun > 0);
515 
516 %!error quad ()
517 %!error quad ("__f", 1, 2, 3, 4, 5)
518 
519 %!test
520 %! quad_options ("absolute tolerance", eps);
521 %! assert (quad_options ("absolute tolerance") == eps);
522 
523 %!error quad_options (1, 2, 3)
524 */
void set_single_precision_relative_tolerance(float val)
Definition: Quad-opts.h:80
virtual double integrate(void)
Definition: Quad.h:56
virtual octave_value_list do_multi_index_op(int nargout, const octave_value_list &idx)
Definition: ov-base.cc:203
#define QUAD_ABORT1(msg)
Definition: quad.cc:159
virtual float float_integrate(void)
Definition: Quad.h:63
OCTINTERP_API void print_usage(void)
Definition: defun.cc:51
void set_relative_tolerance(double val)
Definition: Quad-opts.h:74
octave_idx_type length(void) const
Definition: oct-obj.h:89
double quad_user_function(double x)
Definition: quad.cc:61
void protect_var(T &var)
#define DEFUN(name, args_name, nargout_name, doc)
Definition: defun.h:44
bool xisinf(double x)
Definition: lo-mappers.cc:160
IntegralType
Definition: Quad.h:163
static Quad_options quad_opts
Definition: Quad-opts.cc:20
std::string unique_symbol_name(const std::string &basename)
Definition: variables.cc:495
void clear_function(const std::string &nm)
Definition: variables.cc:77
void set_single_precision_absolute_tolerance(float val)
Definition: Quad-opts.h:77
float quad_float_user_function(float x)
Definition: quad.cc:106
#define QUAD_ABORT()
Definition: quad.cc:150
int error_state
Definition: error.cc:101
#define panic_impossible()
Definition: error.h:33
octave_idx_type capacity(void) const
Number of elements in the array.
Definition: Array.h:256
Definition: Quad.h:120
static bool warned_imaginary
Definition: quad.cc:55
void set_options(const Quad_options &opt)
Definition: Quad-opts.h:60
static octave_function * quad_fcn
Definition: quad.cc:52
void warning(const char *fmt,...)
Definition: error.cc:681
void set_absolute_tolerance(double val)
Definition: Quad-opts.h:71
octave_function * extract_function(const octave_value &arg, const std::string &warn_for, const std::string &fname, const std::string &header, const std::string &trailer)
Definition: variables.cc:140
static int call_depth
Definition: quad.cc:58
virtual octave_function * function_value(bool silent=false)
Definition: ov-base.cc:992
int quad_integration_error
Definition: Quad.cc:39
void gripe_user_supplied_eval(const char *name)
Definition: gripes.cc:87
F77_RET_T const double * x