annotate scripts/general/quadl.m @ 19040:0850b5212619

doc: Add @nospell macro around proper names in documentation. * aspell-octave.en.pws: Remove proper names from Octave's custom dictionary. * doc/interpreter/audio.txi, doc/interpreter/diagperm.txi, doc/interpreter/diffeq.txi, doc/interpreter/emacs.txi, doc/interpreter/func.txi, doc/interpreter/geometry.txi, doc/interpreter/intro.txi, doc/interpreter/linalg.txi, doc/interpreter/numbers.txi, doc/interpreter/preface.txi, doc/interpreter/quad.txi, doc/interpreter/sparse.txi, libinterp/corefcn/colloc.cc, libinterp/corefcn/ellipj.cc, libinterp/corefcn/gammainc.cc, libinterp/corefcn/hess.cc, libinterp/corefcn/quadcc.cc, libinterp/corefcn/rand.cc, libinterp/corefcn/sqrtm.cc, libinterp/dldfcn/ccolamd.cc, libinterp/dldfcn/colamd.cc, libinterp/dldfcn/dmperm.cc, libinterp/dldfcn/symrcm.cc, scripts/general/quadgk.m, scripts/general/quadl.m, scripts/linear-algebra/commutation_matrix.m, scripts/linear-algebra/condest.m, scripts/linear-algebra/duplication_matrix.m, scripts/linear-algebra/expm.m, scripts/linear-algebra/housh.m, scripts/linear-algebra/krylov.m, scripts/linear-algebra/onenormest.m, scripts/linear-algebra/qzhess.m, scripts/optimization/fminsearch.m, scripts/optimization/glpk.m, scripts/signal/arch_fit.m, scripts/signal/bartlett.m, scripts/signal/blackman.m, scripts/signal/diffpara.m, scripts/signal/durbinlevinson.m, scripts/signal/hamming.m, scripts/signal/hanning.m, scripts/sparse/eigs.m, scripts/sparse/pcg.m, scripts/sparse/pcr.m, scripts/specfun/ellipke.m, scripts/special-matrix/gallery.m, scripts/special-matrix/invhilb.m, scripts/statistics/base/quantile.m, scripts/statistics/tests/kruskal_wallis_test.m, scripts/statistics/tests/manova.m, scripts/statistics/tests/mcnemar_test.m: Add @nospell macro around proper names in documentation.
author Rik <rik@octave.org>
date Thu, 21 Aug 2014 15:57:18 -0700
parents d63878346099
children 4197fc428c7d
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
17744
d63878346099 maint: Update copyright notices for release.
John W. Eaton <jwe@octave.org>
parents: 14868
diff changeset
1 ## Copyright (C) 1998-2013 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 ##
19040
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-x1*h, m-alpha*h, m-x2*h, m-beta*h, m-x3*h, m, m+x3*h, ...
376e02b2ce70 [project @ 2006-06-01 20:23:53 by jwe]
jwe
parents: 5837
diff changeset
103 m+beta*h, m+x2*h, m+alpha*h, m+x1*h, 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 (3*x).*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