Mercurial > octave-nkf

annotate scripts/general/quadgk.m @ 11523:fd0a3ac60b0e

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
update copyright notices
author John W. Eaton <jwe@octave.org>
date Fri, 14 Jan 2011 05:47:45 -0500
parents a4f482e66b65
children 1811f4f8b3b5
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
11523
fd0a3ac60b0e update copyright notices
John W. Eaton <jwe@octave.org>
parents: 10846
diff changeset
1 ## Copyright (C) 2008-2011 David Bateman
7771
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
2 ##
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
3 ## This file is part of Octave.
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
4 ##
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
5 ## Octave is free software; you can redistribute it and/or modify it
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
6 ## under the terms of the GNU General Public License as published by
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
7 ## the Free Software Foundation; either version 3 of the License, or (at
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
8 ## your option) any later version.
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
9 ##
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
10 ## Octave is distributed in the hope that it will be useful, but
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
11 ## WITHOUT ANY WARRANTY; without even the implied warranty of
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
12 ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
13 ## General Public License for more details.
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
14 ##
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
15 ## You should have received a copy of the GNU General Public License
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
16 ## along with Octave; see the file COPYING. If not, see
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
17 ## <http://www.gnu.org/licenses/>.
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
18
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
19 ## -*- texinfo -*-
10793
be55736a0783 Grammarcheck the documentation from m-files.
Rik <octave@nomad.inbox5.com>
parents: 10549
diff changeset
20 ## @deftypefn {Function File} {} quadgk (@var{f}, @var{a}, @var{b}, @var{abstol}, @var{trace})
7771
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
21 ## @deftypefnx {Function File} {} quadgk (@var{f}, @var{a}, @var{b}, @var{prop}, @var{val}, @dots{})
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
22 ## @deftypefnx {Function File} {[@var{q}, @var{err}] =} quadgk (@dots{})
9067
8970b4b10e9f Cleanup documentation for quad.texi and diffeq.texi
Rik <rdrider0-list@yahoo.com>
parents: 9051
diff changeset
23 ## Numerically evaluate integral using adaptive Gauss-Konrod quadrature.
7771
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
24 ## The formulation is based on a proposal by L.F. Shampine,
9051
1bf0ce0930be Grammar check TexInfo in all .m files
Rik <rdrider0-list@yahoo.com>
parents: 8920
diff changeset
25 ## @cite{"Vectorized adaptive quadrature in @sc{matlab}", Journal of
7771
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
26 ## Computational and Applied Mathematics, pp131-140, Vol 211, Issue 2,
9067
8970b4b10e9f Cleanup documentation for quad.texi and diffeq.texi
Rik <rdrider0-list@yahoo.com>
parents: 9051
diff changeset
27 ## Feb 2008} where all function evaluations at an iteration are
9051
1bf0ce0930be Grammar check TexInfo in all .m files
Rik <rdrider0-list@yahoo.com>
parents: 8920
diff changeset
28 ## calculated with a single call to @var{f}. Therefore the function
7771
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
29 ## @var{f} must be of the form @code{@var{f} (@var{x})} and accept
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
30 ## vector values of @var{x} and return a vector of the same length
9067
8970b4b10e9f Cleanup documentation for quad.texi and diffeq.texi
Rik <rdrider0-list@yahoo.com>
parents: 9051
diff changeset
31 ## representing the function evaluations at the given values of @var{x}.
7771
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
32 ## The function @var{f} can be defined in terms of a function handle,
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
33 ## inline function or string.
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
34 ##
7989
23c248d415b5 Various doc fixes. Readd cellidx
David Bateman <dbateman@free.fr>
parents: 7795
diff changeset
35 ## The bounds of the quadrature @code{[@var{a}, @var{b}]} can be finite
9051
1bf0ce0930be Grammar check TexInfo in all .m files
Rik <rdrider0-list@yahoo.com>
parents: 8920
diff changeset
36 ## or infinite and contain weak end singularities. Variable
7771
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
37 ## transformation will be used to treat infinite intervals and weaken
10821
693e22af08ae Grammarcheck documentation of m-files
Rik <octave@nomad.inbox5.com>
parents: 10793
diff changeset
38 ## the singularities. For example:
7771
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
39 ##
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
40 ## @example
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
41 ## quadgk(@@(x) 1 ./ (sqrt (x) .* (x + 1)), 0, Inf)
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
42 ## @end example
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
43 ##
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
44 ## @noindent
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
45 ## Note that the formulation of the integrand uses the
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
46 ## element-by-element operator @code{./} and all user functions to
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
47 ## @code{quadgk} should do the same.
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
48 ##
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
49 ## The absolute tolerance can be passed as a fourth argument in a manner
9051
1bf0ce0930be Grammar check TexInfo in all .m files
Rik <rdrider0-list@yahoo.com>
parents: 8920
diff changeset
50 ## compatible with @code{quadv}. Equally the user can request that
7771
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
51 ## information on the convergence can be printed is the fifth argument
9067
8970b4b10e9f Cleanup documentation for quad.texi and diffeq.texi
Rik <rdrider0-list@yahoo.com>
parents: 9051
diff changeset
52 ## is logically true.
7771
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
53 ##
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
54 ## Alternatively, certain properties of @code{quadgk} can be passed as
9051
1bf0ce0930be Grammar check TexInfo in all .m files
Rik <rdrider0-list@yahoo.com>
parents: 8920
diff changeset
55 ## pairs @code{@var{prop}, @var{val}}. Valid properties are
7771
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
56 ##
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
57 ## @table @code
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
58 ## @item AbsTol
9051
1bf0ce0930be Grammar check TexInfo in all .m files
Rik <rdrider0-list@yahoo.com>
parents: 8920
diff changeset
59 ## Defines the absolute error tolerance for the quadrature. The default
7771
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
60 ## absolute tolerance is 1e-10.
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
61 ##
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
62 ## @item RelTol
9051
1bf0ce0930be Grammar check TexInfo in all .m files
Rik <rdrider0-list@yahoo.com>
parents: 8920
diff changeset
63 ## Defines the relative error tolerance for the quadrature. The default
7771
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
64 ## relative tolerance is 1e-5.
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
65 ##
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
66 ## @item MaxIntervalCount
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
67 ## @code{quadgk} initially subdivides the interval on which to perform
9051
1bf0ce0930be Grammar check TexInfo in all .m files
Rik <rdrider0-list@yahoo.com>
parents: 8920
diff changeset
68 ## the quadrature into 10 intervals. Sub-intervals that have an
1bf0ce0930be Grammar check TexInfo in all .m files
Rik <rdrider0-list@yahoo.com>
parents: 8920
diff changeset
69 ## unacceptable error are sub-divided and re-evaluated. If the number of
7771
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
70 ## sub-intervals exceeds at any point 650 sub-intervals then a poor
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
71 ## convergence is signaled and the current estimate of the integral is
9051
1bf0ce0930be Grammar check TexInfo in all .m files
Rik <rdrider0-list@yahoo.com>
parents: 8920
diff changeset
72 ## returned. The property 'MaxIntervalCount' can be used to alter the
7771
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
73 ## number of sub-intervals that can exist before exiting.
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
74 ##
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
75 ## @item WayPoints
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
76 ## If there exists discontinuities in the first derivative of the
8516
e2a179415bac doc fixes
John W. Eaton <jwe@octave.org>
parents: 8507
diff changeset
77 ## function to integrate, then these can be flagged with the
9051
1bf0ce0930be Grammar check TexInfo in all .m files
Rik <rdrider0-list@yahoo.com>
parents: 8920
diff changeset
78 ## @code{"WayPoints"} property. This forces the ends of a sub-interval
8516
e2a179415bac doc fixes
John W. Eaton <jwe@octave.org>
parents: 8507
diff changeset
79 ## to fall on the breakpoints of the function and can result in
9323
8a3afa8ebb06 Fixes for quadgk for doubly infinite interval
David Bateman <dbateman@free.fr>
parents: 9067
diff changeset
80 ## significantly improved estimation of the error in the integral, faster
10846
a4f482e66b65 Grammarcheck more of the documentation.
Rik <octave@nomad.inbox5.com>
parents: 10821
diff changeset
81 ## computation or both. For example,
7771
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
82 ##
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
83 ## @example
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
84 ## quadgk (@@(x) abs (1 - x .^ 2), 0, 2, 'Waypoints', 1)
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
85 ## @end example
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
86 ##
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
87 ## @noindent
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
88 ## signals the breakpoint in the integrand at @code{@var{x} = 1}.
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
89 ##
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
90 ## @item Trace
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
91 ## If logically true, then @code{quadgk} prints information on the
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
92 ## convergence of the quadrature at each iteration.
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
93 ##@end table
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
94 ##
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
95 ## If any of @var{a}, @var{b} or @var{waypoints} is complex, then the
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
96 ## quadrature is treated as a contour integral along a piecewise
9051
1bf0ce0930be Grammar check TexInfo in all .m files
Rik <rdrider0-list@yahoo.com>
parents: 8920
diff changeset
97 ## continuous path defined by the above. In this case the integral is
10846
a4f482e66b65 Grammarcheck more of the documentation.
Rik <octave@nomad.inbox5.com>
parents: 10821
diff changeset
98 ## assumed to have no edge singularities. For example,
7771
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
99 ##
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
100 ## @example
9051
1bf0ce0930be Grammar check TexInfo in all .m files
Rik <rdrider0-list@yahoo.com>
parents: 8920
diff changeset
101 ## @group
8516
e2a179415bac doc fixes
John W. Eaton <jwe@octave.org>
parents: 8507
diff changeset
102 ## quadgk (@@(z) log (z), 1+1i, 1+1i, "WayPoints",
e2a179415bac doc fixes
John W. Eaton <jwe@octave.org>
parents: 8507
diff changeset
103 ## [1-1i, -1,-1i, -1+1i])
9051
1bf0ce0930be Grammar check TexInfo in all .m files
Rik <rdrider0-list@yahoo.com>
parents: 8920
diff changeset
104 ## @end group
7771
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
105 ## @end example
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
106 ##
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
107 ## @noindent
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
108 ## integrates @code{log (z)} along the square defined by @code{[1+1i,
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
109 ## 1-1i, -1-1i, -1+1i]}
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
110 ##
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
111 ## If two output arguments are requested, then @var{err} returns the
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
112 ## approximate bounds on the error in the integral @code{abs (@var{q} -
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
113 ## @var{i})}, where @var{i} is the exact value of the integral.
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
114 ##
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
115 ## @seealso{triplequad, dblquad, quad, quadl, quadv, trapz}
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
116 ## @end deftypefn
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
117
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
118 function [q, err] = quadgk (f, a, b, varargin)
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
119 if (nargin < 3)
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
120 print_usage ();
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
121 endif
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
122
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
123 if (b < a)
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
124 [q, err] = quadgk (f, b, a, varargin{:});
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
125 q = -q;
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
126 else
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
127 abstol = 1e-10;
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
128 reltol = 1e-5;
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
129 waypoints = [];
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
130 maxint = 650;
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
131 trace = false;
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
132
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
133 if (nargin > 3)
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
134 if (! ischar (varargin{1}))
10549
95c3e38098bf Untabify .m scripts
Rik <code@nomad.inbox5.com>
parents: 9590
diff changeset
135 if (!isempty (varargin{1}))
95c3e38098bf Untabify .m scripts
Rik <code@nomad.inbox5.com>
parents: 9590
diff changeset
136 abstol = varargin{1};
95c3e38098bf Untabify .m scripts
Rik <code@nomad.inbox5.com>
parents: 9590
diff changeset
137 reltol = 0;
95c3e38098bf Untabify .m scripts
Rik <code@nomad.inbox5.com>
parents: 9590
diff changeset
138 endif
95c3e38098bf Untabify .m scripts
Rik <code@nomad.inbox5.com>
parents: 9590
diff changeset
139 if (nargin > 4)
95c3e38098bf Untabify .m scripts
Rik <code@nomad.inbox5.com>
parents: 9590
diff changeset
140 trace = varargin{2};
95c3e38098bf Untabify .m scripts
Rik <code@nomad.inbox5.com>
parents: 9590
diff changeset
141 endif
95c3e38098bf Untabify .m scripts
Rik <code@nomad.inbox5.com>
parents: 9590
diff changeset
142 if (nargin > 5)
95c3e38098bf Untabify .m scripts
Rik <code@nomad.inbox5.com>
parents: 9590
diff changeset
143 error ("quadgk: can not pass additional arguments to user function");
7771
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
144 endif
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
145 else
10549
95c3e38098bf Untabify .m scripts
Rik <code@nomad.inbox5.com>
parents: 9590
diff changeset
146 idx = 1;
95c3e38098bf Untabify .m scripts
Rik <code@nomad.inbox5.com>
parents: 9590
diff changeset
147 while (idx < nargin - 3)
95c3e38098bf Untabify .m scripts
Rik <code@nomad.inbox5.com>
parents: 9590
diff changeset
148 if (ischar (varargin{idx}))
95c3e38098bf Untabify .m scripts
Rik <code@nomad.inbox5.com>
parents: 9590
diff changeset
149 str = varargin{idx++};
95c3e38098bf Untabify .m scripts
Rik <code@nomad.inbox5.com>
parents: 9590
diff changeset
150 if (strcmpi (str, "reltol"))
95c3e38098bf Untabify .m scripts
Rik <code@nomad.inbox5.com>
parents: 9590
diff changeset
151 reltol = varargin{idx++};
95c3e38098bf Untabify .m scripts
Rik <code@nomad.inbox5.com>
parents: 9590
diff changeset
152 elseif (strcmpi (str, "abstol"))
95c3e38098bf Untabify .m scripts
Rik <code@nomad.inbox5.com>
parents: 9590
diff changeset
153 abstol = varargin{idx++};
95c3e38098bf Untabify .m scripts
Rik <code@nomad.inbox5.com>
parents: 9590
diff changeset
154 elseif (strcmpi (str, "waypoints"))
95c3e38098bf Untabify .m scripts
Rik <code@nomad.inbox5.com>
parents: 9590
diff changeset
155 waypoints = varargin{idx++} (:);
95c3e38098bf Untabify .m scripts
Rik <code@nomad.inbox5.com>
parents: 9590
diff changeset
156 if (isreal(waypoints))
95c3e38098bf Untabify .m scripts
Rik <code@nomad.inbox5.com>
parents: 9590
diff changeset
157 waypoints (waypoints < a | waypoints > b) = [];
95c3e38098bf Untabify .m scripts
Rik <code@nomad.inbox5.com>
parents: 9590
diff changeset
158 endif
95c3e38098bf Untabify .m scripts
Rik <code@nomad.inbox5.com>
parents: 9590
diff changeset
159 elseif (strcmpi (str, "maxintervalcount"))
95c3e38098bf Untabify .m scripts
Rik <code@nomad.inbox5.com>
parents: 9590
diff changeset
160 maxint = varargin{idx++};
95c3e38098bf Untabify .m scripts
Rik <code@nomad.inbox5.com>
parents: 9590
diff changeset
161 elseif (strcmpi (str, "trace"))
95c3e38098bf Untabify .m scripts
Rik <code@nomad.inbox5.com>
parents: 9590
diff changeset
162 trace = varargin{idx++};
95c3e38098bf Untabify .m scripts
Rik <code@nomad.inbox5.com>
parents: 9590
diff changeset
163 else
95c3e38098bf Untabify .m scripts
Rik <code@nomad.inbox5.com>
parents: 9590
diff changeset
164 error ("quadgk: unknown property %s", str);
95c3e38098bf Untabify .m scripts
Rik <code@nomad.inbox5.com>
parents: 9590
diff changeset
165 endif
95c3e38098bf Untabify .m scripts
Rik <code@nomad.inbox5.com>
parents: 9590
diff changeset
166 else
95c3e38098bf Untabify .m scripts
Rik <code@nomad.inbox5.com>
parents: 9590
diff changeset
167 error ("quadgk: expecting property to be a string");
95c3e38098bf Untabify .m scripts
Rik <code@nomad.inbox5.com>
parents: 9590
diff changeset
168 endif
95c3e38098bf Untabify .m scripts
Rik <code@nomad.inbox5.com>
parents: 9590
diff changeset
169 endwhile
95c3e38098bf Untabify .m scripts
Rik <code@nomad.inbox5.com>
parents: 9590
diff changeset
170 if (idx != nargin - 2)
95c3e38098bf Untabify .m scripts
Rik <code@nomad.inbox5.com>
parents: 9590
diff changeset
171 error ("quadgk: expecting properties in pairs");
95c3e38098bf Untabify .m scripts
Rik <code@nomad.inbox5.com>
parents: 9590
diff changeset
172 endif
7771
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
173 endif
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
174 endif
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
175
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
176 ## Convert function given as a string to a function handle
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
177 if (ischar (f))
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
178 f = @(x) feval (f, x);
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
179 endif
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
180
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
181 ## Use variable subsitution to weaken endpoint singularities and to
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
182 ## perform integration with endpoints at infinity. No transform for
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
183 ## contour integrals
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
184 if (iscomplex (a) || iscomplex (b) || iscomplex(waypoints))
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
185 ## contour integral, no transform
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
186 subs = [a; waypoints; b];
9590
b18a50c56144 More care with contour integrals in quadgk
David Bateman <dbateman@free.fr>
parents: 9326
diff changeset
187 h = sum (abs (diff (subs)));
b18a50c56144 More care with contour integrals in quadgk
David Bateman <dbateman@free.fr>
parents: 9326
diff changeset
188 h0 = h;
7771
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
189 trans = @(t) t;
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
190 elseif (isinf (a) && isinf(b))
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
191 ## Standard Infinite to finite integral transformation.
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
192 ## \int_{-\infinity_^\infinity f(x) dx = \int_-1^1 f (g(t)) g'(t) dt
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
193 ## where
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
194 ## g(t) = t / (1 - t^2)
9323
8a3afa8ebb06 Fixes for quadgk for doubly infinite interval
David Bateman <dbateman@free.fr>
parents: 9067
diff changeset
195 ## g'(t) = (1 + t^2) / (1 - t^2) ^ 2
7771
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
196 ## waypoint transform is then
9325
22a4433b41e3 Better waypoint transform pour quadgk
Marco Caliari <marco.caliari@univr.it>
parents: 9323
diff changeset
197 ## t = (2 * g(t)) ./ (1 + sqrt(1 + 4 * g(t) .^ 2))
7771
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
198 if (!isempty (waypoints))
10549
95c3e38098bf Untabify .m scripts
Rik <code@nomad.inbox5.com>
parents: 9590
diff changeset
199 trans = @(x) (2 * x) ./ (1 + sqrt(1 + 4 * x .^ 2));
95c3e38098bf Untabify .m scripts
Rik <code@nomad.inbox5.com>
parents: 9590
diff changeset
200 subs = [-1; trans(waypoints); 1];
7771
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
201 else
10549
95c3e38098bf Untabify .m scripts
Rik <code@nomad.inbox5.com>
parents: 9590
diff changeset
202 subs = linspace (-1, 1, 11)';
7771
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
203 endif
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
204 h = 2;
9590
b18a50c56144 More care with contour integrals in quadgk
David Bateman <dbateman@free.fr>
parents: 9326
diff changeset
205 h0 = b - a;
7771
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
206 trans = @(t) t ./ (1 - t.^2);
9325
22a4433b41e3 Better waypoint transform pour quadgk
Marco Caliari <marco.caliari@univr.it>
parents: 9323
diff changeset
207 f = @(t) f (t ./ (1 - t .^ 2)) .* (1 + t .^ 2) ./ ((1 - t .^ 2) .^ 2);
7771
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
208 elseif (isinf(a))
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
209 ## Formula defined in Shampine paper as two separate steps. One to
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
210 ## weaken singularity at finite end, then a second to transform to
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
211 ## a finite interval. The singularity weakening transform is
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
212 ## \int_{-\infinity}^b f(x) dx =
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
213 ## - \int_{-\infinity}^0 f (b - t^2) 2 t dt
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
214 ## (note minus sign) and the finite interval transform is
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
215 ## \int_{-\infinity}^0 f(b - t^2) 2 t dt =
8507
cadc73247d65 style fixes
John W. Eaton <jwe@octave.org>
parents: 8122
diff changeset
216 ## \int_{-1}^0 f (b - g(s) ^ 2) 2 g(s) g'(s) ds
7771
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
217 ## where
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
218 ## g(s) = s / (1 + s)
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
219 ## g'(s) = 1 / (1 + s) ^ 2
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
220 ## waypoint transform is then
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
221 ## t = sqrt (b - x)
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
222 ## s = - t / (t + 1)
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
223 if (!isempty (waypoints))
10549
95c3e38098bf Untabify .m scripts
Rik <code@nomad.inbox5.com>
parents: 9590
diff changeset
224 tmp = sqrt (b - waypoints);
95c3e38098bf Untabify .m scripts
Rik <code@nomad.inbox5.com>
parents: 9590
diff changeset
225 trans = @(x) - x ./ (x + 1);
95c3e38098bf Untabify .m scripts
Rik <code@nomad.inbox5.com>
parents: 9590
diff changeset
226 subs = [0; trans(tmp); 1];
7771
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
227 else
10549
95c3e38098bf Untabify .m scripts
Rik <code@nomad.inbox5.com>
parents: 9590
diff changeset
228 subs = linspace (0, 1, 11)';
7771
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
229 endif
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
230 h = 1;
9590
b18a50c56144 More care with contour integrals in quadgk
David Bateman <dbateman@free.fr>
parents: 9326
diff changeset
231 h0 = b - a;
7771
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
232 trans = @(t) b - (t ./ (1 + t)).^2;
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
233 f = @(s) - 2 * s .* f (b - (s ./ (1 + s)) .^ 2) ./ ((1 + s) .^ 3);
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
234 elseif (isinf(b))
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
235 ## Formula defined in Shampine paper as two separate steps. One to
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
236 ## weaken singularity at finite end, then a second to transform to
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
237 ## a finite interval. The singularity weakening transform is
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
238 ## \int_a^\infinity f(x) dx = \int_0^\infinity f (a + t^2) 2 t dt
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
239 ## and the finite interval transform is
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
240 ## \int_0^\infinity f(a + t^2) 2 t dt =
8507
cadc73247d65 style fixes
John W. Eaton <jwe@octave.org>
parents: 8122
diff changeset
241 ## \int_0^1 f (a + g(s) ^ 2) 2 g(s) g'(s) ds
7771
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
242 ## where
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
243 ## g(s) = s / (1 - s)
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
244 ## g'(s) = 1 / (1 - s) ^ 2
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
245 ## waypoint transform is then
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
246 ## t = sqrt (x - a)
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
247 ## s = t / (t + 1)
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
248 if (!isempty (waypoints))
10549
95c3e38098bf Untabify .m scripts
Rik <code@nomad.inbox5.com>
parents: 9590
diff changeset
249 tmp = sqrt (waypoints - a);
95c3e38098bf Untabify .m scripts
Rik <code@nomad.inbox5.com>
parents: 9590
diff changeset
250 trans = @(x) x ./ (x + 1);
95c3e38098bf Untabify .m scripts
Rik <code@nomad.inbox5.com>
parents: 9590
diff changeset
251 subs = [0; trans(tmp); 1];
7771
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
252 else
10549
95c3e38098bf Untabify .m scripts
Rik <code@nomad.inbox5.com>
parents: 9590
diff changeset
253 subs = linspace (0, 1, 11)';
7771
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
254 endif
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
255 h = 1;
9590
b18a50c56144 More care with contour integrals in quadgk
David Bateman <dbateman@free.fr>
parents: 9326
diff changeset
256 h0 = b - a;
7771
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
257 trans = @(t) a + (t ./ (1 - t)).^2;
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
258 f = @(s) 2 * s .* f (a + (s ./ (1 - s)) .^ 2) ./ ((1 - s) .^ 3);
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
259 else
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
260 ## Davis, Rabinowitz, "Methods of Numerical Integration" p441 2ed.
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
261 ## Presented in section 5 of the Shampine paper as
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
262 ## g(t) = ((b - a) / 2) * (t / 2 * (3 - t^2)) + (b + a) / 2
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
263 ## g'(t) = ((b-a)/4) * (3 - 3t^2);
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
264 ## waypoint transform can then be found by solving for t with
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
265 ## Maxima (solve (c + 3*t - 3^3, t);). This gives 3 roots, two of
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
266 ## which are complex for values between a and b and so can be
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
267 ## ignored. The third is
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
268 ## c = (-4*x + 2*(b+a)) / (b-a);
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
269 ## k = ((sqrt(c^2 - 4) + c)/2)^(1/3);
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
270 ## t = (sqrt(3)* 1i * (1 - k^2) - (1 + k^2)) / 2 / k;
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
271 if (! isempty (waypoints))
10549
95c3e38098bf Untabify .m scripts
Rik <code@nomad.inbox5.com>
parents: 9590
diff changeset
272 trans = @__quadgk_finite_waypoint__;
95c3e38098bf Untabify .m scripts
Rik <code@nomad.inbox5.com>
parents: 9590
diff changeset
273 subs = [-1; trans(waypoints, a, b); 1];
7771
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
274 else
10549
95c3e38098bf Untabify .m scripts
Rik <code@nomad.inbox5.com>
parents: 9590
diff changeset
275 subs = linspace(-1, 1, 11)';
7771
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
276 endif
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
277 h = 2;
9590
b18a50c56144 More care with contour integrals in quadgk
David Bateman <dbateman@free.fr>
parents: 9326
diff changeset
278 h0 = b - a;
7771
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
279 trans = @(t) ((b - a) ./ 4) * t .* (3 - t.^2) + (b + a) ./ 2;
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
280 f = @(t) f((b - a) ./ 4 .* t .* (3 - t.^2) + (b + a) ./ 2) .* ...
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
281 3 .* (b - a) ./ 4 .* (1 - t.^2);
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
282 endif
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
283
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
284 ## Split interval into at least 10 sub-interval with a 15 point
9067
8970b4b10e9f Cleanup documentation for quad.texi and diffeq.texi
Rik <rdrider0-list@yahoo.com>
parents: 9051
diff changeset
285 ## Gauss-Kronrod rule giving a minimum of 150 function evaluations
7771
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
286 while (length (subs) < 11)
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
287 subs = [subs' ; subs(1:end-1)' + diff(subs') ./ 2, NaN](:)(1 : end - 1);
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
288 endwhile
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
289 subs = [subs(1:end-1), subs(2:end)];
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
290
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
291 warn_state = warning ("query", "Octave:divide-by-zero");
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
292
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
293 unwind_protect
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
294 ## Singularity will cause divide by zero warnings
8664
e07e93c04080 style fixes
John W. Eaton <jwe@octave.org>
parents: 8516
diff changeset
295 warning ("off", "Octave:divide-by-zero");
7771
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
296
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
297 ## Initial evaluation of the integrand on the sub-intervals
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
298 [q_subs, q_errs] = __quadgk_eval__ (f, subs);
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
299 q0 = sum (q_subs);
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
300 err0 = sum (q_errs);
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
301
7795
df9519e9990c Handle single precision eps values
David Bateman <dbateman@free.fr>
parents: 7771
diff changeset
302 if (isa (a, "single") || isa (b, "single") || isa (waypoints, "single"))
10549
95c3e38098bf Untabify .m scripts
Rik <code@nomad.inbox5.com>
parents: 9590
diff changeset
303 myeps = eps ("single");
7795
df9519e9990c Handle single precision eps values
David Bateman <dbateman@free.fr>
parents: 7771
diff changeset
304 else
10549
95c3e38098bf Untabify .m scripts
Rik <code@nomad.inbox5.com>
parents: 9590
diff changeset
305 myeps = eps;
7795
df9519e9990c Handle single precision eps values
David Bateman <dbateman@free.fr>
parents: 7771
diff changeset
306 endif
df9519e9990c Handle single precision eps values
David Bateman <dbateman@free.fr>
parents: 7771
diff changeset
307
7771
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
308 first = true;
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
309 while (true)
10549
95c3e38098bf Untabify .m scripts
Rik <code@nomad.inbox5.com>
parents: 9590
diff changeset
310 ## Check for sub-intervals that are too small. Test must be
95c3e38098bf Untabify .m scripts
Rik <code@nomad.inbox5.com>
parents: 9590
diff changeset
311 ## performed in untransformed sub-intervals. What is a good
95c3e38098bf Untabify .m scripts
Rik <code@nomad.inbox5.com>
parents: 9590
diff changeset
312 ## value for this test. Shampine suggests 100*eps
95c3e38098bf Untabify .m scripts
Rik <code@nomad.inbox5.com>
parents: 9590
diff changeset
313 if (any (abs (diff (trans (subs), [], 2) / h0) < 100 * myeps))
95c3e38098bf Untabify .m scripts
Rik <code@nomad.inbox5.com>
parents: 9590
diff changeset
314 q = q0;
95c3e38098bf Untabify .m scripts
Rik <code@nomad.inbox5.com>
parents: 9590
diff changeset
315 err = err0;
95c3e38098bf Untabify .m scripts
Rik <code@nomad.inbox5.com>
parents: 9590
diff changeset
316 break;
95c3e38098bf Untabify .m scripts
Rik <code@nomad.inbox5.com>
parents: 9590
diff changeset
317 endif
7771
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
318
10549
95c3e38098bf Untabify .m scripts
Rik <code@nomad.inbox5.com>
parents: 9590
diff changeset
319 ## Quit if any evaluations are not finite (Inf or NaN)
95c3e38098bf Untabify .m scripts
Rik <code@nomad.inbox5.com>
parents: 9590
diff changeset
320 if (any (! isfinite (q_subs)))
95c3e38098bf Untabify .m scripts
Rik <code@nomad.inbox5.com>
parents: 9590
diff changeset
321 warning ("quadgk: non finite integrand encountered");
95c3e38098bf Untabify .m scripts
Rik <code@nomad.inbox5.com>
parents: 9590
diff changeset
322 q = q0;
95c3e38098bf Untabify .m scripts
Rik <code@nomad.inbox5.com>
parents: 9590
diff changeset
323 err = err0;
95c3e38098bf Untabify .m scripts
Rik <code@nomad.inbox5.com>
parents: 9590
diff changeset
324 break;
95c3e38098bf Untabify .m scripts
Rik <code@nomad.inbox5.com>
parents: 9590
diff changeset
325 endif
7771
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
326
10549
95c3e38098bf Untabify .m scripts
Rik <code@nomad.inbox5.com>
parents: 9590
diff changeset
327 tol = max (abstol, reltol .* abs (q0));
7771
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
328
10549
95c3e38098bf Untabify .m scripts
Rik <code@nomad.inbox5.com>
parents: 9590
diff changeset
329 ## If the global error estimate is meet exit
95c3e38098bf Untabify .m scripts
Rik <code@nomad.inbox5.com>
parents: 9590
diff changeset
330 if (err0 < tol)
95c3e38098bf Untabify .m scripts
Rik <code@nomad.inbox5.com>
parents: 9590
diff changeset
331 q = q0;
95c3e38098bf Untabify .m scripts
Rik <code@nomad.inbox5.com>
parents: 9590
diff changeset
332 err = err0;
95c3e38098bf Untabify .m scripts
Rik <code@nomad.inbox5.com>
parents: 9590
diff changeset
333 break;
95c3e38098bf Untabify .m scripts
Rik <code@nomad.inbox5.com>
parents: 9590
diff changeset
334 endif
7771
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
335
10549
95c3e38098bf Untabify .m scripts
Rik <code@nomad.inbox5.com>
parents: 9590
diff changeset
336 ## Accept the sub-intervals that meet the convergence criteria
95c3e38098bf Untabify .m scripts
Rik <code@nomad.inbox5.com>
parents: 9590
diff changeset
337 idx = find (abs (q_errs) < tol .* abs(diff (subs, [], 2)) ./ h);
95c3e38098bf Untabify .m scripts
Rik <code@nomad.inbox5.com>
parents: 9590
diff changeset
338 if (first)
95c3e38098bf Untabify .m scripts
Rik <code@nomad.inbox5.com>
parents: 9590
diff changeset
339 q = sum (q_subs (idx));
95c3e38098bf Untabify .m scripts
Rik <code@nomad.inbox5.com>
parents: 9590
diff changeset
340 err = sum (q_errs(idx));
95c3e38098bf Untabify .m scripts
Rik <code@nomad.inbox5.com>
parents: 9590
diff changeset
341 first = false;
95c3e38098bf Untabify .m scripts
Rik <code@nomad.inbox5.com>
parents: 9590
diff changeset
342 else
95c3e38098bf Untabify .m scripts
Rik <code@nomad.inbox5.com>
parents: 9590
diff changeset
343 q0 = q + sum (q_subs);
95c3e38098bf Untabify .m scripts
Rik <code@nomad.inbox5.com>
parents: 9590
diff changeset
344 err0 = err + sum (q_errs);
95c3e38098bf Untabify .m scripts
Rik <code@nomad.inbox5.com>
parents: 9590
diff changeset
345 q += sum (q_subs (idx));
95c3e38098bf Untabify .m scripts
Rik <code@nomad.inbox5.com>
parents: 9590
diff changeset
346 err += sum (q_errs(idx));
95c3e38098bf Untabify .m scripts
Rik <code@nomad.inbox5.com>
parents: 9590
diff changeset
347 endif
95c3e38098bf Untabify .m scripts
Rik <code@nomad.inbox5.com>
parents: 9590
diff changeset
348 subs(idx,:) = [];
7771
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
349
10549
95c3e38098bf Untabify .m scripts
Rik <code@nomad.inbox5.com>
parents: 9590
diff changeset
350 ## If no remaining sub-intervals exit
95c3e38098bf Untabify .m scripts
Rik <code@nomad.inbox5.com>
parents: 9590
diff changeset
351 if (rows (subs) == 0)
95c3e38098bf Untabify .m scripts
Rik <code@nomad.inbox5.com>
parents: 9590
diff changeset
352 break;
95c3e38098bf Untabify .m scripts
Rik <code@nomad.inbox5.com>
parents: 9590
diff changeset
353 endif
7771
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
354
10549
95c3e38098bf Untabify .m scripts
Rik <code@nomad.inbox5.com>
parents: 9590
diff changeset
355 if (trace)
95c3e38098bf Untabify .m scripts
Rik <code@nomad.inbox5.com>
parents: 9590
diff changeset
356 disp([rows(subs), err, q0]);
95c3e38098bf Untabify .m scripts
Rik <code@nomad.inbox5.com>
parents: 9590
diff changeset
357 endif
7771
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
358
10549
95c3e38098bf Untabify .m scripts
Rik <code@nomad.inbox5.com>
parents: 9590
diff changeset
359 ## Split remaining sub-intervals in two
95c3e38098bf Untabify .m scripts
Rik <code@nomad.inbox5.com>
parents: 9590
diff changeset
360 mid = (subs(:,2) + subs(:,1)) ./ 2;
95c3e38098bf Untabify .m scripts
Rik <code@nomad.inbox5.com>
parents: 9590
diff changeset
361 subs = [subs(:,1), mid; mid, subs(:,2)];
7771
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
362
10549
95c3e38098bf Untabify .m scripts
Rik <code@nomad.inbox5.com>
parents: 9590
diff changeset
363 ## If the maximum sub-interval count is met accept remaining
95c3e38098bf Untabify .m scripts
Rik <code@nomad.inbox5.com>
parents: 9590
diff changeset
364 ## sub-interval and exit
95c3e38098bf Untabify .m scripts
Rik <code@nomad.inbox5.com>
parents: 9590
diff changeset
365 if (rows (subs) > maxint)
95c3e38098bf Untabify .m scripts
Rik <code@nomad.inbox5.com>
parents: 9590
diff changeset
366 warning ("quadgk: maximum interval count (%d) met", maxint);
95c3e38098bf Untabify .m scripts
Rik <code@nomad.inbox5.com>
parents: 9590
diff changeset
367 q += sum (q_subs);
95c3e38098bf Untabify .m scripts
Rik <code@nomad.inbox5.com>
parents: 9590
diff changeset
368 err += sum (q_errs);
95c3e38098bf Untabify .m scripts
Rik <code@nomad.inbox5.com>
parents: 9590
diff changeset
369 break;
95c3e38098bf Untabify .m scripts
Rik <code@nomad.inbox5.com>
parents: 9590
diff changeset
370 endif
7771
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
371
10549
95c3e38098bf Untabify .m scripts
Rik <code@nomad.inbox5.com>
parents: 9590
diff changeset
372 ## Evaluation of the integrand on the remaining sub-intervals
95c3e38098bf Untabify .m scripts
Rik <code@nomad.inbox5.com>
parents: 9590
diff changeset
373 [q_subs, q_errs] = __quadgk_eval__ (f, subs);
7771
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
374 endwhile
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
375
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
376 if (err > max (abstol, reltol * abs(q)))
10549
95c3e38098bf Untabify .m scripts
Rik <code@nomad.inbox5.com>
parents: 9590
diff changeset
377 warning ("quadgk: Error tolerance not met. Estimated error %g", err);
7771
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
378 endif
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
379 unwind_protect_cleanup
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
380 if (strcmp (warn_state.state, "on"))
10549
95c3e38098bf Untabify .m scripts
Rik <code@nomad.inbox5.com>
parents: 9590
diff changeset
381 warning ("on", "Octave:divide-by-zero");
7771
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
382 endif
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
383 end_unwind_protect
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
384 endif
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
385 endfunction
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
386
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
387 function [q, err] = __quadgk_eval__ (f, subs)
9067
8970b4b10e9f Cleanup documentation for quad.texi and diffeq.texi
Rik <rdrider0-list@yahoo.com>
parents: 9051
diff changeset
388 ## A (15,7) point pair of Gauss-Konrod quadrature rules. The abscissa
7771
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
389 ## and weights are copied directly from dqk15w.f from quadpack
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
390
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
391 persistent abscissa = [-0.9914553711208126e+00, -0.9491079123427585e+00, ...
10549
95c3e38098bf Untabify .m scripts
Rik <code@nomad.inbox5.com>
parents: 9590
diff changeset
392 -0.8648644233597691e+00, -0.7415311855993944e+00, ...
95c3e38098bf Untabify .m scripts
Rik <code@nomad.inbox5.com>
parents: 9590
diff changeset
393 -0.5860872354676911e+00, -0.4058451513773972e+00, ...
95c3e38098bf Untabify .m scripts
Rik <code@nomad.inbox5.com>
parents: 9590
diff changeset
394 -0.2077849550078985e+00, 0.0000000000000000e+00, ...
95c3e38098bf Untabify .m scripts
Rik <code@nomad.inbox5.com>
parents: 9590
diff changeset
395 0.2077849550078985e+00, 0.4058451513773972e+00, ...
95c3e38098bf Untabify .m scripts
Rik <code@nomad.inbox5.com>
parents: 9590
diff changeset
396 0.5860872354676911e+00, 0.7415311855993944e+00, ...
95c3e38098bf Untabify .m scripts
Rik <code@nomad.inbox5.com>
parents: 9590
diff changeset
397 0.8648644233597691e+00, 0.9491079123427585e+00, ...
95c3e38098bf Untabify .m scripts
Rik <code@nomad.inbox5.com>
parents: 9590
diff changeset
398 0.9914553711208126e+00];
7771
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
399
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
400 persistent weights15 = ...
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
401 diag ([0.2293532201052922e-01, 0.6309209262997855e-01, ...
10549
95c3e38098bf Untabify .m scripts
Rik <code@nomad.inbox5.com>
parents: 9590
diff changeset
402 0.1047900103222502e+00, 0.1406532597155259e+00, ...
95c3e38098bf Untabify .m scripts
Rik <code@nomad.inbox5.com>
parents: 9590
diff changeset
403 0.1690047266392679e+00, 0.1903505780647854e+00, ...
95c3e38098bf Untabify .m scripts
Rik <code@nomad.inbox5.com>
parents: 9590
diff changeset
404 0.2044329400752989e+00, 0.2094821410847278e+00, ...
95c3e38098bf Untabify .m scripts
Rik <code@nomad.inbox5.com>
parents: 9590
diff changeset
405 0.2044329400752989e+00, 0.1903505780647854e+00, ...
95c3e38098bf Untabify .m scripts
Rik <code@nomad.inbox5.com>
parents: 9590
diff changeset
406 0.1690047266392679e+00, 0.1406532597155259e+00, ...
95c3e38098bf Untabify .m scripts
Rik <code@nomad.inbox5.com>
parents: 9590
diff changeset
407 0.1047900103222502e+00, 0.6309209262997855e-01, ...
95c3e38098bf Untabify .m scripts
Rik <code@nomad.inbox5.com>
parents: 9590
diff changeset
408 0.2293532201052922e-01]);
7771
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
409
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
410 persistent weights7 = ...
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
411 diag ([0.1294849661688697e+00, 0.2797053914892767e+00, ...
10549
95c3e38098bf Untabify .m scripts
Rik <code@nomad.inbox5.com>
parents: 9590
diff changeset
412 0.3818300505051889e+00, 0.4179591836734694e+00, ...
95c3e38098bf Untabify .m scripts
Rik <code@nomad.inbox5.com>
parents: 9590
diff changeset
413 0.3818300505051889e+00, 0.2797053914892767e+00, ...
95c3e38098bf Untabify .m scripts
Rik <code@nomad.inbox5.com>
parents: 9590
diff changeset
414 0.1294849661688697e+00]);
7771
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
415
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
416 halfwidth = diff (subs, [], 2) ./ 2;
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
417 center = sum (subs, 2) ./ 2;;
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
418 x = bsxfun (@plus, halfwidth * abscissa, center);
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
419 y = reshape (f (x(:)), size(x));
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
420
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
421 ## This is faster than using bsxfun as the * operator can use a
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
422 ## single BLAS call, rather than rows(sub) calls to the @times
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
423 ## function.
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
424 q = sum (y * weights15, 2) .* halfwidth;
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
425 err = abs (sum (y(:,2:2:end) * weights7, 2) .* halfwidth - q);
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
426 endfunction
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
427
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
428 function t = __quadgk_finite_waypoint__ (x, a, b)
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
429 c = (-4 .* x + 2.* (b + a)) ./ (b - a);
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
430 k = ((sqrt(c .^ 2 - 4) + c) ./ 2) .^ (1/3);
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
431 t = real ((sqrt(3) .* 1i * (1 - k .^ 2) - (1 + k .^ 2)) ./ 2 ./ k);
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
432 endfunction
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
433
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
434 %error (quadgk (@sin))
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
435 %error (quadgk (@sin, -pi))
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
436 %error (quadgk (@sin, -pi, pi, 'DummyArg'))
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
437
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
438 %!assert (quadgk(@sin,-pi,pi), 0, 1e-6)
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
439 %!assert (quadgk(inline('sin'),-pi,pi), 0, 1e-6)
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
440 %!assert (quadgk('sin',-pi,pi), 0, 1e-6)
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
441 %!assert (quadgk(@sin,-pi,pi,'waypoints', 0, 'MaxIntervalCount', 100, 'reltol', 1e-3, 'abstol', 1e-6, 'trace', false), 0, 1e-6)
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
442 %!assert (quadgk(@sin,-pi,pi,1e-6,false), 0, 1e-6)
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
443
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
444 %!assert (quadgk(@sin,-pi,0), -2, 1e-6)
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
445 %!assert (quadgk(@sin,0,pi), 2, 1e-6)
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
446 %!assert (quadgk(@(x) 1./sqrt(x), 0, 1), 2, 1e-6)
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
447 %!assert (quadgk (@(x) abs (1 - x.^2), 0, 2, 'Waypoints', 1), 2, 1e-6)
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
448 %!assert (quadgk(@(x) 1./(sqrt(x).*(x+1)), 0, Inf), pi, 1e-6)
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
449 %!assert (quadgk (@(z) log (z), 1+1i, 1+1i, 'WayPoints', [1-1i, -1,-1i, -1+1i]), -pi * 1i, 1e-6)
9323
8a3afa8ebb06 Fixes for quadgk for doubly infinite interval
David Bateman <dbateman@free.fr>
parents: 9067
diff changeset
450
10846
a4f482e66b65 Grammarcheck more of the documentation.
Rik <octave@nomad.inbox5.com>
parents: 10821
diff changeset
451 %!assert (quadgk (@(x) exp(-x .^ 2), -Inf, Inf), sqrt(pi), 1e-6)