Mercurial > octave
annotate src/DLD-FUNCTIONS/quad.cc @ 7562:c827f5673321
move tests to individual source files
| author | John W. Eaton <jwe@octave.org> |
|---|---|
| date | Thu, 06 Mar 2008 02:27:55 -0500 |
| parents | a1dbe9d80eee |
| children | 62affb34e648 |
| rev | line source |
|---|---|
| 2928 | 1 /* |
| 2 | |
| 7017 | 3 Copyright (C) 1996, 1997, 1998, 1999, 2000, 2002, 2004, 2005, 2006, |
| 4 2007 John W. Eaton | |
| 2928 | 5 |
| 6 This file is part of Octave. | |
| 7 | |
| 8 Octave is free software; you can redistribute it and/or modify it | |
| 9 under the terms of the GNU General Public License as published by the | |
| 7016 | 10 Free Software Foundation; either version 3 of the License, or (at your |
| 11 option) any later version. | |
| 2928 | 12 |
| 13 Octave is distributed in the hope that it will be useful, but WITHOUT | |
| 14 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or | |
| 15 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License | |
| 16 for more details. | |
| 17 | |
| 18 You should have received a copy of the GNU General Public License | |
| 7016 | 19 along with Octave; see the file COPYING. If not, see |
| 20 <http://www.gnu.org/licenses/>. | |
| 2928 | 21 |
| 22 */ | |
| 23 | |
| 24 #ifdef HAVE_CONFIG_H | |
| 25 #include <config.h> | |
| 26 #endif | |
| 27 | |
| 28 #include <string> | |
| 29 | |
| 3567 | 30 #include <iomanip> |
| 3523 | 31 #include <iostream> |
| 2928 | 32 |
| 33 #include "Quad.h" | |
| 34 #include "lo-mappers.h" | |
| 35 | |
| 36 #include "defun-dld.h" | |
| 37 #include "error.h" | |
| 38 #include "gripes.h" | |
| 39 #include "pager.h" | |
| 40 #include "oct-obj.h" | |
| 2968 | 41 #include "ov-fcn.h" |
| 3243 | 42 #include "unwind-prot.h" |
| 2928 | 43 #include "utils.h" |
| 44 #include "variables.h" | |
| 45 | |
| 3998 | 46 #include "Quad-opts.cc" |
| 47 | |
| 2928 | 48 #if defined (quad) |
| 49 #undef quad | |
| 50 #endif | |
| 51 | |
| 52 // Global pointer for user defined function required by quadrature functions. | |
| 2968 | 53 static octave_function *quad_fcn; |
| 2928 | 54 |
| 4140 | 55 // Have we warned about imaginary values returned from user function? |
| 56 static bool warned_imaginary = false; | |
| 57 | |
| 3243 | 58 // Is this a recursive call? |
| 59 static int call_depth = 0; | |
| 60 | |
| 2928 | 61 double |
| 62 quad_user_function (double x) | |
| 63 { | |
| 64 double retval = 0.0; | |
| 65 | |
| 66 octave_value_list args; | |
| 67 args(0) = x; | |
| 68 | |
| 69 if (quad_fcn) | |
| 70 { | |
| 3544 | 71 octave_value_list tmp = quad_fcn->do_multi_index_op (1, args); |
| 2928 | 72 |
| 73 if (error_state) | |
| 74 { | |
| 5775 | 75 quad_integration_error = 1; // FIXME |
| 2928 | 76 gripe_user_supplied_eval ("quad"); |
| 77 return retval; | |
| 78 } | |
| 79 | |
| 80 if (tmp.length () && tmp(0).is_defined ()) | |
| 81 { | |
| 4140 | 82 if (! warned_imaginary && tmp(0).is_complex_type ()) |
| 83 { | |
| 84 warning ("quad: ignoring imaginary part returned from user-supplied function"); | |
| 85 warned_imaginary = true; | |
| 86 } | |
| 87 | |
| 2928 | 88 retval = tmp(0).double_value (); |
| 89 | |
| 90 if (error_state) | |
| 91 { | |
| 5775 | 92 quad_integration_error = 1; // FIXME |
| 2928 | 93 gripe_user_supplied_eval ("quad"); |
| 94 } | |
| 95 } | |
| 96 else | |
| 97 { | |
| 5775 | 98 quad_integration_error = 1; // FIXME |
| 2928 | 99 gripe_user_supplied_eval ("quad"); |
| 100 } | |
| 101 } | |
| 102 | |
| 103 return retval; | |
| 104 } | |
| 105 | |
| 3323 | 106 #define QUAD_ABORT() \ |
| 107 do \ | |
| 108 { \ | |
| 4954 | 109 if (fcn_name.length()) \ |
| 110 clear_function (fcn_name); \ | |
| 3323 | 111 unwind_protect::run_frame ("Fquad"); \ |
| 112 return retval; \ | |
| 113 } \ | |
| 114 while (0) | |
| 115 | |
| 116 #define QUAD_ABORT1(msg) \ | |
| 117 do \ | |
| 118 { \ | |
| 3747 | 119 ::error ("quad: " msg); \ |
| 3323 | 120 QUAD_ABORT (); \ |
| 121 } \ | |
| 122 while (0) | |
| 123 | |
| 124 #define QUAD_ABORT2(fmt, arg) \ | |
| 125 do \ | |
| 126 { \ | |
| 3747 | 127 ::error ("quad: " fmt, arg); \ |
| 3323 | 128 QUAD_ABORT (); \ |
| 129 } \ | |
| 130 while (0) | |
| 131 | |
| 2928 | 132 DEFUN_DLD (quad, args, nargout, |
| 3368 | 133 "-*- texinfo -*-\n\ |
| 134 @deftypefn {Loadable Function} {[@var{v}, @var{ier}, @var{nfun}, @var{err}] =} quad (@var{f}, @var{a}, @var{b}, @var{tol}, @var{sing})\n\ | |
| 135 Integrate a nonlinear function of one variable using Quadpack.\n\ | |
| 4954 | 136 The first argument is the name of the function, the function handle or\n\ |
| 137 the inline function to call to compute the value of the integrand. It\n\ | |
| 138 must have the form\n\ | |
| 2928 | 139 \n\ |
| 3368 | 140 @example\n\ |
| 141 y = f (x)\n\ | |
| 142 @end example\n\ | |
| 2928 | 143 \n\ |
| 3368 | 144 @noindent\n\ |
| 145 where @var{y} and @var{x} are scalars.\n\ | |
| 2928 | 146 \n\ |
| 147 The second and third arguments are limits of integration. Either or\n\ | |
| 148 both may be infinite.\n\ | |
| 149 \n\ | |
| 3368 | 150 The optional argument @var{tol} is a vector that specifies the desired\n\ |
| 2928 | 151 accuracy of the result. The first element of the vector is the desired\n\ |
| 152 absolute tolerance, and the second element is the desired relative\n\ | |
| 3368 | 153 tolerance. To choose a relative test only, set the absolute\n\ |
| 154 tolerance to zero. To choose an absolute test only, set the relative\n\ | |
| 155 tolerance to zero. \n\ | |
| 2928 | 156 \n\ |
| 157 The optional argument @var{sing} is a vector of values at which the\n\ | |
| 3368 | 158 integrand is known to be singular.\n\ |
| 159 \n\ | |
| 160 The result of the integration is returned in @var{v} and @var{ier}\n\ | |
| 161 contains an integer error code (0 indicates a successful integration).\n\ | |
| 162 The value of @var{nfun} indicates how many function evaluations were\n\ | |
| 163 required, and @var{err} contains an estimate of the error in the\n\ | |
| 164 solution.\n\ | |
| 3964 | 165 \n\ |
| 166 You can use the function @code{quad_options} to set optional\n\ | |
| 167 parameters for @code{quad}.\n\ | |
| 6741 | 168 \n\ |
| 169 It should be noted that since @code{quad} is written in Fortran it\n\ | |
| 170 cannot be called recursively.\n\ | |
| 3368 | 171 @end deftypefn") |
| 2928 | 172 { |
| 173 octave_value_list retval; | |
| 174 | |
| 4954 | 175 std::string fcn_name; |
| 176 | |
| 4140 | 177 warned_imaginary = false; |
| 178 | |
| 3243 | 179 unwind_protect::begin_frame ("Fquad"); |
| 180 | |
| 181 unwind_protect_int (call_depth); | |
| 182 call_depth++; | |
| 183 | |
| 184 if (call_depth > 1) | |
| 3323 | 185 QUAD_ABORT1 ("invalid recursive call"); |
| 3243 | 186 |
| 2928 | 187 int nargin = args.length (); |
| 188 | |
| 3243 | 189 if (nargin > 2 && nargin < 6 && nargout < 5) |
| 190 { | |
| 4954 | 191 if (args(0).is_function_handle () || args(0).is_inline_function ()) |
| 192 quad_fcn = args(0).function_value (); | |
| 193 else | |
| 194 { | |
| 4962 | 195 fcn_name = unique_symbol_name ("__quad_fcn_"); |
| 4954 | 196 std::string fname = "function y = "; |
| 197 fname.append (fcn_name); | |
| 198 fname.append ("(x) y = "); | |
| 199 quad_fcn = extract_function (args(0), "quad", fcn_name, fname, | |
| 200 "; endfunction"); | |
| 201 } | |
| 202 | |
| 3243 | 203 if (! quad_fcn) |
| 3323 | 204 QUAD_ABORT (); |
| 3243 | 205 |
| 206 double a = args(1).double_value (); | |
| 207 | |
| 208 if (error_state) | |
| 3323 | 209 QUAD_ABORT1 ("expecting second argument to be a scalar"); |
| 3243 | 210 |
| 211 double b = args(2).double_value (); | |
| 212 | |
| 213 if (error_state) | |
| 3323 | 214 QUAD_ABORT1 ("expecting third argument to be a scalar"); |
| 3243 | 215 |
| 216 int indefinite = 0; | |
| 217 IndefQuad::IntegralType indef_type = IndefQuad::doubly_infinite; | |
| 218 double bound = 0.0; | |
| 219 if (xisinf (a) && xisinf (b)) | |
| 220 { | |
| 221 indefinite = 1; | |
| 222 indef_type = IndefQuad::doubly_infinite; | |
| 223 } | |
| 224 else if (xisinf (a)) | |
| 225 { | |
| 226 indefinite = 1; | |
| 227 bound = b; | |
| 228 indef_type = IndefQuad::neg_inf_to_bound; | |
| 229 } | |
| 230 else if (xisinf (b)) | |
| 231 { | |
| 232 indefinite = 1; | |
| 233 bound = a; | |
| 234 indef_type = IndefQuad::bound_to_inf; | |
| 235 } | |
| 236 | |
| 5275 | 237 octave_idx_type ier = 0; |
| 238 octave_idx_type nfun = 0; | |
| 3243 | 239 double abserr = 0.0; |
| 240 double val = 0.0; | |
| 3998 | 241 bool have_sing = false; |
| 3243 | 242 ColumnVector sing; |
| 3998 | 243 ColumnVector tol; |
| 244 | |
| 3243 | 245 switch (nargin) |
| 246 { | |
| 247 case 5: | |
| 248 if (indefinite) | |
| 3323 | 249 QUAD_ABORT1 ("singularities not allowed on infinite intervals"); |
| 3243 | 250 |
| 3998 | 251 have_sing = true; |
| 3243 | 252 |
| 3419 | 253 sing = ColumnVector (args(4).vector_value ()); |
| 3243 | 254 |
| 255 if (error_state) | |
| 3323 | 256 QUAD_ABORT1 ("expecting vector of singularities as fourth argument"); |
| 3243 | 257 |
| 258 case 4: | |
| 3419 | 259 tol = ColumnVector (args(3).vector_value ()); |
| 3243 | 260 |
| 261 if (error_state) | |
| 3323 | 262 QUAD_ABORT1 ("expecting vector of tolerances as fifth argument"); |
| 3243 | 263 |
| 264 switch (tol.capacity ()) | |
| 265 { | |
| 266 case 2: | |
| 3998 | 267 quad_opts.set_relative_tolerance (tol (1)); |
| 3243 | 268 |
| 269 case 1: | |
| 3998 | 270 quad_opts.set_absolute_tolerance (tol (0)); |
| 3243 | 271 break; |
| 272 | |
| 273 default: | |
| 3323 | 274 QUAD_ABORT1 ("expecting tol to contain no more than two values"); |
| 3243 | 275 } |
| 276 | |
| 277 case 3: | |
| 278 if (indefinite) | |
| 279 { | |
| 3998 | 280 IndefQuad iq (quad_user_function, bound, indef_type); |
| 4122 | 281 iq.set_options (quad_opts); |
| 3243 | 282 val = iq.integrate (ier, nfun, abserr); |
| 283 } | |
| 284 else | |
| 285 { | |
| 286 if (have_sing) | |
| 287 { | |
| 3998 | 288 DefQuad dq (quad_user_function, a, b, sing); |
| 4122 | 289 dq.set_options (quad_opts); |
| 3243 | 290 val = dq.integrate (ier, nfun, abserr); |
| 291 } | |
| 292 else | |
| 293 { | |
| 3998 | 294 DefQuad dq (quad_user_function, a, b); |
| 4122 | 295 dq.set_options (quad_opts); |
| 3243 | 296 val = dq.integrate (ier, nfun, abserr); |
| 297 } | |
| 298 } | |
| 299 break; | |
| 300 | |
| 301 default: | |
| 302 panic_impossible (); | |
| 303 break; | |
| 304 } | |
| 305 | |
| 306 retval(3) = abserr; | |
| 307 retval(2) = static_cast<double> (nfun); | |
| 308 retval(1) = static_cast<double> (ier); | |
| 309 retval(0) = val; | |
| 4954 | 310 |
| 311 if (fcn_name.length()) | |
| 312 clear_function (fcn_name); | |
| 3243 | 313 } |
| 314 else | |
| 5823 | 315 print_usage (); |
| 2928 | 316 |
| 3243 | 317 unwind_protect::run_frame ("Fquad"); |
| 2928 | 318 |
| 319 return retval; | |
| 320 } | |
| 321 | |
| 322 /* | |
|
7562
c827f5673321
move tests to individual source files
John W. Eaton <jwe@octave.org>
parents:
7017
diff
changeset
|
323 |
|
c827f5673321
move tests to individual source files
John W. Eaton <jwe@octave.org>
parents:
7017
diff
changeset
|
324 %!function y = f (x) |
|
c827f5673321
move tests to individual source files
John W. Eaton <jwe@octave.org>
parents:
7017
diff
changeset
|
325 %! y = x + 1; |
|
c827f5673321
move tests to individual source files
John W. Eaton <jwe@octave.org>
parents:
7017
diff
changeset
|
326 %!test |
|
c827f5673321
move tests to individual source files
John W. Eaton <jwe@octave.org>
parents:
7017
diff
changeset
|
327 %! [v, ier, nfun, err] = quad ("f", 0, 5); |
|
c827f5673321
move tests to individual source files
John W. Eaton <jwe@octave.org>
parents:
7017
diff
changeset
|
328 %! assert(ier == 0 && abs (v - 17.5) < sqrt (eps) && nfun > 0 && |
|
c827f5673321
move tests to individual source files
John W. Eaton <jwe@octave.org>
parents:
7017
diff
changeset
|
329 %! err < sqrt (eps)) |
|
c827f5673321
move tests to individual source files
John W. Eaton <jwe@octave.org>
parents:
7017
diff
changeset
|
330 |
|
c827f5673321
move tests to individual source files
John W. Eaton <jwe@octave.org>
parents:
7017
diff
changeset
|
331 %!function y = f (x) |
|
c827f5673321
move tests to individual source files
John W. Eaton <jwe@octave.org>
parents:
7017
diff
changeset
|
332 %! y = x .* sin (1 ./ x) .* sqrt (abs (1 - x)); |
|
c827f5673321
move tests to individual source files
John W. Eaton <jwe@octave.org>
parents:
7017
diff
changeset
|
333 %!test |
|
c827f5673321
move tests to individual source files
John W. Eaton <jwe@octave.org>
parents:
7017
diff
changeset
|
334 %! [v, ier, nfun, err] = quad ("f", 0.001, 3); |
|
c827f5673321
move tests to individual source files
John W. Eaton <jwe@octave.org>
parents:
7017
diff
changeset
|
335 %! assert((ier == 0 || ier == 1) && abs (v - 1.98194120273598) < sqrt (eps) && nfun > 0); |
|
c827f5673321
move tests to individual source files
John W. Eaton <jwe@octave.org>
parents:
7017
diff
changeset
|
336 |
|
c827f5673321
move tests to individual source files
John W. Eaton <jwe@octave.org>
parents:
7017
diff
changeset
|
337 %!error <Invalid call to quad.*> quad (); |
|
c827f5673321
move tests to individual source files
John W. Eaton <jwe@octave.org>
parents:
7017
diff
changeset
|
338 |
|
c827f5673321
move tests to individual source files
John W. Eaton <jwe@octave.org>
parents:
7017
diff
changeset
|
339 %!error <Invalid call to quad.*> quad ("f", 1, 2, 3, 4, 5); |
|
c827f5673321
move tests to individual source files
John W. Eaton <jwe@octave.org>
parents:
7017
diff
changeset
|
340 |
|
c827f5673321
move tests to individual source files
John W. Eaton <jwe@octave.org>
parents:
7017
diff
changeset
|
341 %!test |
|
c827f5673321
move tests to individual source files
John W. Eaton <jwe@octave.org>
parents:
7017
diff
changeset
|
342 %! quad_options ("absolute tolerance", eps); |
|
c827f5673321
move tests to individual source files
John W. Eaton <jwe@octave.org>
parents:
7017
diff
changeset
|
343 %! assert(quad_options ("absolute tolerance") == eps); |
|
c827f5673321
move tests to individual source files
John W. Eaton <jwe@octave.org>
parents:
7017
diff
changeset
|
344 |
|
c827f5673321
move tests to individual source files
John W. Eaton <jwe@octave.org>
parents:
7017
diff
changeset
|
345 %!error <Invalid call to quad_options.*> quad_options (1, 2, 3); |
|
c827f5673321
move tests to individual source files
John W. Eaton <jwe@octave.org>
parents:
7017
diff
changeset
|
346 |
|
c827f5673321
move tests to individual source files
John W. Eaton <jwe@octave.org>
parents:
7017
diff
changeset
|
347 */ |
|
c827f5673321
move tests to individual source files
John W. Eaton <jwe@octave.org>
parents:
7017
diff
changeset
|
348 |
|
c827f5673321
move tests to individual source files
John W. Eaton <jwe@octave.org>
parents:
7017
diff
changeset
|
349 /* |
| 2928 | 350 ;;; Local Variables: *** |
| 351 ;;; mode: C++ *** | |
| 352 ;;; End: *** | |
| 353 */ |