Mercurial > octave-nkf

annotate scripts/general/quadgk.m @ 9067:8970b4b10e9f

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Cleanup documentation for quad.texi and diffeq.texi Grammarcheck input .txi files Spellcheck .texi files
author Rik <rdrider0-list@yahoo.com>
date Sat, 28 Mar 2009 22:18:51 -0700
parents 1bf0ce0930be
children 8a3afa8ebb06
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
8920
eb63fbe60fab update copyright notices
John W. Eaton <jwe@octave.org>
parents: 8664
diff changeset
1 ## Copyright (C) 2008, 2009 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 -*-
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
20 ## @deftypefn {Function File} {} quadgk (@var{f}, @var{a}, @var{b}, @var{abstol}, @var{trace})
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
9051
1bf0ce0930be Grammar check TexInfo in all .m files
Rik <rdrider0-list@yahoo.com>
parents: 8920
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
e2a179415bac doc fixes
John W. Eaton <jwe@octave.org>
parents: 8507
diff changeset
80 ## significantly improved estimated of the error in the integral, faster
e2a179415bac doc fixes
John W. Eaton <jwe@octave.org>
parents: 8507
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
9067
8970b4b10e9f Cleanup documentation for quad.texi and diffeq.texi
Rik <rdrider0-list@yahoo.com>
parents: 9051
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}))
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
135 if (!isempty (varargin{1}))
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
136 abstol = varargin{1};
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
137 reltol = 0;
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
138 endif
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
139 if (nargin > 4)
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
140 trace = varargin{2};
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
141 endif
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
142 if (nargin > 5)
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
143 error ("quadgk: can not pass additional arguments to user function");
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
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
146 idx = 1;
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
147 while (idx < nargin - 3)
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
148 if (ischar (varargin{idx}))
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
149 str = varargin{idx++};
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
150 if (strcmpi (str, "reltol"))
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
151 reltol = varargin{idx++};
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
152 elseif (strcmpi (str, "abstol"))
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
153 abstol = varargin{idx++};
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
154 elseif (strcmpi (str, "waypoints"))
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
155 waypoints = varargin{idx++} (:);
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
156 if (isreal(waypoints))
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
157 waypoints (waypoints < a | waypoints > b) = [];
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
158 endif
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
159 elseif (strcmpi (str, "maxintervalcount"))
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
160 maxint = varargin{idx++};
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
161 elseif (strcmpi (str, "trace"))
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
162 trace = varargin{idx++};
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
163 else
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
164 error ("quadgk: unknown property %s", str);
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
165 endif
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
166 else
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
167 error ("quadgk: expecting property to be a string");
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
168 endif
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
169 endwhile
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
170 if (idx != nargin - 2)
8664
e07e93c04080 style fixes
John W. Eaton <jwe@octave.org>
parents: 8516
diff changeset
171 error ("quadgk: expecting properties in pairs");
7771
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
172 endif
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];
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
187 h = b - a;
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
188 trans = @(t) t;
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
189 elseif (isinf (a) && isinf(b))
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
190 ## Standard Infinite to finite integral transformation.
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
191 ## \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
192 ## where
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
193 ## g(t) = t / (1 - t^2)
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
194 ## g'(t) = 1 / (1 + t^2) ^ 2
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
195 ## waypoint transform is then
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
196 ## t = (-1 + sqrt(1 + 4 * g(t) .^ 2)) ./ (2 * g(t))
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
197 if (!isempty (waypoints))
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
198 trans = @(x) (-1 + sqrt(1 + 4 * x .^ 2)) ./ (2 * x);
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
199 subs = [-1; trans(waypoints); 1];
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
200 else
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
201 subs = linspace (-1, 1, 11)';
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
202 endif
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
203 h = 2;
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
204 trans = @(t) t ./ (1 - t.^2);
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
205 f = @(t) f (t ./ (1 + t .^ 2)) ./ ((1 + t .^ 2) .^ 2);
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
206 elseif (isinf(a))
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
207 ## 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
208 ## 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
209 ## a finite interval. The singularity weakening transform is
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
210 ## \int_{-\infinity}^b f(x) dx =
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
211 ## - \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
212 ## (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
213 ## \int_{-\infinity}^0 f(b - t^2) 2 t dt =
8507
cadc73247d65 style fixes
John W. Eaton <jwe@octave.org>
parents: 8122
diff changeset
214 ## \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
215 ## where
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
216 ## g(s) = s / (1 + s)
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
217 ## g'(s) = 1 / (1 + s) ^ 2
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
218 ## waypoint transform is then
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
219 ## t = sqrt (b - x)
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
220 ## s = - t / (t + 1)
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
221 if (!isempty (waypoints))
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
222 tmp = sqrt (b - waypoints);
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
223 trans = @(x) - x ./ (x + 1);
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
224 subs = [0; trans(tmp); 1];
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
225 else
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
226 subs = linspace (0, 1, 11)';
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
227 endif
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
228 h = 1;
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
229 trans = @(t) b - (t ./ (1 + t)).^2;
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
230 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
231 elseif (isinf(b))
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
232 ## 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
233 ## 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
234 ## a finite interval. The singularity weakening transform is
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
235 ## \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
236 ## and the finite interval transform is
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
237 ## \int_0^\infinity f(a + t^2) 2 t dt =
8507
cadc73247d65 style fixes
John W. Eaton <jwe@octave.org>
parents: 8122
diff changeset
238 ## \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
239 ## where
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
240 ## g(s) = s / (1 - s)
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
241 ## g'(s) = 1 / (1 - s) ^ 2
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
242 ## waypoint transform is then
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
243 ## t = sqrt (x - a)
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
244 ## s = t / (t + 1)
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
245 if (!isempty (waypoints))
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
246 tmp = sqrt (waypoints - a);
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
247 trans = @(x) x ./ (x + 1);
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
248 subs = [0; trans(tmp); 1];
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
249 else
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
250 subs = linspace (0, 1, 11)';
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
251 endif
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
252 h = 1;
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
253 trans = @(t) a + (t ./ (1 - t)).^2;
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
254 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
255 else
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
256 ## Davis, Rabinowitz, "Methods of Numerical Integration" p441 2ed.
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
257 ## 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
258 ## 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
259 ## g'(t) = ((b-a)/4) * (3 - 3t^2);
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
260 ## 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
261 ## 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
262 ## 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
263 ## ignored. The third is
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
264 ## c = (-4*x + 2*(b+a)) / (b-a);
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
265 ## 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
266 ## 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
267 if (! isempty (waypoints))
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
268 trans = @__quadgk_finite_waypoint__;
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
269 subs = [-1; trans(waypoints, a, b); 1];
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
270 else
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
271 subs = linspace(-1, 1, 11)';
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
272 endif
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
273 h = 2;
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
274 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
275 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
276 3 .* (b - a) ./ 4 .* (1 - t.^2);
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
277 endif
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
278
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
279 ## 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
280 ## 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
281 while (length (subs) < 11)
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
282 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
283 endwhile
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
284 subs = [subs(1:end-1), subs(2:end)];
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
285
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
286 warn_state = warning ("query", "Octave:divide-by-zero");
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
287
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
288 unwind_protect
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
289 ## Singularity will cause divide by zero warnings
8664
e07e93c04080 style fixes
John W. Eaton <jwe@octave.org>
parents: 8516
diff changeset
290 warning ("off", "Octave:divide-by-zero");
7771
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
291
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
292 ## 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
293 [q_subs, q_errs] = __quadgk_eval__ (f, subs);
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
294 q0 = sum (q_subs);
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
295 err0 = sum (q_errs);
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
296
7795
df9519e9990c Handle single precision eps values
David Bateman <dbateman@free.fr>
parents: 7771
diff changeset
297 if (isa (a, "single") || isa (b, "single") || isa (waypoints, "single"))
df9519e9990c Handle single precision eps values
David Bateman <dbateman@free.fr>
parents: 7771
diff changeset
298 myeps = eps ("single");
df9519e9990c Handle single precision eps values
David Bateman <dbateman@free.fr>
parents: 7771
diff changeset
299 else
df9519e9990c Handle single precision eps values
David Bateman <dbateman@free.fr>
parents: 7771
diff changeset
300 myeps = eps;
df9519e9990c Handle single precision eps values
David Bateman <dbateman@free.fr>
parents: 7771
diff changeset
301 endif
df9519e9990c Handle single precision eps values
David Bateman <dbateman@free.fr>
parents: 7771
diff changeset
302
7771
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
303 first = true;
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
304 while (true)
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
305 ## Check for sub-intervals that are too small. Test must be
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
306 ## performed in untransformed sub-intervals. What is a good
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
307 ## value for this test. Shampine suggests 100*eps
7795
df9519e9990c Handle single precision eps values
David Bateman <dbateman@free.fr>
parents: 7771
diff changeset
308 if (any (diff (trans (subs), [], 2) / (b - a) < 100 * myeps))
7771
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
309 q = q0;
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
310 err = err0;
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
311 break;
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
312 endif
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
313
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
314 ## Quit if any evaluations are not finite (Inf or NaN)
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
315 if (any (! isfinite (q_subs)))
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
316 warning ("quadgk: non finite integrand encountered");
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
317 q = q0;
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
318 err = err0;
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
319 break;
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
320 endif
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
321
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
322 tol = max (abstol, reltol .* abs (q0));
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
323
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
324 ## If the global error estimate is meet exit
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
325 if (err0 < tol)
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
326 q = q0;
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
327 err = err0;
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
328 break;
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
329 endif
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
330
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
331 ## Accept the sub-intervals that meet the convergence criteria
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
332 idx = find (abs (q_errs) < tol .* diff (subs, [], 2) ./ h);
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
333 if (first)
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
334 q = sum (q_subs (idx));
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
335 err = sum (q_errs(idx));
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
336 first = false;
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
337 else
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
338 q0 = q + sum (q_subs);
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
339 err0 = err + sum (q_errs);
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
340 q += sum (q_subs (idx));
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
341 err += sum (q_errs(idx));
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
342 endif
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
343 subs(idx,:) = [];
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
344
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
345 ## If no remaining sub-intervals exit
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
346 if (rows (subs) == 0)
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
347 break;
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
348 endif
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
349
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
350 if (trace)
8122
99602635172a Trivial patch to quadgk for absence of trace function count
David Bateman <dbateman@free.fr>
parents: 7989
diff changeset
351 disp([rows(subs), err, q0]);
7771
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
352 endif
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
353
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
354 ## Split remaining sub-intervals in two
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
355 mid = (subs(:,2) + subs(:,1)) ./ 2;
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
356 subs = [subs(:,1), mid; mid, subs(:,2)];
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
357
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
358 ## If the maximum sub-interval count is met accept remaining
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
359 ## sub-interval and exit
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
360 if (rows (subs) > maxint)
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
361 warning ("quadgk: maximum interval count (%d) met", maxint);
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
362 q += sum (q_subs);
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
363 err += sum (q_errs);
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
364 break;
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
365 endif
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
366
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
367 ## Evaluation of the integrand on the remaining sub-intervals
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
368 [q_subs, q_errs] = __quadgk_eval__ (f, subs);
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
369 endwhile
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
370
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
371 if (err > max (abstol, reltol * abs(q)))
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
372 warning ("quadgk: Error tolerance not met. Estimated error %g", err);
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
373 endif
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
374 unwind_protect_cleanup
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
375 if (strcmp (warn_state.state, "on"))
8664
e07e93c04080 style fixes
John W. Eaton <jwe@octave.org>
parents: 8516
diff changeset
376 warning ("on", "Octave:divide-by-zero");
7771
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
377 endif
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
378 end_unwind_protect
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
379 endif
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
380 endfunction
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
381
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
382 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
383 ## 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
384 ## 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
385
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
386 persistent abscissa = [-0.9914553711208126e+00, -0.9491079123427585e+00, ...
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
387 -0.8648644233597691e+00, -0.7415311855993944e+00, ...
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
388 -0.5860872354676911e+00, -0.4058451513773972e+00, ...
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
389 -0.2077849550078985e+00, 0.0000000000000000e+00, ...
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
390 0.2077849550078985e+00, 0.4058451513773972e+00, ...
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
391 0.5860872354676911e+00, 0.7415311855993944e+00, ...
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
392 0.8648644233597691e+00, 0.9491079123427585e+00, ...
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
393 0.9914553711208126e+00];
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
394
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
395 persistent weights15 = ...
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
396 diag ([0.2293532201052922e-01, 0.6309209262997855e-01, ...
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
397 0.1047900103222502e+00, 0.1406532597155259e+00, ...
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
398 0.1690047266392679e+00, 0.1903505780647854e+00, ...
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
399 0.2044329400752989e+00, 0.2094821410847278e+00, ...
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
400 0.2044329400752989e+00, 0.1903505780647854e+00, ...
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
401 0.1690047266392679e+00, 0.1406532597155259e+00, ...
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
402 0.1047900103222502e+00, 0.6309209262997855e-01, ...
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
403 0.2293532201052922e-01]);
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
404
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
405 persistent weights7 = ...
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
406 diag ([0.1294849661688697e+00, 0.2797053914892767e+00, ...
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
407 0.3818300505051889e+00, 0.4179591836734694e+00, ...
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
408 0.3818300505051889e+00, 0.2797053914892767e+00, ...
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
409 0.1294849661688697e+00]);
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
410
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
411 halfwidth = diff (subs, [], 2) ./ 2;
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
412 center = sum (subs, 2) ./ 2;;
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
413 x = bsxfun (@plus, halfwidth * abscissa, center);
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
414 y = reshape (f (x(:)), size(x));
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 ## 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
417 ## 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
418 ## function.
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
419 q = sum (y * weights15, 2) .* halfwidth;
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
420 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
421 endfunction
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
422
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
423 function t = __quadgk_finite_waypoint__ (x, a, b)
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
424 c = (-4 .* x + 2.* (b + a)) ./ (b - a);
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
425 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
426 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
427 endfunction
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
428
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
429 %error (quadgk (@sin))
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
430 %error (quadgk (@sin, -pi))
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
431 %error (quadgk (@sin, -pi, pi, 'DummyArg'))
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
432
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
433 %!assert (quadgk(@sin,-pi,pi), 0, 1e-6)
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
434 %!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
435 %!assert (quadgk('sin',-pi,pi), 0, 1e-6)
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
436 %!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
437 %!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
438
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
439 %!assert (quadgk(@sin,-pi,0), -2, 1e-6)
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
440 %!assert (quadgk(@sin,0,pi), 2, 1e-6)
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
441 %!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
442 %!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
443 %!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
444 %!assert (quadgk (@(z) log (z), 1+1i, 1+1i, 'WayPoints', [1-1i, -1,-1i, -1+1i]), -pi * 1i, 1e-6)