annotate scripts/general/quadl.m @ 27919:1891570abac8

update Octave Project Developers copyright for the new year In files that have the "Octave Project Developers" copyright notice, update for 2020.
author John W. Eaton <jwe@octave.org>
date Mon, 06 Jan 2020 22:29:51 -0500
parents b442ec6dda5c
children bd51beb6205e
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
27919
1891570abac8 update Octave Project Developers copyright for the new year
John W. Eaton <jwe@octave.org>
parents: 27918
diff changeset
1 ## Copyright (C) 1998-2020 The Octave Project Developers
27918
b442ec6dda5c use centralized file for copyright info for individual contributors
John W. Eaton <jwe@octave.org>
parents: 27800
diff changeset
2 ##
b442ec6dda5c use centralized file for copyright info for individual contributors
John W. Eaton <jwe@octave.org>
parents: 27800
diff changeset
3 ## See the file COPYRIGHT.md in the top-level directory of this distribution
b442ec6dda5c use centralized file for copyright info for individual contributors
John W. Eaton <jwe@octave.org>
parents: 27800
diff changeset
4 ## or <https://octave.org/COPYRIGHT.html/>.
b442ec6dda5c use centralized file for copyright info for individual contributors
John W. Eaton <jwe@octave.org>
parents: 27800
diff changeset
5 ##
5837
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe]
jwe
parents:
diff changeset
6 ##
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe]
jwe
parents:
diff changeset
7 ## This file is part of Octave.
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe]
jwe
parents:
diff changeset
8 ##
24534
194eb4bd202b maint: Update punctuation for GPL v3 license text.
Rik <rik@octave.org>
parents: 24208
diff changeset
9 ## Octave is free software: you can redistribute it and/or modify it
5837
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe]
jwe
parents:
diff changeset
10 ## under the terms of the GNU General Public License as published by
24534
194eb4bd202b maint: Update punctuation for GPL v3 license text.
Rik <rik@octave.org>
parents: 24208
diff changeset
11 ## the Free Software Foundation, either version 3 of the License, or
22755
3a2b891d0b33 maint: Standardize Copyright formatting.
Rik <rik@octave.org>
parents: 22323
diff changeset
12 ## (at your option) any later version.
5837
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe]
jwe
parents:
diff changeset
13 ##
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe]
jwe
parents:
diff changeset
14 ## 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
15 ## WITHOUT ANY WARRANTY; without even the implied warranty of
22755
3a2b891d0b33 maint: Standardize Copyright formatting.
Rik <rik@octave.org>
parents: 22323
diff changeset
16 ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
3a2b891d0b33 maint: Standardize Copyright formatting.
Rik <rik@octave.org>
parents: 22323
diff changeset
17 ## GNU General Public License for more details.
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 ## 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
20 ## along with Octave; see the file COPYING. If not, see
24534
194eb4bd202b maint: Update punctuation for GPL v3 license text.
Rik <rik@octave.org>
parents: 24208
diff changeset
21 ## <https://www.gnu.org/licenses/>.
5837
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe]
jwe
parents:
diff changeset
22
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe]
jwe
parents:
diff changeset
23 ## -*- texinfo -*-
20852
516bb87ea72e 2015 Code Sprint: remove class of function from docstring for all m-files.
Rik <rik@octave.org>
parents: 20689
diff changeset
24 ## @deftypefn {} {@var{q} =} quadl (@var{f}, @var{a}, @var{b})
516bb87ea72e 2015 Code Sprint: remove class of function from docstring for all m-files.
Rik <rik@octave.org>
parents: 20689
diff changeset
25 ## @deftypefnx {} {@var{q} =} quadl (@var{f}, @var{a}, @var{b}, @var{tol})
516bb87ea72e 2015 Code Sprint: remove class of function from docstring for all m-files.
Rik <rik@octave.org>
parents: 20689
diff changeset
26 ## @deftypefnx {} {@var{q} =} quadl (@var{f}, @var{a}, @var{b}, @var{tol}, @var{trace})
516bb87ea72e 2015 Code Sprint: remove class of function from docstring for all m-files.
Rik <rik@octave.org>
parents: 20689
diff changeset
27 ## @deftypefnx {} {@var{q} =} quadl (@var{f}, @var{a}, @var{b}, @var{tol}, @var{trace}, @var{p1}, @var{p2}, @dots{})
516bb87ea72e 2015 Code Sprint: remove class of function from docstring for all m-files.
Rik <rik@octave.org>
parents: 20689
diff changeset
28 ## @deftypefnx {} {[@var{q}, @var{nfun}] =} quadl (@dots{})
5837
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe]
jwe
parents:
diff changeset
29 ##
20158
7503499a252b doc: Update docstrings to have one sentence summary as first line.
Rik <rik@octave.org>
parents: 19697
diff changeset
30 ## Numerically evaluate the integral of @var{f} from @var{a} to @var{b} using
25003
2365c2661b3c doc: Spellcheck documentation ahead of 4.4 release.
Rik <rik@octave.org>
parents: 24546
diff changeset
31 ## an adaptive @nospell{Lobatto} rule.
20158
7503499a252b doc: Update docstrings to have one sentence summary as first line.
Rik <rik@octave.org>
parents: 19697
diff changeset
32 ##
7503499a252b doc: Update docstrings to have one sentence summary as first line.
Rik <rik@octave.org>
parents: 19697
diff changeset
33 ## @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
34 ## 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
35 ## 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
36 ##
16cca721117b doc: Update all documentation for chapter on Numerical Integration
Rik <octave@nomad.inbox5.com>
parents: 12575
diff changeset
37 ## @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
38 ## limits must be finite.
5837
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe]
jwe
parents:
diff changeset
39 ##
20689
655816377845 quadl.m: Overhal function and switch to absolute tolerance.
Rik <rik@octave.org>
parents: 20231
diff changeset
40 ## The optional argument @var{tol} defines the absolute tolerance with which
24545
ba8b828ee4f2 doc: Replace @math{1e^{XXX}} sequences with raw 1eXXX (bug #52827).
Rik <rik@octave.org>
parents: 23219
diff changeset
41 ## to perform the integration. The default value is 1e-6.
5837
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe]
jwe
parents:
diff changeset
42 ##
12612
16cca721117b doc: Update all documentation for chapter on Numerical Integration
Rik <octave@nomad.inbox5.com>
parents: 12575
diff changeset
43 ## 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
44 ## integration interval. If @var{trace} is defined then for each subinterval
20689
655816377845 quadl.m: Overhal function and switch to absolute tolerance.
Rik <rik@octave.org>
parents: 20231
diff changeset
45 ## display: (1) the total number of function evaluations, (2) the left end of
655816377845 quadl.m: Overhal function and switch to absolute tolerance.
Rik <rik@octave.org>
parents: 20231
diff changeset
46 ## the subinterval, (3) the length of the subinterval, (4) the approximation of
655816377845 quadl.m: Overhal function and switch to absolute tolerance.
Rik <rik@octave.org>
parents: 20231
diff changeset
47 ## 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
48 ##
16cca721117b doc: Update all documentation for chapter on Numerical Integration
Rik <octave@nomad.inbox5.com>
parents: 12575
diff changeset
49 ## 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
50 ## @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
51 ## empty matrices ([]).
5837
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe]
jwe
parents:
diff changeset
52 ##
20689
655816377845 quadl.m: Overhal function and switch to absolute tolerance.
Rik <rik@octave.org>
parents: 20231
diff changeset
53 ## The result of the integration is returned in @var{q}.
655816377845 quadl.m: Overhal function and switch to absolute tolerance.
Rik <rik@octave.org>
parents: 20231
diff changeset
54 ##
655816377845 quadl.m: Overhal function and switch to absolute tolerance.
Rik <rik@octave.org>
parents: 20231
diff changeset
55 ## The optional output @var{nfun} indicates the total number of function
655816377845 quadl.m: Overhal function and switch to absolute tolerance.
Rik <rik@octave.org>
parents: 20231
diff changeset
56 ## evaluations performed.
655816377845 quadl.m: Overhal function and switch to absolute tolerance.
Rik <rik@octave.org>
parents: 20231
diff changeset
57 ##
19040
0850b5212619 doc: Add @nospell macro around proper names in documentation.
Rik <rik@octave.org>
parents: 17744
diff changeset
58 ## Reference: @nospell{W. Gander and W. Gautschi}, @cite{Adaptive Quadrature -
27800
5a6a19a4e3da doc: Use Texinfo non-sentence ending periods in citations.
Rik <rik@octave.org>
parents: 26376
diff changeset
59 ## Revisited}, BIT Vol.@: 40, No.@: 1, March 2000, pp.@: 84--101.
25143
13fd0610480f doc: Use https whenever possible in @url entries.
Rik <rik@octave.org>
parents: 25054
diff changeset
60 ## @url{https://www.inf.ethz.ch/personal/gander/}
24208
eec262017c6a maint: Strip trailing whitespace from source files.
John W. Eaton <jwe@octave.org>
parents: 24156
diff changeset
61 ## @seealso{quad, quadv, quadgk, quadcc, trapz, dblquad, triplequad, integral,
24102
c723faa56ab4 Add missing functions integral2.m and integral3.m (bug #52074).
Nicholas R. Jankowski <jankowskin@asme.org>
parents: 24066
diff changeset
62 ## integral2, integral3}
5837
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe]
jwe
parents:
diff changeset
63 ## @end deftypefn
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe]
jwe
parents:
diff changeset
64
20689
655816377845 quadl.m: Overhal function and switch to absolute tolerance.
Rik <rik@octave.org>
parents: 20231
diff changeset
65 ## Original Author: Walter Gautschi
655816377845 quadl.m: Overhal function and switch to absolute tolerance.
Rik <rik@octave.org>
parents: 20231
diff changeset
66 ## Date: 08/03/98
655816377845 quadl.m: Overhal function and switch to absolute tolerance.
Rik <rik@octave.org>
parents: 20231
diff changeset
67 ## Reference: Gander, Computermathematik, Birkhaeuser, 1992.
5837
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe]
jwe
parents:
diff changeset
68
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe]
jwe
parents:
diff changeset
69 ## 2003-08-05 Shai Ayal
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe]
jwe
parents:
diff changeset
70 ## * permission from author to release as GPL
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe]
jwe
parents:
diff changeset
71
20689
655816377845 quadl.m: Overhal function and switch to absolute tolerance.
Rik <rik@octave.org>
parents: 20231
diff changeset
72 function [q, nfun] = quadl (f, a, b, tol = [], trace = false, varargin)
13008
85dac13a911b quadl.m: Fix integration with large error tolerances (Bug #33792)
Rik <octave@nomad.inbox5.com>
parents: 12802
diff changeset
73
85dac13a911b quadl.m: Fix integration with large error tolerances (Bug #33792)
Rik <octave@nomad.inbox5.com>
parents: 12802
diff changeset
74 if (nargin < 3)
85dac13a911b quadl.m: Fix integration with large error tolerances (Bug #33792)
Rik <octave@nomad.inbox5.com>
parents: 12802
diff changeset
75 print_usage ();
5837
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe]
jwe
parents:
diff changeset
76 endif
13008
85dac13a911b quadl.m: Fix integration with large error tolerances (Bug #33792)
Rik <octave@nomad.inbox5.com>
parents: 12802
diff changeset
77
7795
df9519e9990c Handle single precision eps values
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
78 if (isa (a, "single") || isa (b, "single"))
20689
655816377845 quadl.m: Overhal function and switch to absolute tolerance.
Rik <rik@octave.org>
parents: 20231
diff changeset
79 eps = eps ("single");
7795
df9519e9990c Handle single precision eps values
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
80 else
20689
655816377845 quadl.m: Overhal function and switch to absolute tolerance.
Rik <rik@octave.org>
parents: 20231
diff changeset
81 eps = eps ("double");
7795
df9519e9990c Handle single precision eps values
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
82 endif
5838
376e02b2ce70 [project @ 2006-06-01 20:23:53 by jwe]
jwe
parents: 5837
diff changeset
83 if (isempty (tol))
20689
655816377845 quadl.m: Overhal function and switch to absolute tolerance.
Rik <rik@octave.org>
parents: 20231
diff changeset
84 tol = 1e-6;
655816377845 quadl.m: Overhal function and switch to absolute tolerance.
Rik <rik@octave.org>
parents: 20231
diff changeset
85 elseif (! isscalar (tol) || tol < 0)
655816377845 quadl.m: Overhal function and switch to absolute tolerance.
Rik <rik@octave.org>
parents: 20231
diff changeset
86 error ("quadl: TOL must be a scalar >=0");
655816377845 quadl.m: Overhal function and switch to absolute tolerance.
Rik <rik@octave.org>
parents: 20231
diff changeset
87 elseif (tol < eps)
24066
d9b0d8ae734f Update tolerances in BIST tests for quad functions.
Rik <rik@octave.org>
parents: 23220
diff changeset
88 warning ("quadl: TOL specified is smaller than machine precision, using %g",
d9b0d8ae734f Update tolerances in BIST tests for quad functions.
Rik <rik@octave.org>
parents: 23220
diff changeset
89 tol);
20689
655816377845 quadl.m: Overhal function and switch to absolute tolerance.
Rik <rik@octave.org>
parents: 20231
diff changeset
90 tol = eps;
5837
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe]
jwe
parents:
diff changeset
91 endif
5838
376e02b2ce70 [project @ 2006-06-01 20:23:53 by jwe]
jwe
parents: 5837
diff changeset
92 if (isempty (trace))
13008
85dac13a911b quadl.m: Fix integration with large error tolerances (Bug #33792)
Rik <octave@nomad.inbox5.com>
parents: 12802
diff changeset
93 trace = false;
5837
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe]
jwe
parents:
diff changeset
94 endif
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe]
jwe
parents:
diff changeset
95
20689
655816377845 quadl.m: Overhal function and switch to absolute tolerance.
Rik <rik@octave.org>
parents: 20231
diff changeset
96 y = feval (f, [a, b], varargin{:});
655816377845 quadl.m: Overhal function and switch to absolute tolerance.
Rik <rik@octave.org>
parents: 20231
diff changeset
97 nfun = 1;
5838
376e02b2ce70 [project @ 2006-06-01 20:23:53 by jwe]
jwe
parents: 5837
diff changeset
98
11587
c792872f8942 all script files: untabify and strip trailing whitespace
John W. Eaton <jwe@octave.org>
parents: 11562
diff changeset
99 fa = y(1);
20689
655816377845 quadl.m: Overhal function and switch to absolute tolerance.
Rik <rik@octave.org>
parents: 20231
diff changeset
100 fb = y(2);
5838
376e02b2ce70 [project @ 2006-06-01 20:23:53 by jwe]
jwe
parents: 5837
diff changeset
101
20689
655816377845 quadl.m: Overhal function and switch to absolute tolerance.
Rik <rik@octave.org>
parents: 20231
diff changeset
102 h = b - a;
5838
376e02b2ce70 [project @ 2006-06-01 20:23:53 by jwe]
jwe
parents: 5837
diff changeset
103
20689
655816377845 quadl.m: Overhal function and switch to absolute tolerance.
Rik <rik@octave.org>
parents: 20231
diff changeset
104 [q, nfun, hmin] = adaptlobstp (f, a, b, fa, fb, Inf, nfun, abs (h),
655816377845 quadl.m: Overhal function and switch to absolute tolerance.
Rik <rik@octave.org>
parents: 20231
diff changeset
105 tol, trace, varargin{:});
5838
376e02b2ce70 [project @ 2006-06-01 20:23:53 by jwe]
jwe
parents: 5837
diff changeset
106
20689
655816377845 quadl.m: Overhal function and switch to absolute tolerance.
Rik <rik@octave.org>
parents: 20231
diff changeset
107 if (nfun > 10_000)
655816377845 quadl.m: Overhal function and switch to absolute tolerance.
Rik <rik@octave.org>
parents: 20231
diff changeset
108 warning ("quadl: maximum iteration count reached -- possible singular integral");
655816377845 quadl.m: Overhal function and switch to absolute tolerance.
Rik <rik@octave.org>
parents: 20231
diff changeset
109 elseif (any (! isfinite (q(:))))
655816377845 quadl.m: Overhal function and switch to absolute tolerance.
Rik <rik@octave.org>
parents: 20231
diff changeset
110 warning ("quadl: infinite or NaN function evaluations were returned");
655816377845 quadl.m: Overhal function and switch to absolute tolerance.
Rik <rik@octave.org>
parents: 20231
diff changeset
111 elseif (hmin < (b - a) * eps)
655816377845 quadl.m: Overhal function and switch to absolute tolerance.
Rik <rik@octave.org>
parents: 20231
diff changeset
112 warning ("quadl: minimum step size reached -- possible singular integral");
5837
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe]
jwe
parents:
diff changeset
113 endif
13008
85dac13a911b quadl.m: Fix integration with large error tolerances (Bug #33792)
Rik <octave@nomad.inbox5.com>
parents: 12802
diff changeset
114
5837
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe]
jwe
parents:
diff changeset
115 endfunction
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe]
jwe
parents:
diff changeset
116
20689
655816377845 quadl.m: Overhal function and switch to absolute tolerance.
Rik <rik@octave.org>
parents: 20231
diff changeset
117 function [q, nfun, hmin] = adaptlobstp (f, a, b, fa, fb, q0, nfun, hmin,
655816377845 quadl.m: Overhal function and switch to absolute tolerance.
Rik <rik@octave.org>
parents: 20231
diff changeset
118 tol, trace, varargin)
655816377845 quadl.m: Overhal function and switch to absolute tolerance.
Rik <rik@octave.org>
parents: 20231
diff changeset
119 persistent alpha = sqrt (2/3);
655816377845 quadl.m: Overhal function and switch to absolute tolerance.
Rik <rik@octave.org>
parents: 20231
diff changeset
120 persistent beta = 1 / sqrt (5);
5837
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe]
jwe
parents:
diff changeset
121
20689
655816377845 quadl.m: Overhal function and switch to absolute tolerance.
Rik <rik@octave.org>
parents: 20231
diff changeset
122 if (nfun > 10_000)
655816377845 quadl.m: Overhal function and switch to absolute tolerance.
Rik <rik@octave.org>
parents: 20231
diff changeset
123 q = q0;
655816377845 quadl.m: Overhal function and switch to absolute tolerance.
Rik <rik@octave.org>
parents: 20231
diff changeset
124 return;
655816377845 quadl.m: Overhal function and switch to absolute tolerance.
Rik <rik@octave.org>
parents: 20231
diff changeset
125 endif
13008
85dac13a911b quadl.m: Fix integration with large error tolerances (Bug #33792)
Rik <octave@nomad.inbox5.com>
parents: 12802
diff changeset
126
20689
655816377845 quadl.m: Overhal function and switch to absolute tolerance.
Rik <rik@octave.org>
parents: 20231
diff changeset
127 h = (b - a) / 2;
655816377845 quadl.m: Overhal function and switch to absolute tolerance.
Rik <rik@octave.org>
parents: 20231
diff changeset
128 m = (a + b) / 2;
655816377845 quadl.m: Overhal function and switch to absolute tolerance.
Rik <rik@octave.org>
parents: 20231
diff changeset
129 mll = m - alpha*h;
655816377845 quadl.m: Overhal function and switch to absolute tolerance.
Rik <rik@octave.org>
parents: 20231
diff changeset
130 ml = m - beta*h;
655816377845 quadl.m: Overhal function and switch to absolute tolerance.
Rik <rik@octave.org>
parents: 20231
diff changeset
131 mr = m + beta*h;
655816377845 quadl.m: Overhal function and switch to absolute tolerance.
Rik <rik@octave.org>
parents: 20231
diff changeset
132 mrr = m + alpha*h;
5838
376e02b2ce70 [project @ 2006-06-01 20:23:53 by jwe]
jwe
parents: 5837
diff changeset
133 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
134 y = feval (f, x, varargin{:});
20689
655816377845 quadl.m: Overhal function and switch to absolute tolerance.
Rik <rik@octave.org>
parents: 20231
diff changeset
135 nfun += 1;
11587
c792872f8942 all script files: untabify and strip trailing whitespace
John W. Eaton <jwe@octave.org>
parents: 11562
diff changeset
136 fmll = y(1);
13008
85dac13a911b quadl.m: Fix integration with large error tolerances (Bug #33792)
Rik <octave@nomad.inbox5.com>
parents: 12802
diff changeset
137 fml = y(2);
85dac13a911b quadl.m: Fix integration with large error tolerances (Bug #33792)
Rik <octave@nomad.inbox5.com>
parents: 12802
diff changeset
138 fm = y(3);
85dac13a911b quadl.m: Fix integration with large error tolerances (Bug #33792)
Rik <octave@nomad.inbox5.com>
parents: 12802
diff changeset
139 fmr = y(4);
5837
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe]
jwe
parents:
diff changeset
140 fmrr = y(5);
5838
376e02b2ce70 [project @ 2006-06-01 20:23:53 by jwe]
jwe
parents: 5837
diff changeset
141 i2 = (h/6)*(fa + fb + 5*(fml+fmr));
376e02b2ce70 [project @ 2006-06-01 20:23:53 by jwe]
jwe
parents: 5837
diff changeset
142 i1 = (h/1470)*(77*(fa+fb) + 432*(fmll+fmrr) + 625*(fml+fmr) + 672*fm);
20689
655816377845 quadl.m: Overhal function and switch to absolute tolerance.
Rik <rik@octave.org>
parents: 20231
diff changeset
143
655816377845 quadl.m: Overhal function and switch to absolute tolerance.
Rik <rik@octave.org>
parents: 20231
diff changeset
144 if (abs (b - a) < hmin)
655816377845 quadl.m: Overhal function and switch to absolute tolerance.
Rik <rik@octave.org>
parents: 20231
diff changeset
145 hmin = abs (b - a);
655816377845 quadl.m: Overhal function and switch to absolute tolerance.
Rik <rik@octave.org>
parents: 20231
diff changeset
146 endif
655816377845 quadl.m: Overhal function and switch to absolute tolerance.
Rik <rik@octave.org>
parents: 20231
diff changeset
147
655816377845 quadl.m: Overhal function and switch to absolute tolerance.
Rik <rik@octave.org>
parents: 20231
diff changeset
148 if (trace)
655816377845 quadl.m: Overhal function and switch to absolute tolerance.
Rik <rik@octave.org>
parents: 20231
diff changeset
149 disp ([nfun, a, b-a, i1]);
655816377845 quadl.m: Overhal function and switch to absolute tolerance.
Rik <rik@octave.org>
parents: 20231
diff changeset
150 endif
655816377845 quadl.m: Overhal function and switch to absolute tolerance.
Rik <rik@octave.org>
parents: 20231
diff changeset
151
655816377845 quadl.m: Overhal function and switch to absolute tolerance.
Rik <rik@octave.org>
parents: 20231
diff changeset
152 ## Force at least one adaptive step (nfun > 2 test).
655816377845 quadl.m: Overhal function and switch to absolute tolerance.
Rik <rik@octave.org>
parents: 20231
diff changeset
153 if ((abs (i1-i2) < tol || mll <= a || b <= mrr) && nfun > 2)
12612
16cca721117b doc: Update all documentation for chapter on Numerical Integration
Rik <octave@nomad.inbox5.com>
parents: 12575
diff changeset
154 q = i1;
5837
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe]
jwe
parents:
diff changeset
155 else
24156
ce435b70fd94 quadl.m: Return single output when inputs are single (bug #52243).
Rik <rik@octave.org>
parents: 24102
diff changeset
156 q = zeros (6, 1, class (x));
20689
655816377845 quadl.m: Overhal function and switch to absolute tolerance.
Rik <rik@octave.org>
parents: 20231
diff changeset
157 [q(1), nfun, hmin] = adaptlobstp (f, a , mll, fa , fmll, q0/6, nfun, hmin,
655816377845 quadl.m: Overhal function and switch to absolute tolerance.
Rik <rik@octave.org>
parents: 20231
diff changeset
158 tol, trace, varargin{:});
655816377845 quadl.m: Overhal function and switch to absolute tolerance.
Rik <rik@octave.org>
parents: 20231
diff changeset
159 [q(2), nfun, hmin] = adaptlobstp (f, mll, ml , fmll, fml , q0/6, nfun, hmin,
655816377845 quadl.m: Overhal function and switch to absolute tolerance.
Rik <rik@octave.org>
parents: 20231
diff changeset
160 tol, trace, varargin{:});
655816377845 quadl.m: Overhal function and switch to absolute tolerance.
Rik <rik@octave.org>
parents: 20231
diff changeset
161 [q(3), nfun, hmin] = adaptlobstp (f, ml , m , fml , fm , q0/6, nfun, hmin,
655816377845 quadl.m: Overhal function and switch to absolute tolerance.
Rik <rik@octave.org>
parents: 20231
diff changeset
162 tol, trace, varargin{:});
655816377845 quadl.m: Overhal function and switch to absolute tolerance.
Rik <rik@octave.org>
parents: 20231
diff changeset
163 [q(4), nfun, hmin] = adaptlobstp (f, m , mr , fm , fmr , q0/6, nfun, hmin,
655816377845 quadl.m: Overhal function and switch to absolute tolerance.
Rik <rik@octave.org>
parents: 20231
diff changeset
164 tol, trace, varargin{:});
655816377845 quadl.m: Overhal function and switch to absolute tolerance.
Rik <rik@octave.org>
parents: 20231
diff changeset
165 [q(5), nfun, hmin] = adaptlobstp (f, mr , mrr, fmr , fmrr, q0/6, nfun, hmin,
655816377845 quadl.m: Overhal function and switch to absolute tolerance.
Rik <rik@octave.org>
parents: 20231
diff changeset
166 tol, trace, varargin{:});
655816377845 quadl.m: Overhal function and switch to absolute tolerance.
Rik <rik@octave.org>
parents: 20231
diff changeset
167 [q(6), nfun, hmin] = adaptlobstp (f, mrr, b , fmrr, fb , q0/6, nfun, hmin,
655816377845 quadl.m: Overhal function and switch to absolute tolerance.
Rik <rik@octave.org>
parents: 20231
diff changeset
168 tol, trace, varargin{:});
655816377845 quadl.m: Overhal function and switch to absolute tolerance.
Rik <rik@octave.org>
parents: 20231
diff changeset
169 q = sum (q);
5837
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe]
jwe
parents:
diff changeset
170 endif
20689
655816377845 quadl.m: Overhal function and switch to absolute tolerance.
Rik <rik@octave.org>
parents: 20231
diff changeset
171
5837
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe]
jwe
parents:
diff changeset
172 endfunction
12802
412882f498b4 codesprint: Wrote 5 tests for quadl.m
David Wells <drwells@vt.edu>
parents: 12642
diff changeset
173
412882f498b4 codesprint: Wrote 5 tests for quadl.m
David Wells <drwells@vt.edu>
parents: 12642
diff changeset
174
412882f498b4 codesprint: Wrote 5 tests for quadl.m
David Wells <drwells@vt.edu>
parents: 12642
diff changeset
175 ## basic functionality
24066
d9b0d8ae734f Update tolerances in BIST tests for quad functions.
Rik <rik@octave.org>
parents: 23220
diff changeset
176 %!assert (quadl (@(x) sin (x), 0, pi), 2, 1e-6)
12802
412882f498b4 codesprint: Wrote 5 tests for quadl.m
David Wells <drwells@vt.edu>
parents: 12642
diff changeset
177
412882f498b4 codesprint: Wrote 5 tests for quadl.m
David Wells <drwells@vt.edu>
parents: 12642
diff changeset
178 ## the values here are very high so it may be unavoidable that this fails
20689
655816377845 quadl.m: Overhal function and switch to absolute tolerance.
Rik <rik@octave.org>
parents: 20231
diff changeset
179 %!assert (quadl (@(x) sin (3*x).*cosh (x).*sinh (x),10,15, 1e-3),
655816377845 quadl.m: Overhal function and switch to absolute tolerance.
Rik <rik@octave.org>
parents: 20231
diff changeset
180 %! 2.588424538641647e+10, 1e-3)
12802
412882f498b4 codesprint: Wrote 5 tests for quadl.m
David Wells <drwells@vt.edu>
parents: 12642
diff changeset
181
412882f498b4 codesprint: Wrote 5 tests for quadl.m
David Wells <drwells@vt.edu>
parents: 12642
diff changeset
182 ## extra parameters
412882f498b4 codesprint: Wrote 5 tests for quadl.m
David Wells <drwells@vt.edu>
parents: 12642
diff changeset
183 %!assert (quadl (@(x,a,b) sin (a + b*x), 0, 1, [], [], 2, 3),
24066
d9b0d8ae734f Update tolerances in BIST tests for quad functions.
Rik <rik@octave.org>
parents: 23220
diff changeset
184 %! cos(2)/3 - cos(5)/3, 1e-6)
12802
412882f498b4 codesprint: Wrote 5 tests for quadl.m
David Wells <drwells@vt.edu>
parents: 12642
diff changeset
185
13141
e81ddf9cacd5 maint: untabify and remove trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents: 13008
diff changeset
186 ## test different tolerances.
20689
655816377845 quadl.m: Overhal function and switch to absolute tolerance.
Rik <rik@octave.org>
parents: 20231
diff changeset
187 %!test
655816377845 quadl.m: Overhal function and switch to absolute tolerance.
Rik <rik@octave.org>
parents: 20231
diff changeset
188 %! [q, nfun1] = quadl (@(x) sin (2 + 3*x).^2, 0, 10, 0.5, []);
655816377845 quadl.m: Overhal function and switch to absolute tolerance.
Rik <rik@octave.org>
parents: 20231
diff changeset
189 %! assert (q, (60 + sin(4) - sin(64))/12, 0.5);
655816377845 quadl.m: Overhal function and switch to absolute tolerance.
Rik <rik@octave.org>
parents: 20231
diff changeset
190 %! [q, nfun2] = quadl (@(x) sin (2 + 3*x).^2, 0, 10, 0.1, []);
655816377845 quadl.m: Overhal function and switch to absolute tolerance.
Rik <rik@octave.org>
parents: 20231
diff changeset
191 %! assert (q, (60 + sin(4) - sin(64))/12, 0.1);
655816377845 quadl.m: Overhal function and switch to absolute tolerance.
Rik <rik@octave.org>
parents: 20231
diff changeset
192 %! assert (nfun2 > nfun1);
12802
412882f498b4 codesprint: Wrote 5 tests for quadl.m
David Wells <drwells@vt.edu>
parents: 12642
diff changeset
193
24156
ce435b70fd94 quadl.m: Return single output when inputs are single (bug #52243).
Rik <rik@octave.org>
parents: 24102
diff changeset
194 %!test # test single input/output
ce435b70fd94 quadl.m: Return single output when inputs are single (bug #52243).
Rik <rik@octave.org>
parents: 24102
diff changeset
195 %! assert (class (quadl (@sin, 0, 1)), "double");
ce435b70fd94 quadl.m: Return single output when inputs are single (bug #52243).
Rik <rik@octave.org>
parents: 24102
diff changeset
196 %! assert (class (quadl (@sin, single (0), 1)), "single");
ce435b70fd94 quadl.m: Return single output when inputs are single (bug #52243).
Rik <rik@octave.org>
parents: 24102
diff changeset
197 %! assert (class (quadl (@sin, 0, single (1))), "single");
ce435b70fd94 quadl.m: Return single output when inputs are single (bug #52243).
Rik <rik@octave.org>
parents: 24102
diff changeset
198
20689
655816377845 quadl.m: Overhal function and switch to absolute tolerance.
Rik <rik@octave.org>
parents: 20231
diff changeset
199 ## Test input validation
655816377845 quadl.m: Overhal function and switch to absolute tolerance.
Rik <rik@octave.org>
parents: 20231
diff changeset
200 %!error quadl ()
655816377845 quadl.m: Overhal function and switch to absolute tolerance.
Rik <rik@octave.org>
parents: 20231
diff changeset
201 %!error quadl (@sin)
655816377845 quadl.m: Overhal function and switch to absolute tolerance.
Rik <rik@octave.org>
parents: 20231
diff changeset
202 %!error quadl (@sin,1)
24066
d9b0d8ae734f Update tolerances in BIST tests for quad functions.
Rik <rik@octave.org>
parents: 23220
diff changeset
203 %!error <TOL must be a scalar> quadl (@sin,0,1, ones (2,2))
d9b0d8ae734f Update tolerances in BIST tests for quad functions.
Rik <rik@octave.org>
parents: 23220
diff changeset
204 %!error <TOL must be .* .=0> quadl (@sin,0,1, -1)