octave-nkf: scripts/general/quadl.m annotate

annotate scripts/general/quadl.m @ 19898:4197fc428c7d

maint: Update copyright notices for 2015.

author	John W. Eaton <jwe@octave.org>
date	Wed, 11 Feb 2015 14:19:08 -0500
parents	0850b5212619
children	7503499a252b

rev	line source
19898 4197fc428c7d maint: Update copyright notices for 2015. John W. Eaton <jwe@octave.org> parents: 19232 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 ##
12612 16cca721117b doc: Update all documentation for chapter on Numerical Integration Rik <octave@nomad.inbox5.com> parents: 12575 diff changeset	25 ## Numerically evaluate the integral of @var{f} from @var{a} to @var{b}
16cca721117b doc: Update all documentation for chapter on Numerical Integration Rik <octave@nomad.inbox5.com> parents: 12575 diff changeset	26 ## using an adaptive Lobatto rule.
16cca721117b doc: Update all documentation for chapter on Numerical Integration Rik <octave@nomad.inbox5.com> parents: 12575 diff changeset	27 ## @var{f} is a function handle, inline function, or string
16cca721117b doc: Update all documentation for chapter on Numerical Integration Rik <octave@nomad.inbox5.com> parents: 12575 diff changeset	28 ## containing the name of the function to evaluate.
16cca721117b doc: Update all documentation for chapter on Numerical Integration Rik <octave@nomad.inbox5.com> parents: 12575 diff changeset	29 ## The function @var{f} must be vectorized and return a vector of output values
16cca721117b doc: Update all documentation for chapter on Numerical Integration Rik <octave@nomad.inbox5.com> parents: 12575 diff changeset	30 ## if given a vector of input values.
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
16cca721117b doc: Update all documentation for chapter on Numerical Integration Rik <octave@nomad.inbox5.com> parents: 12575 diff changeset	39 ## integration interval.
16cca721117b doc: Update all documentation for chapter on Numerical Integration Rik <octave@nomad.inbox5.com> parents: 12575 diff changeset	40 ## If @var{trace} is defined then for each subinterval display: (1) the left
16cca721117b doc: Update all documentation for chapter on Numerical Integration Rik <octave@nomad.inbox5.com> parents: 12575 diff changeset	41 ## end of the subinterval, (2) the length of the subinterval, (3) the
12642 f96b9b9f141b doc: Periodic grammarcheck and spellcheck of documentation. Rik <octave@nomad.inbox5.com> parents: 12612 diff changeset	42 ## 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	43 ##
16cca721117b doc: Update all documentation for chapter on Numerical Integration Rik <octave@nomad.inbox5.com> parents: 12575 diff changeset	44 ## 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	45 ## @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	46 ## empty matrices ([]).
5837 55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	47 ##
19232 0850b5212619 doc: Add @nospell macro around proper names in documentation. Rik <rik@octave.org> parents: 17744 diff changeset	48 ## 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	49 ## 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	50 ## @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	51 ## @seealso{quad, quadv, quadgk, quadcc, trapz, dblquad, triplequad}
5837 55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	52 ## @end deftypefn
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	53
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	54 ## Author: Walter Gautschi
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	55 ## Date: 08/03/98
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	56 ## Reference: Gander, Computermathematik, Birkhaeuser, 1992.
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	57
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	58 ## 2003-08-05 Shai Ayal
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	59 ## * permission from author to release as GPL
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	60 ## 2004-02-10 Paul Kienzle
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	61 ## * renamed to quadl for compatibility
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	62 ## * replace global variable terminate2 with local function need_warning
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	63 ## * add paper ref to docs
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	64
13008 85dac13a911b quadl.m: Fix integration with large error tolerances (Bug #33792) Rik <octave@nomad.inbox5.com> parents: 12802 diff changeset	65 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	66
85dac13a911b quadl.m: Fix integration with large error tolerances (Bug #33792) Rik <octave@nomad.inbox5.com> parents: 12802 diff changeset	67 if (nargin < 3)
85dac13a911b quadl.m: Fix integration with large error tolerances (Bug #33792) Rik <octave@nomad.inbox5.com> parents: 12802 diff changeset	68 print_usage ();
5837 55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	69 endif
13008 85dac13a911b quadl.m: Fix integration with large error tolerances (Bug #33792) Rik <octave@nomad.inbox5.com> parents: 12802 diff changeset	70
7795 df9519e9990c Handle single precision eps values David Bateman <dbateman@free.fr> parents: 7017 diff changeset	71 if (isa (a, "single") \|\| isa (b, "single"))
df9519e9990c Handle single precision eps values David Bateman <dbateman@free.fr> parents: 7017 diff changeset	72 myeps = eps ("single");
df9519e9990c Handle single precision eps values David Bateman <dbateman@free.fr> parents: 7017 diff changeset	73 else
df9519e9990c Handle single precision eps values David Bateman <dbateman@free.fr> parents: 7017 diff changeset	74 myeps = eps;
df9519e9990c Handle single precision eps values David Bateman <dbateman@free.fr> parents: 7017 diff changeset	75 endif
5838 376e02b2ce70 [project @ 2006-06-01 20:23:53 by jwe] jwe parents: 5837 diff changeset	76 if (isempty (tol))
11587 c792872f8942 all script files: untabify and strip trailing whitespace John W. Eaton <jwe@octave.org> parents: 11562 diff changeset	77 tol = myeps;
5837 55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	78 endif
5838 376e02b2ce70 [project @ 2006-06-01 20:23:53 by jwe] jwe parents: 5837 diff changeset	79 if (isempty (trace))
13008 85dac13a911b quadl.m: Fix integration with large error tolerances (Bug #33792) Rik <octave@nomad.inbox5.com> parents: 12802 diff changeset	80 trace = false;
5837 55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	81 endif
7795 df9519e9990c Handle single precision eps values David Bateman <dbateman@free.fr> parents: 7017 diff changeset	82 if (tol < myeps)
df9519e9990c Handle single precision eps values David Bateman <dbateman@free.fr> parents: 7017 diff changeset	83 tol = myeps;
5837 55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	84 endif
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	85
13008 85dac13a911b quadl.m: Fix integration with large error tolerances (Bug #33792) Rik <octave@nomad.inbox5.com> parents: 12802 diff changeset	86 ## 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	87 global __quadl_recurse_done__;
85dac13a911b quadl.m: Fix integration with large error tolerances (Bug #33792) Rik <octave@nomad.inbox5.com> parents: 12802 diff changeset	88 __quadl_recurse_done__ = false;
85dac13a911b quadl.m: Fix integration with large error tolerances (Bug #33792) Rik <octave@nomad.inbox5.com> parents: 12802 diff changeset	89 ## 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	90 global __quadl_need_warning__;
85dac13a911b quadl.m: Fix integration with large error tolerances (Bug #33792) Rik <octave@nomad.inbox5.com> parents: 12802 diff changeset	91 __quadl_need_warning__ = true;
85dac13a911b quadl.m: Fix integration with large error tolerances (Bug #33792) Rik <octave@nomad.inbox5.com> parents: 12802 diff changeset	92
11587 c792872f8942 all script files: untabify and strip trailing whitespace John W. Eaton <jwe@octave.org> parents: 11562 diff changeset	93 m = (a+b)/2;
5837 55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	94 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	95 alpha = sqrt (2/3);
85dac13a911b quadl.m: Fix integration with large error tolerances (Bug #33792) Rik <octave@nomad.inbox5.com> parents: 12802 diff changeset	96 beta = 1/sqrt (5);
5838 376e02b2ce70 [project @ 2006-06-01 20:23:53 by jwe] jwe parents: 5837 diff changeset	97
11587 c792872f8942 all script files: untabify and strip trailing whitespace John W. Eaton <jwe@octave.org> parents: 11562 diff changeset	98 x1 = .942882415695480;
5837 55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	99 x2 = .641853342345781;
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	100 x3 = .236383199662150;
5838 376e02b2ce70 [project @ 2006-06-01 20:23:53 by jwe] jwe parents: 5837 diff changeset	101
376e02b2ce70 [project @ 2006-06-01 20:23:53 by jwe] jwe parents: 5837 diff changeset	102 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	103 m+betah, m+x2h, m+alphah, m+x1h, b];
376e02b2ce70 [project @ 2006-06-01 20:23:53 by jwe] jwe parents: 5837 diff changeset	104
376e02b2ce70 [project @ 2006-06-01 20:23:53 by jwe] jwe parents: 5837 diff changeset	105 y = feval (f, x, varargin{:});
376e02b2ce70 [project @ 2006-06-01 20:23:53 by jwe] jwe parents: 5837 diff changeset	106
11587 c792872f8942 all script files: untabify and strip trailing whitespace John W. Eaton <jwe@octave.org> parents: 11562 diff changeset	107 fa = y(1);
5837 55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	108 fb = y(13);
5838 376e02b2ce70 [project @ 2006-06-01 20:23:53 by jwe] jwe parents: 5837 diff changeset	109
376e02b2ce70 [project @ 2006-06-01 20:23:53 by jwe] jwe parents: 5837 diff changeset	110 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	111
13008 85dac13a911b quadl.m: Fix integration with large error tolerances (Bug #33792) Rik <octave@nomad.inbox5.com> parents: 12802 diff changeset	112 i1 = (h/1470)( 77(y(1)+y(13))
10549 95c3e38098bf Untabify .m scripts Rik <code@nomad.inbox5.com> parents: 9245 diff changeset	113 + 432*(y(3)+y(11))
95c3e38098bf Untabify .m scripts Rik <code@nomad.inbox5.com> parents: 9245 diff changeset	114 + 625*(y(5)+y(9))
95c3e38098bf Untabify .m scripts Rik <code@nomad.inbox5.com> parents: 9245 diff changeset	115 + 672*y(7));
5838 376e02b2ce70 [project @ 2006-06-01 20:23:53 by jwe] jwe parents: 5837 diff changeset	116
13008 85dac13a911b quadl.m: Fix integration with large error tolerances (Bug #33792) Rik <octave@nomad.inbox5.com> parents: 12802 diff changeset	117 is = h( .0158271919734802(y(1)+y(13))
10549 95c3e38098bf Untabify .m scripts Rik <code@nomad.inbox5.com> parents: 9245 diff changeset	118 +.0942738402188500*(y(2)+y(12))
95c3e38098bf Untabify .m scripts Rik <code@nomad.inbox5.com> parents: 9245 diff changeset	119 + .155071987336585*(y(3)+y(11))
95c3e38098bf Untabify .m scripts Rik <code@nomad.inbox5.com> parents: 9245 diff changeset	120 + .188821573960182*(y(4)+y(10))
95c3e38098bf Untabify .m scripts Rik <code@nomad.inbox5.com> parents: 9245 diff changeset	121 + .199773405226859*(y(5)+y(9))
95c3e38098bf Untabify .m scripts Rik <code@nomad.inbox5.com> parents: 9245 diff changeset	122 + .224926465333340*(y(6)+y(8))
95c3e38098bf Untabify .m scripts Rik <code@nomad.inbox5.com> parents: 9245 diff changeset	123 + .242611071901408*y(7));
5838 376e02b2ce70 [project @ 2006-06-01 20:23:53 by jwe] jwe parents: 5837 diff changeset	124
13008 85dac13a911b quadl.m: Fix integration with large error tolerances (Bug #33792) Rik <octave@nomad.inbox5.com> parents: 12802 diff changeset	125 s = sign (is);
5838 376e02b2ce70 [project @ 2006-06-01 20:23:53 by jwe] jwe parents: 5837 diff changeset	126 if (s == 0)
376e02b2ce70 [project @ 2006-06-01 20:23:53 by jwe] jwe parents: 5837 diff changeset	127 s = 1;
5837 55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	128 endif
13008 85dac13a911b quadl.m: Fix integration with large error tolerances (Bug #33792) Rik <octave@nomad.inbox5.com> parents: 12802 diff changeset	129 erri1 = abs (i1-is);
85dac13a911b quadl.m: Fix integration with large error tolerances (Bug #33792) Rik <octave@nomad.inbox5.com> parents: 12802 diff changeset	130 erri2 = abs (i2-is);
5837 55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	131 if (erri2 != 0)
11587 c792872f8942 all script files: untabify and strip trailing whitespace John W. Eaton <jwe@octave.org> parents: 11562 diff changeset	132 R = erri1/erri2;
13008 85dac13a911b quadl.m: Fix integration with large error tolerances (Bug #33792) Rik <octave@nomad.inbox5.com> parents: 12802 diff changeset	133 else
85dac13a911b quadl.m: Fix integration with large error tolerances (Bug #33792) Rik <octave@nomad.inbox5.com> parents: 12802 diff changeset	134 R = 1;
5837 55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	135 endif
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	136 if (R > 0 && R < 1)
11587 c792872f8942 all script files: untabify and strip trailing whitespace John W. Eaton <jwe@octave.org> parents: 11562 diff changeset	137 tol = tol/R;
5837 55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	138 endif
14868 5d3a684236b0 maint: Use Octave coding conventions for cuddling parentheses in scripts directory Rik <octave@nomad.inbox5.com> parents: 14138 diff changeset	139 is = s * abs (is) * tol/myeps;
5837 55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	140 if (is == 0)
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	141 is = b-a;
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	142 endif
13008 85dac13a911b quadl.m: Fix integration with large error tolerances (Bug #33792) Rik <octave@nomad.inbox5.com> parents: 12802 diff changeset	143
12612 16cca721117b doc: Update all documentation for chapter on Numerical Integration Rik <octave@nomad.inbox5.com> parents: 12575 diff changeset	144 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	145
5837 55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	146 endfunction
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	147
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	148 ## ADAPTLOBSTP Recursive function used by QUADL.
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	149 ##
5838 376e02b2ce70 [project @ 2006-06-01 20:23:53 by jwe] jwe parents: 5837 diff changeset	150 ## 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	151 ## 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	152 ## an appropriate relative error. The argument 'F' is
5837 55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	153 ## a string containing the name of f. The remaining
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	154 ## arguments are generated by ADAPTLOB or by recursion.
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	155 ##
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	156 ## Walter Gautschi, 08/03/98
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	157
12612 16cca721117b doc: Update all documentation for chapter on Numerical Integration Rik <octave@nomad.inbox5.com> parents: 12575 diff changeset	158 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	159 global __quadl_recurse_done__;
85dac13a911b quadl.m: Fix integration with large error tolerances (Bug #33792) Rik <octave@nomad.inbox5.com> parents: 12802 diff changeset	160 global __quadl_need_warning__;
85dac13a911b quadl.m: Fix integration with large error tolerances (Bug #33792) Rik <octave@nomad.inbox5.com> parents: 12802 diff changeset	161
11587 c792872f8942 all script files: untabify and strip trailing whitespace John W. Eaton <jwe@octave.org> parents: 11562 diff changeset	162 h = (b-a)/2;
5837 55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	163 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	164 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	165 beta = 1 / sqrt (5);
11587 c792872f8942 all script files: untabify and strip trailing whitespace John W. Eaton <jwe@octave.org> parents: 11562 diff changeset	166 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	167 ml = m-beta*h;
85dac13a911b quadl.m: Fix integration with large error tolerances (Bug #33792) Rik <octave@nomad.inbox5.com> parents: 12802 diff changeset	168 mr = m+beta*h;
5837 55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	169 mrr = m+alpha*h;
5838 376e02b2ce70 [project @ 2006-06-01 20:23:53 by jwe] jwe parents: 5837 diff changeset	170 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	171 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	172 fmll = y(1);
13008 85dac13a911b quadl.m: Fix integration with large error tolerances (Bug #33792) Rik <octave@nomad.inbox5.com> parents: 12802 diff changeset	173 fml = y(2);
85dac13a911b quadl.m: Fix integration with large error tolerances (Bug #33792) Rik <octave@nomad.inbox5.com> parents: 12802 diff changeset	174 fm = y(3);
85dac13a911b quadl.m: Fix integration with large error tolerances (Bug #33792) Rik <octave@nomad.inbox5.com> parents: 12802 diff changeset	175 fmr = y(4);
5837 55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	176 fmrr = y(5);
5838 376e02b2ce70 [project @ 2006-06-01 20:23:53 by jwe] jwe parents: 5837 diff changeset	177 i2 = (h/6)(fa + fb + 5(fml+fmr));
376e02b2ce70 [project @ 2006-06-01 20:23:53 by jwe] jwe parents: 5837 diff changeset	178 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	179 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	180 if ((m <= a \|\| b <= m) && __quadl_need_warning__)
5838 376e02b2ce70 [project @ 2006-06-01 20:23:53 by jwe] jwe parents: 5837 diff changeset	181 warning ("quadl: interval contains no more machine number");
376e02b2ce70 [project @ 2006-06-01 20:23:53 by jwe] jwe parents: 5837 diff changeset	182 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	183 __quadl_need_warning__ = false;
5837 55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	184 endif
12612 16cca721117b doc: Update all documentation for chapter on Numerical Integration Rik <octave@nomad.inbox5.com> parents: 12575 diff changeset	185 q = i1;
5837 55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	186 if (trace)
12612 16cca721117b doc: Update all documentation for chapter on Numerical Integration Rik <octave@nomad.inbox5.com> parents: 12575 diff changeset	187 disp ([a, b-a, q]);
5837 55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	188 endif
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	189 else
13008 85dac13a911b quadl.m: Fix integration with large error tolerances (Bug #33792) Rik <octave@nomad.inbox5.com> parents: 12802 diff changeset	190 __quadl_recurse_done__ = true;
85dac13a911b quadl.m: Fix integration with large error tolerances (Bug #33792) Rik <octave@nomad.inbox5.com> parents: 12802 diff changeset	191 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	192 + 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	193 + 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	194 + 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	195 + 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	196 + adaptlobstp (f, mrr, b , fmrr, fb , is, trace, varargin{:}));
5837 55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	197 endif
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	198 endfunction
12802 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
412882f498b4 codesprint: Wrote 5 tests for quadl.m David Wells <drwells@vt.edu> parents: 12642 diff changeset	201 ## basic functionality
13008 85dac13a911b quadl.m: Fix integration with large error tolerances (Bug #33792) Rik <octave@nomad.inbox5.com> parents: 12802 diff changeset	202 %!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	203
412882f498b4 codesprint: Wrote 5 tests for quadl.m David Wells <drwells@vt.edu> parents: 12642 diff changeset	204 ## 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	205 %!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	206 %! 2.588424538641647e+10, -1.1e-14)
12802 412882f498b4 codesprint: Wrote 5 tests for quadl.m David Wells <drwells@vt.edu> parents: 12642 diff changeset	207
412882f498b4 codesprint: Wrote 5 tests for quadl.m David Wells <drwells@vt.edu> parents: 12642 diff changeset	208 ## extra parameters
412882f498b4 codesprint: Wrote 5 tests for quadl.m David Wells <drwells@vt.edu> parents: 12642 diff changeset	209 %!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	210 %! 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	211
13141 e81ddf9cacd5 maint: untabify and remove trailing whitespace from source files John W. Eaton <jwe@octave.org> parents: 13008 diff changeset	212 ## test different tolerances.
13008 85dac13a911b quadl.m: Fix integration with large error tolerances (Bug #33792) Rik <octave@nomad.inbox5.com> parents: 12802 diff changeset	213 %!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	214 %! (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	215 %!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	216 %! (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	217

Mercurial > hg > octave-nkf

annotate scripts/general/quadl.m @ 19898:4197fc428c7d