octave: scripts/general/quadl.m annotate

annotate scripts/general/quadl.m @ 20231:83792dd9bcc1

Use in-place operators in m-files where possible. * scripts/audio/@audioplayer/set.m, scripts/audio/@audiorecorder/set.m, scripts/audio/mu2lin.m, scripts/elfun/cosd.m, scripts/general/del2.m, scripts/general/profexplore.m, scripts/general/quadl.m, scripts/general/rat.m, scripts/general/rotdim.m, scripts/help/get_first_help_sentence.m, scripts/help/private/__strip_html_tags__.m, scripts/image/cubehelix.m, scripts/io/textread.m, scripts/linear-algebra/duplication_matrix.m, scripts/linear-algebra/housh.m, scripts/linear-algebra/krylov.m, scripts/linear-algebra/logm.m, scripts/linear-algebra/normest.m, scripts/linear-algebra/onenormest.m, scripts/optimization/fminsearch.m, scripts/optimization/lsqnonneg.m, scripts/optimization/qp.m, scripts/plot/appearance/annotation.m, scripts/plot/appearance/axis.m, scripts/plot/appearance/legend.m, scripts/plot/appearance/specular.m, scripts/plot/draw/colorbar.m, scripts/plot/draw/hist.m, scripts/plot/draw/plotmatrix.m, scripts/plot/draw/private/__stem__.m, scripts/plot/util/__actual_axis_position__.m, scripts/plot/util/__gnuplot_drawnow__.m, scripts/plot/util/findobj.m, scripts/plot/util/print.m, scripts/plot/util/private/__go_draw_axes__.m, scripts/plot/util/private/__print_parse_opts__.m, scripts/plot/util/rotate.m, scripts/polynomial/pchip.m, scripts/polynomial/polyaffine.m, scripts/polynomial/polyder.m, scripts/polynomial/private/__splinefit__.m, scripts/polynomial/residue.m, scripts/signal/arch_fit.m, scripts/signal/arch_rnd.m, scripts/signal/bartlett.m, scripts/signal/blackman.m, scripts/signal/freqz.m, scripts/signal/hamming.m, scripts/signal/hanning.m, scripts/signal/spectral_adf.m, scripts/signal/spectral_xdf.m, scripts/signal/stft.m, scripts/sparse/bicgstab.m, scripts/sparse/cgs.m, scripts/sparse/private/__sprand_impl__.m, scripts/sparse/qmr.m, scripts/sparse/sprandsym.m, scripts/sparse/svds.m, scripts/specfun/legendre.m, scripts/special-matrix/gallery.m, scripts/statistics/base/gls.m, scripts/statistics/models/logistic_regression.m, scripts/statistics/tests/kruskal_wallis_test.m, scripts/statistics/tests/manova.m, scripts/statistics/tests/wilcoxon_test.m, scripts/time/datevec.m: Use in-place operators in m-files where possible.

author	Rik <rik@octave.org>
date	Tue, 26 May 2015 21:07:42 -0700
parents	7503499a252b
children	655816377845

