octave: scripts/general/quadl.m annotate

annotate scripts/general/quadl.m @ 12802:412882f498b4

codesprint: Wrote 5 tests for quadl.m

author	David Wells <drwells@vt.edu>
date	Sat, 16 Jul 2011 12:52:47 -0500
parents	f96b9b9f141b
children	85dac13a911b

rev	line source
11523 fd0a3ac60b0e update copyright notices John W. Eaton <jwe@octave.org> parents: 10793 diff changeset	1 ## Copyright (C) 1998-2011 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 ##
11587 c792872f8942 all script files: untabify and strip trailing whitespace John W. Eaton <jwe@octave.org> parents: 11562 diff changeset	48 ## Reference: 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
12612 16cca721117b doc: Update all documentation for chapter on Numerical Integration Rik <octave@nomad.inbox5.com> parents: 12575 diff changeset	65 function q = quadl (f, a, b, tol, trace, varargin)
5838 376e02b2ce70 [project @ 2006-06-01 20:23:53 by jwe] jwe parents: 5837 diff changeset	66 need_warning (1);
5837 55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	67 if (nargin < 4)
11587 c792872f8942 all script files: untabify and strip trailing whitespace John W. Eaton <jwe@octave.org> parents: 11562 diff changeset	68 tol = [];
5837 55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	69 endif
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	70 if (nargin < 5)
11587 c792872f8942 all script files: untabify and strip trailing whitespace John W. Eaton <jwe@octave.org> parents: 11562 diff changeset	71 trace = [];
5837 55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	72 endif
7795 df9519e9990c Handle single precision eps values David Bateman <dbateman@free.fr> parents: 7017 diff changeset	73 if (isa (a, "single") \|\| isa (b, "single"))
df9519e9990c Handle single precision eps values David Bateman <dbateman@free.fr> parents: 7017 diff changeset	74 myeps = eps ("single");
df9519e9990c Handle single precision eps values David Bateman <dbateman@free.fr> parents: 7017 diff changeset	75 else
df9519e9990c Handle single precision eps values David Bateman <dbateman@free.fr> parents: 7017 diff changeset	76 myeps = eps;
df9519e9990c Handle single precision eps values David Bateman <dbateman@free.fr> parents: 7017 diff changeset	77 endif
5838 376e02b2ce70 [project @ 2006-06-01 20:23:53 by jwe] jwe parents: 5837 diff changeset	78 if (isempty (tol))
11587 c792872f8942 all script files: untabify and strip trailing whitespace John W. Eaton <jwe@octave.org> parents: 11562 diff changeset	79 tol = myeps;
5837 55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	80 endif
5838 376e02b2ce70 [project @ 2006-06-01 20:23:53 by jwe] jwe parents: 5837 diff changeset	81 if (isempty (trace))
11587 c792872f8942 all script files: untabify and strip trailing whitespace John W. Eaton <jwe@octave.org> parents: 11562 diff changeset	82 trace = 0;
5837 55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	83 endif
7795 df9519e9990c Handle single precision eps values David Bateman <dbateman@free.fr> parents: 7017 diff changeset	84 if (tol < myeps)
df9519e9990c Handle single precision eps values David Bateman <dbateman@free.fr> parents: 7017 diff changeset	85 tol = myeps;
5837 55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	86 endif
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	87
11587 c792872f8942 all script files: untabify and strip trailing whitespace John W. Eaton <jwe@octave.org> parents: 11562 diff changeset	88 m = (a+b)/2;
5837 55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	89 h = (b-a)/2;
11587 c792872f8942 all script files: untabify and strip trailing whitespace John W. Eaton <jwe@octave.org> parents: 11562 diff changeset	90 alpha = sqrt(2/3);
5837 55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	91 beta = 1/sqrt(5);
5838 376e02b2ce70 [project @ 2006-06-01 20:23:53 by jwe] jwe parents: 5837 diff changeset	92
11587 c792872f8942 all script files: untabify and strip trailing whitespace John W. Eaton <jwe@octave.org> parents: 11562 diff changeset	93 x1 = .942882415695480;
5837 55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	94 x2 = .641853342345781;
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	95 x3 = .236383199662150;
5838 376e02b2ce70 [project @ 2006-06-01 20:23:53 by jwe] jwe parents: 5837 diff changeset	96
376e02b2ce70 [project @ 2006-06-01 20:23:53 by jwe] jwe parents: 5837 diff changeset	97 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	98 m+betah, m+x2h, m+alphah, m+x1h, b];
376e02b2ce70 [project @ 2006-06-01 20:23:53 by jwe] jwe parents: 5837 diff changeset	99
376e02b2ce70 [project @ 2006-06-01 20:23:53 by jwe] jwe parents: 5837 diff changeset	100 y = feval (f, x, varargin{:});
376e02b2ce70 [project @ 2006-06-01 20:23:53 by jwe] jwe parents: 5837 diff changeset	101
11587 c792872f8942 all script files: untabify and strip trailing whitespace John W. Eaton <jwe@octave.org> parents: 11562 diff changeset	102 fa = y(1);
5837 55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	103 fb = y(13);
5838 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 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	106
376e02b2ce70 [project @ 2006-06-01 20:23:53 by jwe] jwe parents: 5837 diff changeset	107 i1 = (h/1470)(77(y(1)+y(13))
10549 95c3e38098bf Untabify .m scripts Rik <code@nomad.inbox5.com> parents: 9245 diff changeset	108 + 432*(y(3)+y(11))
95c3e38098bf Untabify .m scripts Rik <code@nomad.inbox5.com> parents: 9245 diff changeset	109 + 625*(y(5)+y(9))
95c3e38098bf Untabify .m scripts Rik <code@nomad.inbox5.com> parents: 9245 diff changeset	110 + 672*y(7));
5838 376e02b2ce70 [project @ 2006-06-01 20:23:53 by jwe] jwe parents: 5837 diff changeset	111
376e02b2ce70 [project @ 2006-06-01 20:23:53 by jwe] jwe parents: 5837 diff changeset	112 is = h(.0158271919734802(y(1)+y(13))
10549 95c3e38098bf Untabify .m scripts Rik <code@nomad.inbox5.com> parents: 9245 diff changeset	113 +.0942738402188500*(y(2)+y(12))
95c3e38098bf Untabify .m scripts Rik <code@nomad.inbox5.com> parents: 9245 diff changeset	114 + .155071987336585*(y(3)+y(11))
95c3e38098bf Untabify .m scripts Rik <code@nomad.inbox5.com> parents: 9245 diff changeset	115 + .188821573960182*(y(4)+y(10))
95c3e38098bf Untabify .m scripts Rik <code@nomad.inbox5.com> parents: 9245 diff changeset	116 + .199773405226859*(y(5)+y(9))
95c3e38098bf Untabify .m scripts Rik <code@nomad.inbox5.com> parents: 9245 diff changeset	117 + .224926465333340*(y(6)+y(8))
95c3e38098bf Untabify .m scripts Rik <code@nomad.inbox5.com> parents: 9245 diff changeset	118 + .242611071901408*y(7));
5838 376e02b2ce70 [project @ 2006-06-01 20:23:53 by jwe] jwe parents: 5837 diff changeset	119
11587 c792872f8942 all script files: untabify and strip trailing whitespace John W. Eaton <jwe@octave.org> parents: 11562 diff changeset	120 s = sign(is);
5838 376e02b2ce70 [project @ 2006-06-01 20:23:53 by jwe] jwe parents: 5837 diff changeset	121
376e02b2ce70 [project @ 2006-06-01 20:23:53 by jwe] jwe parents: 5837 diff changeset	122 if (s == 0)
376e02b2ce70 [project @ 2006-06-01 20:23:53 by jwe] jwe parents: 5837 diff changeset	123 s = 1;
5837 55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	124 endif
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	125 erri1 = abs(i1-is);
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	126 erri2 = abs(i2-is);
11587 c792872f8942 all script files: untabify and strip trailing whitespace John W. Eaton <jwe@octave.org> parents: 11562 diff changeset	127 R = 1;
5837 55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	128 if (erri2 != 0)
11587 c792872f8942 all script files: untabify and strip trailing whitespace John W. Eaton <jwe@octave.org> parents: 11562 diff changeset	129 R = erri1/erri2;
5837 55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	130 endif
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	131 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	132 tol = tol/R;
5837 55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	133 endif
7795 df9519e9990c Handle single precision eps values David Bateman <dbateman@free.fr> parents: 7017 diff changeset	134 is = sabs(is)tol/myeps;
5837 55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	135 if (is == 0)
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	136 is = b-a;
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	137 endif
12612 16cca721117b doc: Update all documentation for chapter on Numerical Integration Rik <octave@nomad.inbox5.com> parents: 12575 diff changeset	138 q = adaptlobstp (f, a, b, fa, fb, is, trace, varargin{:});
5837 55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	139 endfunction
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	140
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	141 ## ADAPTLOBSTP Recursive function used by QUADL.
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	142 ##
5838 376e02b2ce70 [project @ 2006-06-01 20:23:53 by jwe] jwe parents: 5837 diff changeset	143 ## 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	144 ## approximate the integral of F(X) from A to B to
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	145 ## an appropriate relative error. The argument 'F' is
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	146 ## a string containing the name of f. The remaining
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	147 ## arguments are generated by ADAPTLOB or by recursion.
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	148 ##
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	149 ## Walter Gautschi, 08/03/98
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	150
12612 16cca721117b doc: Update all documentation for chapter on Numerical Integration Rik <octave@nomad.inbox5.com> parents: 12575 diff changeset	151 function q = adaptlobstp (f, a, b, fa, fb, is, trace, varargin)
11587 c792872f8942 all script files: untabify and strip trailing whitespace John W. Eaton <jwe@octave.org> parents: 11562 diff changeset	152 h = (b-a)/2;
5837 55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	153 m = (a+b)/2;
11587 c792872f8942 all script files: untabify and strip trailing whitespace John W. Eaton <jwe@octave.org> parents: 11562 diff changeset	154 alpha = sqrt(2/3);
5837 55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	155 beta = 1/sqrt(5);
11587 c792872f8942 all script files: untabify and strip trailing whitespace John W. Eaton <jwe@octave.org> parents: 11562 diff changeset	156 mll = m-alpha*h;
c792872f8942 all script files: untabify and strip trailing whitespace John W. Eaton <jwe@octave.org> parents: 11562 diff changeset	157 ml = m-beta*h;
c792872f8942 all script files: untabify and strip trailing whitespace John W. Eaton <jwe@octave.org> parents: 11562 diff changeset	158 mr = m+beta*h;
5837 55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	159 mrr = m+alpha*h;
5838 376e02b2ce70 [project @ 2006-06-01 20:23:53 by jwe] jwe parents: 5837 diff changeset	160 x = [mll, ml, m, mr, mrr];
376e02b2ce70 [project @ 2006-06-01 20:23:53 by jwe] jwe parents: 5837 diff changeset	161 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	162 fmll = y(1);
c792872f8942 all script files: untabify and strip trailing whitespace John W. Eaton <jwe@octave.org> parents: 11562 diff changeset	163 fml = y(2);
c792872f8942 all script files: untabify and strip trailing whitespace John W. Eaton <jwe@octave.org> parents: 11562 diff changeset	164 fm = y(3);
c792872f8942 all script files: untabify and strip trailing whitespace John W. Eaton <jwe@octave.org> parents: 11562 diff changeset	165 fmr = y(4);
5837 55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	166 fmrr = y(5);
5838 376e02b2ce70 [project @ 2006-06-01 20:23:53 by jwe] jwe parents: 5837 diff changeset	167 i2 = (h/6)(fa + fb + 5(fml+fmr));
376e02b2ce70 [project @ 2006-06-01 20:23:53 by jwe] jwe parents: 5837 diff changeset	168 i1 = (h/1470)(77(fa+fb) + 432(fmll+fmrr) + 625(fml+fmr) + 672*fm);
376e02b2ce70 [project @ 2006-06-01 20:23:53 by jwe] jwe parents: 5837 diff changeset	169 if (is+(i1-i2) == is \|\| mll <= a \|\| b <= mrr)
376e02b2ce70 [project @ 2006-06-01 20:23:53 by jwe] jwe parents: 5837 diff changeset	170 if ((m <= a \|\| b <= m) && need_warning ())
376e02b2ce70 [project @ 2006-06-01 20:23:53 by jwe] jwe parents: 5837 diff changeset	171 warning ("quadl: interval contains no more machine number");
376e02b2ce70 [project @ 2006-06-01 20:23:53 by jwe] jwe parents: 5837 diff changeset	172 warning ("quadl: required tolerance may not be met");
376e02b2ce70 [project @ 2006-06-01 20:23:53 by jwe] jwe parents: 5837 diff changeset	173 need_warning (0);
5837 55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	174 endif
12612 16cca721117b doc: Update all documentation for chapter on Numerical Integration Rik <octave@nomad.inbox5.com> parents: 12575 diff changeset	175 q = i1;
5837 55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	176 if (trace)
12612 16cca721117b doc: Update all documentation for chapter on Numerical Integration Rik <octave@nomad.inbox5.com> parents: 12575 diff changeset	177 disp ([a, b-a, q]);
5837 55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	178 endif
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	179 else
12612 16cca721117b doc: Update all documentation for chapter on Numerical Integration Rik <octave@nomad.inbox5.com> parents: 12575 diff changeset	180 q = (adaptlobstp (f, a, mll, fa, fmll, is, trace, varargin{:})
10549 95c3e38098bf Untabify .m scripts Rik <code@nomad.inbox5.com> parents: 9245 diff changeset	181 + adaptlobstp (f, mll, ml, fmll, fml, is, trace, varargin{:})
95c3e38098bf Untabify .m scripts Rik <code@nomad.inbox5.com> parents: 9245 diff changeset	182 + adaptlobstp (f, ml, m, fml, fm, is, trace, varargin{:})
95c3e38098bf Untabify .m scripts Rik <code@nomad.inbox5.com> parents: 9245 diff changeset	183 + adaptlobstp (f, m, mr, fm, fmr, is, trace, varargin{:})
95c3e38098bf Untabify .m scripts Rik <code@nomad.inbox5.com> parents: 9245 diff changeset	184 + adaptlobstp (f, mr, mrr, fmr, fmrr, is, trace, varargin{:})
95c3e38098bf Untabify .m scripts Rik <code@nomad.inbox5.com> parents: 9245 diff changeset	185 + adaptlobstp (f, mrr, b, fmrr, fb, is, trace, varargin{:}));
5837 55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	186 endif
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	187 endfunction
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	188
5838 376e02b2ce70 [project @ 2006-06-01 20:23:53 by jwe] jwe parents: 5837 diff changeset	189 function r = need_warning (v)
5837 55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	190 persistent w = [];
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	191 if (nargin == 0)
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	192 r = w;
11587 c792872f8942 all script files: untabify and strip trailing whitespace John W. Eaton <jwe@octave.org> parents: 11562 diff changeset	193 else
c792872f8942 all script files: untabify and strip trailing whitespace John W. Eaton <jwe@octave.org> parents: 11562 diff changeset	194 w = v;
5837 55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	195 endif
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	196 endfunction
12802 412882f498b4 codesprint: Wrote 5 tests for quadl.m David Wells <drwells@vt.edu> parents: 12642 diff changeset	197
412882f498b4 codesprint: Wrote 5 tests for quadl.m David Wells <drwells@vt.edu> parents: 12642 diff changeset	198
412882f498b4 codesprint: Wrote 5 tests for quadl.m David Wells <drwells@vt.edu> parents: 12642 diff changeset	199 ## basic functionality
412882f498b4 codesprint: Wrote 5 tests for quadl.m David Wells <drwells@vt.edu> parents: 12642 diff changeset	200 %!assert( quadl (@(x) sin (x), 0, pi, [], []), 2, -3e-16)
412882f498b4 codesprint: Wrote 5 tests for quadl.m David Wells <drwells@vt.edu> parents: 12642 diff changeset	201
412882f498b4 codesprint: Wrote 5 tests for quadl.m David Wells <drwells@vt.edu> parents: 12642 diff changeset	202 ## the values here are very high so it may be unavoidable that this fails
412882f498b4 codesprint: Wrote 5 tests for quadl.m David Wells <drwells@vt.edu> parents: 12642 diff changeset	203 %!assert ( quadl (@(x) sin (3x).cosh (x).*sinh (x),10,15),
412882f498b4 codesprint: Wrote 5 tests for quadl.m David Wells <drwells@vt.edu> parents: 12642 diff changeset	204 %! 2.588424538641647e+10, -9e-15)
412882f498b4 codesprint: Wrote 5 tests for quadl.m David Wells <drwells@vt.edu> parents: 12642 diff changeset	205
412882f498b4 codesprint: Wrote 5 tests for quadl.m David Wells <drwells@vt.edu> parents: 12642 diff changeset	206 ## extra parameters
412882f498b4 codesprint: Wrote 5 tests for quadl.m David Wells <drwells@vt.edu> parents: 12642 diff changeset	207 %!assert (quadl (@(x,a,b) sin (a + b*x), 0, 1, [], [], 2, 3),
412882f498b4 codesprint: Wrote 5 tests for quadl.m David Wells <drwells@vt.edu> parents: 12642 diff changeset	208 %! cos(2)/3 - cos(5)/3, - 3e-16)
412882f498b4 codesprint: Wrote 5 tests for quadl.m David Wells <drwells@vt.edu> parents: 12642 diff changeset	209
412882f498b4 codesprint: Wrote 5 tests for quadl.m David Wells <drwells@vt.edu> parents: 12642 diff changeset	210 ## test different tolerances. This test currently fails for a very high
412882f498b4 codesprint: Wrote 5 tests for quadl.m David Wells <drwells@vt.edu> parents: 12642 diff changeset	211 ## tolerances.
412882f498b4 codesprint: Wrote 5 tests for quadl.m David Wells <drwells@vt.edu> parents: 12642 diff changeset	212 %!assert ( quadl (@(x) sin (2 + 3*x).^2, 0, 10, 0.3, []),
412882f498b4 codesprint: Wrote 5 tests for quadl.m David Wells <drwells@vt.edu> parents: 12642 diff changeset	213 %! (60 + sin(4) - sin(64))/12, -0.3)
412882f498b4 codesprint: Wrote 5 tests for quadl.m David Wells <drwells@vt.edu> parents: 12642 diff changeset	214
412882f498b4 codesprint: Wrote 5 tests for quadl.m David Wells <drwells@vt.edu> parents: 12642 diff changeset	215
412882f498b4 codesprint: Wrote 5 tests for quadl.m David Wells <drwells@vt.edu> parents: 12642 diff changeset	216 ## for lower tolerances the test passes.
412882f498b4 codesprint: Wrote 5 tests for quadl.m David Wells <drwells@vt.edu> parents: 12642 diff changeset	217 %!assert ( quadl (@(x) sin (2 + 3*x).^2, 0, 10, 0.1, []),
412882f498b4 codesprint: Wrote 5 tests for quadl.m David Wells <drwells@vt.edu> parents: 12642 diff changeset	218 %! (60 + sin(4) - sin(64))/12, -0.1)

Mercurial > octave

annotate scripts/general/quadl.m @ 12802:412882f498b4