rev	line source
19697 4197fc428c7d maint: Update copyright notices for 2015. John W. Eaton <jwe@octave.org> parents: 19040 diff changeset	1 ## Copyright (C) 1998-2015 Walter Gautschi
5837 55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	2 ##
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	3 ## This file is part of Octave.
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	4 ##
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	5 ## Octave is free software; you can redistribute it and/or modify it
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	6 ## under the terms of the GNU General Public License as published by
7016 93c65f2a5668 [project @ 2007-10-12 06:40:56 by jwe] jwe parents: 5838 diff changeset	7 ## the Free Software Foundation; either version 3 of the License, or (at
93c65f2a5668 [project @ 2007-10-12 06:40:56 by jwe] jwe parents: 5838 diff changeset	8 ## your option) any later version.
5837 55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	9 ##
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	10 ## Octave is distributed in the hope that it will be useful, but
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	11 ## WITHOUT ANY WARRANTY; without even the implied warranty of
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	12 ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	13 ## General Public License for more details.
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	14 ##
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	15 ## You should have received a copy of the GNU General Public License
7016 93c65f2a5668 [project @ 2007-10-12 06:40:56 by jwe] jwe parents: 5838 diff changeset	16 ## along with Octave; see the file COPYING. If not, see
93c65f2a5668 [project @ 2007-10-12 06:40:56 by jwe] jwe parents: 5838 diff changeset	17 ## <http://www.gnu.org/licenses/>.
5837 55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	18
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	19 ## -- texinfo --
10793 be55736a0783 Grammarcheck the documentation from m-files. Rik <octave@nomad.inbox5.com> parents: 10791 diff changeset	20 ## @deftypefn {Function File} {@var{q} =} quadl (@var{f}, @var{a}, @var{b})
5837 55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	21 ## @deftypefnx {Function File} {@var{q} =} quadl (@var{f}, @var{a}, @var{b}, @var{tol})
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	22 ## @deftypefnx {Function File} {@var{q} =} quadl (@var{f}, @var{a}, @var{b}, @var{tol}, @var{trace})
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	23 ## @deftypefnx {Function File} {@var{q} =} quadl (@var{f}, @var{a}, @var{b}, @var{tol}, @var{trace}, @var{p1}, @var{p2}, @dots{})
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	24 ##
20158 7503499a252b doc: Update docstrings to have one sentence summary as first line. Rik <rik@octave.org> parents: 19697 diff changeset	25 ## Numerically evaluate the integral of @var{f} from @var{a} to @var{b} using
7503499a252b doc: Update docstrings to have one sentence summary as first line. Rik <rik@octave.org> parents: 19697 diff changeset	26 ## an adaptive Lobatto rule.
7503499a252b doc: Update docstrings to have one sentence summary as first line. Rik <rik@octave.org> parents: 19697 diff changeset	27 ##
7503499a252b doc: Update docstrings to have one sentence summary as first line. Rik <rik@octave.org> parents: 19697 diff changeset	28 ## @var{f} is a function handle, inline function, or string containing the name
7503499a252b doc: Update docstrings to have one sentence summary as first line. Rik <rik@octave.org> parents: 19697 diff changeset	29 ## of the function to evaluate. The function @var{f} must be vectorized and
7503499a252b doc: Update docstrings to have one sentence summary as first line. Rik <rik@octave.org> parents: 19697 diff changeset	30 ## return a vector of output values when given a vector of input values.
12612 16cca721117b doc: Update all documentation for chapter on Numerical Integration Rik <octave@nomad.inbox5.com> parents: 12575 diff changeset	31 ##
16cca721117b doc: Update all documentation for chapter on Numerical Integration Rik <octave@nomad.inbox5.com> parents: 12575 diff changeset	32 ## @var{a} and @var{b} are the lower and upper limits of integration. Both
16cca721117b doc: Update all documentation for chapter on Numerical Integration Rik <octave@nomad.inbox5.com> parents: 12575 diff changeset	33 ## limits must be finite.
5837 55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	34 ##
12612 16cca721117b doc: Update all documentation for chapter on Numerical Integration Rik <octave@nomad.inbox5.com> parents: 12575 diff changeset	35 ## The optional argument @var{tol} defines the relative tolerance with which
16cca721117b doc: Update all documentation for chapter on Numerical Integration Rik <octave@nomad.inbox5.com> parents: 12575 diff changeset	36 ## to perform the integration. The default value is @code{eps}.
5837 55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	37 ##
12612 16cca721117b doc: Update all documentation for chapter on Numerical Integration Rik <octave@nomad.inbox5.com> parents: 12575 diff changeset	38 ## The algorithm used by @code{quadl} involves recursively subdividing the
20158 7503499a252b doc: Update docstrings to have one sentence summary as first line. Rik <rik@octave.org> parents: 19697 diff changeset	39 ## integration interval. If @var{trace} is defined then for each subinterval
7503499a252b doc: Update docstrings to have one sentence summary as first line. Rik <rik@octave.org> parents: 19697 diff changeset	40 ## display: (1) the left end of the subinterval, (2) the length of the
7503499a252b doc: Update docstrings to have one sentence summary as first line. Rik <rik@octave.org> parents: 19697 diff changeset	41 ## subinterval, (3) the approximation of the integral over the subinterval.
12612 16cca721117b doc: Update all documentation for chapter on Numerical Integration Rik <octave@nomad.inbox5.com> parents: 12575 diff changeset	42 ##
16cca721117b doc: Update all documentation for chapter on Numerical Integration Rik <octave@nomad.inbox5.com> parents: 12575 diff changeset	43 ## Additional arguments @var{p1}, etc., are passed directly to the function
16cca721117b doc: Update all documentation for chapter on Numerical Integration Rik <octave@nomad.inbox5.com> parents: 12575 diff changeset	44 ## @var{f}. To use default values for @var{tol} and @var{trace}, one may pass
16cca721117b doc: Update all documentation for chapter on Numerical Integration Rik <octave@nomad.inbox5.com> parents: 12575 diff changeset	45 ## empty matrices ([]).
5837 55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	46 ##
19040 0850b5212619 doc: Add @nospell macro around proper names in documentation. Rik <rik@octave.org> parents: 17744 diff changeset	47 ## Reference: @nospell{W. Gander and W. Gautschi}, @cite{Adaptive Quadrature -
10791 3140cb7a05a1 Add spellchecker scripts for Octave and run spellcheck of documentation Rik <octave@nomad.inbox5.com> parents: 10549 diff changeset	48 ## Revisited}, BIT Vol. 40, No. 1, March 2000, pp. 84--101.
5837 55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	49 ## @url{http://www.inf.ethz.ch/personal/gander/}
12575 d0b799dafede Grammarcheck files for 3.4.1 release. Rik <octave@nomad.inbox5.com> parents: 11587 diff changeset	50 ## @seealso{quad, quadv, quadgk, quadcc, trapz, dblquad, triplequad}
5837 55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	51 ## @end deftypefn
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	52
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	53 ## Author: Walter Gautschi
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	54 ## Date: 08/03/98
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	55 ## Reference: Gander, Computermathematik, Birkhaeuser, 1992.
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	56
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	57 ## 2003-08-05 Shai Ayal
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	58 ## * permission from author to release as GPL
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	59 ## 2004-02-10 Paul Kienzle
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	60 ## * renamed to quadl for compatibility
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	61 ## * replace global variable terminate2 with local function need_warning
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	62 ## * add paper ref to docs
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	63
13008 85dac13a911b quadl.m: Fix integration with large error tolerances (Bug #33792) Rik <octave@nomad.inbox5.com> parents: 12802 diff changeset	64 function q = quadl (f, a, b, tol = [], trace = false, varargin)
85dac13a911b quadl.m: Fix integration with large error tolerances (Bug #33792) Rik <octave@nomad.inbox5.com> parents: 12802 diff changeset	65
85dac13a911b quadl.m: Fix integration with large error tolerances (Bug #33792) Rik <octave@nomad.inbox5.com> parents: 12802 diff changeset	66 if (nargin < 3)
85dac13a911b quadl.m: Fix integration with large error tolerances (Bug #33792) Rik <octave@nomad.inbox5.com> parents: 12802 diff changeset	67 print_usage ();
5837 55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	68 endif
13008 85dac13a911b quadl.m: Fix integration with large error tolerances (Bug #33792) Rik <octave@nomad.inbox5.com> parents: 12802 diff changeset	69
7795 df9519e9990c Handle single precision eps values David Bateman <dbateman@free.fr> parents: 7017 diff changeset	70 if (isa (a, "single") \|\| isa (b, "single"))
df9519e9990c Handle single precision eps values David Bateman <dbateman@free.fr> parents: 7017 diff changeset	71 myeps = eps ("single");
df9519e9990c Handle single precision eps values David Bateman <dbateman@free.fr> parents: 7017 diff changeset	72 else
df9519e9990c Handle single precision eps values David Bateman <dbateman@free.fr> parents: 7017 diff changeset	73 myeps = eps;
df9519e9990c Handle single precision eps values David Bateman <dbateman@free.fr> parents: 7017 diff changeset	74 endif
5838 376e02b2ce70 [project @ 2006-06-01 20:23:53 by jwe] jwe parents: 5837 diff changeset	75 if (isempty (tol))
11587 c792872f8942 all script files: untabify and strip trailing whitespace John W. Eaton <jwe@octave.org> parents: 11562 diff changeset	76 tol = myeps;
5837 55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	77 endif
5838 376e02b2ce70 [project @ 2006-06-01 20:23:53 by jwe] jwe parents: 5837 diff changeset	78 if (isempty (trace))
13008 85dac13a911b quadl.m: Fix integration with large error tolerances (Bug #33792) Rik <octave@nomad.inbox5.com> parents: 12802 diff changeset	79 trace = false;
5837 55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	80 endif
7795 df9519e9990c Handle single precision eps values David Bateman <dbateman@free.fr> parents: 7017 diff changeset	81 if (tol < myeps)
df9519e9990c Handle single precision eps values David Bateman <dbateman@free.fr> parents: 7017 diff changeset	82 tol = myeps;
5837 55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	83 endif
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	84
13008 85dac13a911b quadl.m: Fix integration with large error tolerances (Bug #33792) Rik <octave@nomad.inbox5.com> parents: 12802 diff changeset	85 ## Track whether recursion has occurred
85dac13a911b quadl.m: Fix integration with large error tolerances (Bug #33792) Rik <octave@nomad.inbox5.com> parents: 12802 diff changeset	86 global __quadl_recurse_done__;
85dac13a911b quadl.m: Fix integration with large error tolerances (Bug #33792) Rik <octave@nomad.inbox5.com> parents: 12802 diff changeset	87 __quadl_recurse_done__ = false;
85dac13a911b quadl.m: Fix integration with large error tolerances (Bug #33792) Rik <octave@nomad.inbox5.com> parents: 12802 diff changeset	88 ## Track whether warning about machine precision has been issued
85dac13a911b quadl.m: Fix integration with large error tolerances (Bug #33792) Rik <octave@nomad.inbox5.com> parents: 12802 diff changeset	89 global __quadl_need_warning__;
85dac13a911b quadl.m: Fix integration with large error tolerances (Bug #33792) Rik <octave@nomad.inbox5.com> parents: 12802 diff changeset	90 __quadl_need_warning__ = true;
85dac13a911b quadl.m: Fix integration with large error tolerances (Bug #33792) Rik <octave@nomad.inbox5.com> parents: 12802 diff changeset	91
11587 c792872f8942 all script files: untabify and strip trailing whitespace John W. Eaton <jwe@octave.org> parents: 11562 diff changeset	92 m = (a+b)/2;
5837 55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	93 h = (b-a)/2;
13008 85dac13a911b quadl.m: Fix integration with large error tolerances (Bug #33792) Rik <octave@nomad.inbox5.com> parents: 12802 diff changeset	94 alpha = sqrt (2/3);
85dac13a911b quadl.m: Fix integration with large error tolerances (Bug #33792) Rik <octave@nomad.inbox5.com> parents: 12802 diff changeset	95 beta = 1/sqrt (5);
5838 376e02b2ce70 [project @ 2006-06-01 20:23:53 by jwe] jwe parents: 5837 diff changeset	96
11587 c792872f8942 all script files: untabify and strip trailing whitespace John W. Eaton <jwe@octave.org> parents: 11562 diff changeset	97 x1 = .942882415695480;
5837 55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	98 x2 = .641853342345781;
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	99 x3 = .236383199662150;
5838 376e02b2ce70 [project @ 2006-06-01 20:23:53 by jwe] jwe parents: 5837 diff changeset	100
376e02b2ce70 [project @ 2006-06-01 20:23:53 by jwe] jwe parents: 5837 diff changeset	101 x = [a, m-x1h, m-alphah, m-x2h, m-betah, m-x3h, m, m+x3h, ...
376e02b2ce70 [project @ 2006-06-01 20:23:53 by jwe] jwe parents: 5837 diff changeset	102 m+betah, m+x2h, m+alphah, m+x1h, b];
376e02b2ce70 [project @ 2006-06-01 20:23:53 by jwe] jwe parents: 5837 diff changeset	103
376e02b2ce70 [project @ 2006-06-01 20:23:53 by jwe] jwe parents: 5837 diff changeset	104 y = feval (f, x, varargin{:});
376e02b2ce70 [project @ 2006-06-01 20:23:53 by jwe] jwe parents: 5837 diff changeset	105
11587 c792872f8942 all script files: untabify and strip trailing whitespace John W. Eaton <jwe@octave.org> parents: 11562 diff changeset	106 fa = y(1);
5837 55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	107 fb = y(13);
5838 376e02b2ce70 [project @ 2006-06-01 20:23:53 by jwe] jwe parents: 5837 diff changeset	108
376e02b2ce70 [project @ 2006-06-01 20:23:53 by jwe] jwe parents: 5837 diff changeset	109 i2 = (h/6)(y(1) + y(13) + 5(y(5)+y(9)));
376e02b2ce70 [project @ 2006-06-01 20:23:53 by jwe] jwe parents: 5837 diff changeset	110
13008 85dac13a911b quadl.m: Fix integration with large error tolerances (Bug #33792) Rik <octave@nomad.inbox5.com> parents: 12802 diff changeset	111 i1 = (h/1470)( 77(y(1)+y(13))
10549 95c3e38098bf Untabify .m scripts Rik <code@nomad.inbox5.com> parents: 9245 diff changeset	112 + 432*(y(3)+y(11))
95c3e38098bf Untabify .m scripts Rik <code@nomad.inbox5.com> parents: 9245 diff changeset	113 + 625*(y(5)+y(9))
95c3e38098bf Untabify .m scripts Rik <code@nomad.inbox5.com> parents: 9245 diff changeset	114 + 672*y(7));
5838 376e02b2ce70 [project @ 2006-06-01 20:23:53 by jwe] jwe parents: 5837 diff changeset	115
13008 85dac13a911b quadl.m: Fix integration with large error tolerances (Bug #33792) Rik <octave@nomad.inbox5.com> parents: 12802 diff changeset	116 is = h( .0158271919734802(y(1)+y(13))
10549 95c3e38098bf Untabify .m scripts Rik <code@nomad.inbox5.com> parents: 9245 diff changeset	117 +.0942738402188500*(y(2)+y(12))
95c3e38098bf Untabify .m scripts Rik <code@nomad.inbox5.com> parents: 9245 diff changeset	118 + .155071987336585*(y(3)+y(11))
95c3e38098bf Untabify .m scripts Rik <code@nomad.inbox5.com> parents: 9245 diff changeset	119 + .188821573960182*(y(4)+y(10))
95c3e38098bf Untabify .m scripts Rik <code@nomad.inbox5.com> parents: 9245 diff changeset	120 + .199773405226859*(y(5)+y(9))
95c3e38098bf Untabify .m scripts Rik <code@nomad.inbox5.com> parents: 9245 diff changeset	121 + .224926465333340*(y(6)+y(8))
95c3e38098bf Untabify .m scripts Rik <code@nomad.inbox5.com> parents: 9245 diff changeset	122 + .242611071901408*y(7));
5838 376e02b2ce70 [project @ 2006-06-01 20:23:53 by jwe] jwe parents: 5837 diff changeset	123
13008 85dac13a911b quadl.m: Fix integration with large error tolerances (Bug #33792) Rik <octave@nomad.inbox5.com> parents: 12802 diff changeset	124 s = sign (is);
5838 376e02b2ce70 [project @ 2006-06-01 20:23:53 by jwe] jwe parents: 5837 diff changeset	125 if (s == 0)
376e02b2ce70 [project @ 2006-06-01 20:23:53 by jwe] jwe parents: 5837 diff changeset	126 s = 1;
5837 55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	127 endif
13008 85dac13a911b quadl.m: Fix integration with large error tolerances (Bug #33792) Rik <octave@nomad.inbox5.com> parents: 12802 diff changeset	128 erri1 = abs (i1-is);
85dac13a911b quadl.m: Fix integration with large error tolerances (Bug #33792) Rik <octave@nomad.inbox5.com> parents: 12802 diff changeset	129 erri2 = abs (i2-is);
5837 55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	130 if (erri2 != 0)
11587 c792872f8942 all script files: untabify and strip trailing whitespace John W. Eaton <jwe@octave.org> parents: 11562 diff changeset	131 R = erri1/erri2;
13008 85dac13a911b quadl.m: Fix integration with large error tolerances (Bug #33792) Rik <octave@nomad.inbox5.com> parents: 12802 diff changeset	132 else
85dac13a911b quadl.m: Fix integration with large error tolerances (Bug #33792) Rik <octave@nomad.inbox5.com> parents: 12802 diff changeset	133 R = 1;
5837 55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	134 endif
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	135 if (R > 0 && R < 1)
20231 83792dd9bcc1 Use in-place operators in m-files where possible. Rik <rik@octave.org> parents: 20158 diff changeset	136 tol /= R;
5837 55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	137 endif
14868 5d3a684236b0 maint: Use Octave coding conventions for cuddling parentheses in scripts directory Rik <octave@nomad.inbox5.com> parents: 14138 diff changeset	138 is = s * abs (is) * tol/myeps;
5837 55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	139 if (is == 0)
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	140 is = b-a;
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	141 endif
13008 85dac13a911b quadl.m: Fix integration with large error tolerances (Bug #33792) Rik <octave@nomad.inbox5.com> parents: 12802 diff changeset	142
12612 16cca721117b doc: Update all documentation for chapter on Numerical Integration Rik <octave@nomad.inbox5.com> parents: 12575 diff changeset	143 q = adaptlobstp (f, a, b, fa, fb, is, trace, varargin{:});
13008 85dac13a911b quadl.m: Fix integration with large error tolerances (Bug #33792) Rik <octave@nomad.inbox5.com> parents: 12802 diff changeset	144
5837 55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	145 endfunction
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	146
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	147 ## ADAPTLOBSTP Recursive function used by QUADL.
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	148 ##
5838 376e02b2ce70 [project @ 2006-06-01 20:23:53 by jwe] jwe parents: 5837 diff changeset	149 ## Q = ADAPTLOBSTP('F', A, B, FA, FB, IS, TRACE) tries to
5837 55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	150 ## approximate the integral of F(X) from A to B to
13008 85dac13a911b quadl.m: Fix integration with large error tolerances (Bug #33792) Rik <octave@nomad.inbox5.com> parents: 12802 diff changeset	151 ## an appropriate relative error. The argument 'F' is
5837 55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	152 ## a string containing the name of f. The remaining
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	153 ## arguments are generated by ADAPTLOB or by recursion.
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	154 ##
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	155 ## Walter Gautschi, 08/03/98
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	156
12612 16cca721117b doc: Update all documentation for chapter on Numerical Integration Rik <octave@nomad.inbox5.com> parents: 12575 diff changeset	157 function q = adaptlobstp (f, a, b, fa, fb, is, trace, varargin)
13008 85dac13a911b quadl.m: Fix integration with large error tolerances (Bug #33792) Rik <octave@nomad.inbox5.com> parents: 12802 diff changeset	158 global __quadl_recurse_done__;
85dac13a911b quadl.m: Fix integration with large error tolerances (Bug #33792) Rik <octave@nomad.inbox5.com> parents: 12802 diff changeset	159 global __quadl_need_warning__;
85dac13a911b quadl.m: Fix integration with large error tolerances (Bug #33792) Rik <octave@nomad.inbox5.com> parents: 12802 diff changeset	160
11587 c792872f8942 all script files: untabify and strip trailing whitespace John W. Eaton <jwe@octave.org> parents: 11562 diff changeset	161 h = (b-a)/2;
5837 55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	162 m = (a+b)/2;
13008 85dac13a911b quadl.m: Fix integration with large error tolerances (Bug #33792) Rik <octave@nomad.inbox5.com> parents: 12802 diff changeset	163 alpha = sqrt (2/3);
14868 5d3a684236b0 maint: Use Octave coding conventions for cuddling parentheses in scripts directory Rik <octave@nomad.inbox5.com> parents: 14138 diff changeset	164 beta = 1 / sqrt (5);
11587 c792872f8942 all script files: untabify and strip trailing whitespace John W. Eaton <jwe@octave.org> parents: 11562 diff changeset	165 mll = m-alpha*h;
13008 85dac13a911b quadl.m: Fix integration with large error tolerances (Bug #33792) Rik <octave@nomad.inbox5.com> parents: 12802 diff changeset	166 ml = m-beta*h;
85dac13a911b quadl.m: Fix integration with large error tolerances (Bug #33792) Rik <octave@nomad.inbox5.com> parents: 12802 diff changeset	167 mr = m+beta*h;
5837 55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	168 mrr = m+alpha*h;
5838 376e02b2ce70 [project @ 2006-06-01 20:23:53 by jwe] jwe parents: 5837 diff changeset	169 x = [mll, ml, m, mr, mrr];
13008 85dac13a911b quadl.m: Fix integration with large error tolerances (Bug #33792) Rik <octave@nomad.inbox5.com> parents: 12802 diff changeset	170 y = feval (f, x, varargin{:});
11587 c792872f8942 all script files: untabify and strip trailing whitespace John W. Eaton <jwe@octave.org> parents: 11562 diff changeset	171 fmll = y(1);
13008 85dac13a911b quadl.m: Fix integration with large error tolerances (Bug #33792) Rik <octave@nomad.inbox5.com> parents: 12802 diff changeset	172 fml = y(2);
85dac13a911b quadl.m: Fix integration with large error tolerances (Bug #33792) Rik <octave@nomad.inbox5.com> parents: 12802 diff changeset	173 fm = y(3);
85dac13a911b quadl.m: Fix integration with large error tolerances (Bug #33792) Rik <octave@nomad.inbox5.com> parents: 12802 diff changeset	174 fmr = y(4);
5837 55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	175 fmrr = y(5);
5838 376e02b2ce70 [project @ 2006-06-01 20:23:53 by jwe] jwe parents: 5837 diff changeset	176 i2 = (h/6)(fa + fb + 5(fml+fmr));
376e02b2ce70 [project @ 2006-06-01 20:23:53 by jwe] jwe parents: 5837 diff changeset	177 i1 = (h/1470)(77(fa+fb) + 432(fmll+fmrr) + 625(fml+fmr) + 672*fm);
13008 85dac13a911b quadl.m: Fix integration with large error tolerances (Bug #33792) Rik <octave@nomad.inbox5.com> parents: 12802 diff changeset	178 if ((is+(i1-i2) == is \|\| mll <= a \|\| b <= mrr) && __quadl_recurse_done__)
85dac13a911b quadl.m: Fix integration with large error tolerances (Bug #33792) Rik <octave@nomad.inbox5.com> parents: 12802 diff changeset	179 if ((m <= a \|\| b <= m) && __quadl_need_warning__)
5838 376e02b2ce70 [project @ 2006-06-01 20:23:53 by jwe] jwe parents: 5837 diff changeset	180 warning ("quadl: interval contains no more machine number");
376e02b2ce70 [project @ 2006-06-01 20:23:53 by jwe] jwe parents: 5837 diff changeset	181 warning ("quadl: required tolerance may not be met");
13008 85dac13a911b quadl.m: Fix integration with large error tolerances (Bug #33792) Rik <octave@nomad.inbox5.com> parents: 12802 diff changeset	182 __quadl_need_warning__ = false;
5837 55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	183 endif
12612 16cca721117b doc: Update all documentation for chapter on Numerical Integration Rik <octave@nomad.inbox5.com> parents: 12575 diff changeset	184 q = i1;
5837 55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	185 if (trace)
12612 16cca721117b doc: Update all documentation for chapter on Numerical Integration Rik <octave@nomad.inbox5.com> parents: 12575 diff changeset	186 disp ([a, b-a, q]);
5837 55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	187 endif
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	188 else
13008 85dac13a911b quadl.m: Fix integration with large error tolerances (Bug #33792) Rik <octave@nomad.inbox5.com> parents: 12802 diff changeset	189 __quadl_recurse_done__ = true;
85dac13a911b quadl.m: Fix integration with large error tolerances (Bug #33792) Rik <octave@nomad.inbox5.com> parents: 12802 diff changeset	190 q = ( adaptlobstp (f, a , mll, fa , fmll, is, trace, varargin{:})
85dac13a911b quadl.m: Fix integration with large error tolerances (Bug #33792) Rik <octave@nomad.inbox5.com> parents: 12802 diff changeset	191 + adaptlobstp (f, mll, ml , fmll, fml , is, trace, varargin{:})
85dac13a911b quadl.m: Fix integration with large error tolerances (Bug #33792) Rik <octave@nomad.inbox5.com> parents: 12802 diff changeset	192 + adaptlobstp (f, ml , m , fml , fm , is, trace, varargin{:})
85dac13a911b quadl.m: Fix integration with large error tolerances (Bug #33792) Rik <octave@nomad.inbox5.com> parents: 12802 diff changeset	193 + adaptlobstp (f, m , mr , fm , fmr , is, trace, varargin{:})
85dac13a911b quadl.m: Fix integration with large error tolerances (Bug #33792) Rik <octave@nomad.inbox5.com> parents: 12802 diff changeset	194 + adaptlobstp (f, mr , mrr, fmr , fmrr, is, trace, varargin{:})
85dac13a911b quadl.m: Fix integration with large error tolerances (Bug #33792) Rik <octave@nomad.inbox5.com> parents: 12802 diff changeset	195 + adaptlobstp (f, mrr, b , fmrr, fb , is, trace, varargin{:}));
5837 55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	196 endif
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	197 endfunction
12802 412882f498b4 codesprint: Wrote 5 tests for quadl.m David Wells <drwells@vt.edu> parents: 12642 diff changeset	198
412882f498b4 codesprint: Wrote 5 tests for quadl.m David Wells <drwells@vt.edu> parents: 12642 diff changeset	199
412882f498b4 codesprint: Wrote 5 tests for quadl.m David Wells <drwells@vt.edu> parents: 12642 diff changeset	200 ## basic functionality
13008 85dac13a911b quadl.m: Fix integration with large error tolerances (Bug #33792) Rik <octave@nomad.inbox5.com> parents: 12802 diff changeset	201 %!assert (quadl (@(x) sin (x), 0, pi, [], []), 2, -3e-16)
12802 412882f498b4 codesprint: Wrote 5 tests for quadl.m David Wells <drwells@vt.edu> parents: 12642 diff changeset	202
412882f498b4 codesprint: Wrote 5 tests for quadl.m David Wells <drwells@vt.edu> parents: 12642 diff changeset	203 ## the values here are very high so it may be unavoidable that this fails
13008 85dac13a911b quadl.m: Fix integration with large error tolerances (Bug #33792) Rik <octave@nomad.inbox5.com> parents: 12802 diff changeset	204 %!assert (quadl (@(x) sin (3x).cosh (x).*sinh (x),10,15),
14068 6eeb9e8e63cf quadl.m: Relax tolerance on %!test to pass on all platforms (Bug #33863) Rik <octave@nomad.inbox5.com> parents: 13141 diff changeset	205 %! 2.588424538641647e+10, -1.1e-14)
12802 412882f498b4 codesprint: Wrote 5 tests for quadl.m David Wells <drwells@vt.edu> parents: 12642 diff changeset	206
412882f498b4 codesprint: Wrote 5 tests for quadl.m David Wells <drwells@vt.edu> parents: 12642 diff changeset	207 ## extra parameters
412882f498b4 codesprint: Wrote 5 tests for quadl.m David Wells <drwells@vt.edu> parents: 12642 diff changeset	208 %!assert (quadl (@(x,a,b) sin (a + b*x), 0, 1, [], [], 2, 3),
14068 6eeb9e8e63cf quadl.m: Relax tolerance on %!test to pass on all platforms (Bug #33863) Rik <octave@nomad.inbox5.com> parents: 13141 diff changeset	209 %! cos(2)/3 - cos(5)/3, -3e-16)
12802 412882f498b4 codesprint: Wrote 5 tests for quadl.m David Wells <drwells@vt.edu> parents: 12642 diff changeset	210
13141 e81ddf9cacd5 maint: untabify and remove trailing whitespace from source files John W. Eaton <jwe@octave.org> parents: 13008 diff changeset	211 ## test different tolerances.
13008 85dac13a911b quadl.m: Fix integration with large error tolerances (Bug #33792) Rik <octave@nomad.inbox5.com> parents: 12802 diff changeset	212 %!assert (quadl (@(x) sin (2 + 3*x).^2, 0, 10, 0.3, []),
12802 412882f498b4 codesprint: Wrote 5 tests for quadl.m David Wells <drwells@vt.edu> parents: 12642 diff changeset	213 %! (60 + sin(4) - sin(64))/12, -0.3)
13008 85dac13a911b quadl.m: Fix integration with large error tolerances (Bug #33792) Rik <octave@nomad.inbox5.com> parents: 12802 diff changeset	214 %!assert (quadl (@(x) sin (2 + 3*x).^2, 0, 10, 0.1, []),
85dac13a911b quadl.m: Fix integration with large error tolerances (Bug #33792) Rik <octave@nomad.inbox5.com> parents: 12802 diff changeset	215 %! (60 + sin(4) - sin(64))/12, -0.1)
12802 412882f498b4 codesprint: Wrote 5 tests for quadl.m David Wells <drwells@vt.edu> parents: 12642 diff changeset	216

Mercurial > octave

annotate scripts/general/quadl.m @ 20231:83792dd9bcc1