Mercurial > hg > octave-lyh

annotate scripts/general/quadgk.m @ 8522:d65b33e55d40

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
nargchk.m: improve compatibility; new tests
author Peter L. Sondergaard <peter@sonderport.dk>
date Thu, 15 Jan 2009 01:06:06 -0500
parents e2a179415bac
children e07e93c04080
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
7771
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
1 ## Copyright (C) 2008 David Bateman
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{})
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
23 ## Numerically evaluate integral using adaptive Guass-Konrod quadrature.
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,
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
25 ## @cite{"Vectorized adaptive quadrature in MATLAB", Journal of
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,
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
27 ## Feb 2008} where all function evalutions at an iteration are
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
28 ## calculated with a single call to @var{f}. Therefore the function
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
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
31 ## representing the function evalutaions at the given values of @var{x}.
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
7771
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
36 ## or infinite and contain weak end singularities. Variable
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
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
38 ## the singularities. For example
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
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
50 ## compatible with @code{quadv}. Equally the user can request that
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
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
52 ## is logicallly true.
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
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
55 ## pairs @code{@var{prop}, @var{val}}. Valid properties are
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
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
59 ## Defines the absolute error tolerance for the quadrature. The default
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
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
63 ## Defines the relative error tolerance for the quadrature. The default
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
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
68 ## the quadrature into 10 intervals. Sub-intervals that have an
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
69 ## unacceptable error are sub-divided and re-evaluated. If the number of
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
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
72 ## returned. The property 'MaxIntervalCount' can be used to alter the
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
e2a179415bac doc fixes
John W. Eaton <jwe@octave.org>
parents: 8507
diff changeset
78 ## @code{"WayPoints"} property. This forces the ends of a sub-interval
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
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
97 ## continuous path defined by the above. In this case the integral is
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
98 ## assuemd to have no edge singularities. For example
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
8516
e2a179415bac doc fixes
John W. Eaton <jwe@octave.org>
parents: 8507
diff changeset
101 ## quadgk (@@(z) log (z), 1+1i, 1+1i, "WayPoints",
e2a179415bac doc fixes
John W. Eaton <jwe@octave.org>
parents: 8507
diff changeset
102 ## [1-1i, -1,-1i, -1+1i])
7771
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
103 ## @end example
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
104 ##
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
105 ## @noindent
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
106 ## 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
107 ## 1-1i, -1-1i, -1+1i]}
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
108 ##
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
109 ## 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
110 ## 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
111 ## @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
112 ##
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
113 ## @seealso{triplequad, dblquad, quad, quadl, quadv, trapz}
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
114 ## @end deftypefn
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
115
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
116 function [q, err] = quadgk (f, a, b, varargin)
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
117 if (nargin < 3)
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
118 print_usage ();
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
119 endif
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
120
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
121 if (b < a)
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
122 [q, err] = quadgk (f, b, a, varargin{:});
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
123 q = -q;
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
124 else
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
125 abstol = 1e-10;
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
126 reltol = 1e-5;
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
127 waypoints = [];
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
128 maxint = 650;
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
129 trace = false;
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
130
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
131 if (nargin > 3)
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
132 if (! ischar (varargin{1}))
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
133 if (!isempty (varargin{1}))
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
134 abstol = varargin{1};
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
135 reltol = 0;
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
136 endif
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
137 if (nargin > 4)
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
138 trace = varargin{2};
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
139 endif
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
140 if (nargin > 5)
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
141 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
142 endif
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
143 else
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
144 idx = 1;
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
145 while (idx < nargin - 3)
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
146 if (ischar (varargin{idx}))
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
147 str = varargin{idx++};
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
148 if (strcmpi (str, "reltol"))
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
149 reltol = varargin{idx++};
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
150 elseif (strcmpi (str, "abstol"))
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
151 abstol = varargin{idx++};
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
152 elseif (strcmpi (str, "waypoints"))
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
153 waypoints = varargin{idx++} (:);
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
154 if (isreal(waypoints))
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
155 waypoints (waypoints < a | waypoints > b) = [];
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
156 endif
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
157 elseif (strcmpi (str, "maxintervalcount"))
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
158 maxint = varargin{idx++};
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
159 elseif (strcmpi (str, "trace"))
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
160 trace = varargin{idx++};
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
161 else
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
162 error ("quadgk: unknown property %s", str);
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
163 endif
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
164 else
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
165 error ("quadgk: expecting property to be a string");
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
166 endif
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
167 endwhile
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
168 if (idx != nargin - 2)
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
169 error ("quadgk: expecting properties in pairs")
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
170 endif
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
171 endif
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
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
174 ## 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
175 if (ischar (f))
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
176 f = @(x) feval (f, x);
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
177 endif
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
178
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
179 ## 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
180 ## 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
181 ## contour integrals
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
182 if (iscomplex (a) || iscomplex (b) || iscomplex(waypoints))
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
183 ## contour integral, no transform
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
184 subs = [a; waypoints; b];
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
185 h = b - a;
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
186 trans = @(t) t;
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
187 elseif (isinf (a) && isinf(b))
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
188 ## Standard Infinite to finite integral transformation.
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
189 ## \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
190 ## where
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
191 ## g(t) = t / (1 - t^2)
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
192 ## g'(t) = 1 / (1 + t^2) ^ 2
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
193 ## waypoint transform is then
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
194 ## 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
195 if (!isempty (waypoints))
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
196 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
197 subs = [-1; trans(waypoints); 1];
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
198 else
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
199 subs = linspace (-1, 1, 11)';
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
200 endif
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
201 h = 2;
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
202 trans = @(t) t ./ (1 - t.^2);
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
203 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
204 elseif (isinf(a))
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
205 ## 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
206 ## 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
207 ## a finite interval. The singularity weakening transform is
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
208 ## \int_{-\infinity}^b f(x) dx =
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
209 ## - \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
210 ## (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
211 ## \int_{-\infinity}^0 f(b - t^2) 2 t dt =
8507
cadc73247d65 style fixes
John W. Eaton <jwe@octave.org>
parents: 8122
diff changeset
212 ## \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
213 ## where
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
214 ## g(s) = s / (1 + s)
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
215 ## g'(s) = 1 / (1 + s) ^ 2
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
216 ## waypoint transform is then
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
217 ## t = sqrt (b - x)
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
218 ## s = - t / (t + 1)
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
219 if (!isempty (waypoints))
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
220 tmp = sqrt (b - waypoints);
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
221 trans = @(x) - x ./ (x + 1);
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
222 subs = [0; trans(tmp); 1];
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
223 else
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
224 subs = linspace (0, 1, 11)';
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
225 endif
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
226 h = 1;
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
227 trans = @(t) b - (t ./ (1 + t)).^2;
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
228 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
229 elseif (isinf(b))
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
230 ## 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
231 ## 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
232 ## a finite interval. The singularity weakening transform is
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
233 ## \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
234 ## and the finite interval transform is
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
235 ## \int_0^\infinity f(a + t^2) 2 t dt =
8507
cadc73247d65 style fixes
John W. Eaton <jwe@octave.org>
parents: 8122
diff changeset
236 ## \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
237 ## where
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
238 ## g(s) = s / (1 - s)
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
239 ## g'(s) = 1 / (1 - s) ^ 2
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
240 ## waypoint transform is then
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
241 ## t = sqrt (x - a)
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
242 ## s = t / (t + 1)
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
243 if (!isempty (waypoints))
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
244 tmp = sqrt (waypoints - a);
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
245 trans = @(x) x ./ (x + 1);
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
246 subs = [0; trans(tmp); 1];
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
247 else
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
248 subs = linspace (0, 1, 11)';
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
249 endif
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
250 h = 1;
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
251 trans = @(t) a + (t ./ (1 - t)).^2;
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
252 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
253 else
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
254 ## Davis, Rabinowitz, "Methods of Numerical Integration" p441 2ed.
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
255 ## 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
256 ## 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
257 ## g'(t) = ((b-a)/4) * (3 - 3t^2);
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
258 ## 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
259 ## 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
260 ## 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
261 ## ignored. The third is
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
262 ## c = (-4*x + 2*(b+a)) / (b-a);
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
263 ## 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
264 ## 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
265 if (! isempty (waypoints))
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
266 trans = @__quadgk_finite_waypoint__;
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
267 subs = [-1; trans(waypoints, a, b); 1];
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
268 else
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
269 subs = linspace(-1, 1, 11)';
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
270 endif
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
271 h = 2;
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
272 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
273 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
274 3 .* (b - a) ./ 4 .* (1 - t.^2);
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
275 endif
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
276
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
277 ## Split interval into at least 10 sub-interval with a 15 point
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
278 ## Guass-Kronrod rule giving a minimum of 150 function evaluations
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
279 while (length (subs) < 11)
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
280 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
281 endwhile
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
282 subs = [subs(1:end-1), subs(2:end)];
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 warn_state = warning ("query", "Octave:divide-by-zero");
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 unwind_protect
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
287 ## Singularity will cause divide by zero warnings
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
288 warning ("off", "Octave:divide-by-zero")
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
289
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
290 ## 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
291 [q_subs, q_errs] = __quadgk_eval__ (f, subs);
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
292 q0 = sum (q_subs);
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
293 err0 = sum (q_errs);
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
294
7795
df9519e9990c Handle single precision eps values
David Bateman <dbateman@free.fr>
parents: 7771
diff changeset
295 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
296 myeps = eps ("single");
df9519e9990c Handle single precision eps values
David Bateman <dbateman@free.fr>
parents: 7771
diff changeset
297 else
df9519e9990c Handle single precision eps values
David Bateman <dbateman@free.fr>
parents: 7771
diff changeset
298 myeps = eps;
df9519e9990c Handle single precision eps values
David Bateman <dbateman@free.fr>
parents: 7771
diff changeset
299 endif
df9519e9990c Handle single precision eps values
David Bateman <dbateman@free.fr>
parents: 7771
diff changeset
300
7771
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
301 first = true;
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
302 while (true)
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
303 ## 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
304 ## 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
305 ## value for this test. Shampine suggests 100*eps
7795
df9519e9990c Handle single precision eps values
David Bateman <dbateman@free.fr>
parents: 7771
diff changeset
306 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
307 q = q0;
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
308 err = err0;
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
309 break;
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
310 endif
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
311
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
312 ## 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
313 if (any (! isfinite (q_subs)))
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
314 warning ("quadgk: non finite integrand encountered");
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
315 q = q0;
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
316 err = err0;
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
317 break;
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
318 endif
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
319
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
320 tol = max (abstol, reltol .* abs (q0));
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 ## If the global error estimate is meet exit
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
323 if (err0 < tol)
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
324 q = q0;
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
325 err = err0;
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
326 break;
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
327 endif
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
328
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
329 ## 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
330 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
331 if (first)
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
332 q = sum (q_subs (idx));
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
333 err = sum (q_errs(idx));
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
334 first = false;
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
335 else
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
336 q0 = q + sum (q_subs);
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
337 err0 = err + sum (q_errs);
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
338 q += sum (q_subs (idx));
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
339 err += sum (q_errs(idx));
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
340 endif
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
341 subs(idx,:) = [];
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
342
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
343 ## If no remaining sub-intervals exit
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
344 if (rows (subs) == 0)
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
345 break;
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
346 endif
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
347
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
348 if (trace)
8122
99602635172a Trivial patch to quadgk for absence of trace function count
David Bateman <dbateman@free.fr>
parents: 7989
diff changeset
349 disp([rows(subs), err, q0]);
7771
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
350 endif
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
351
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
352 ## Split remaining sub-intervals in two
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
353 mid = (subs(:,2) + subs(:,1)) ./ 2;
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
354 subs = [subs(:,1), mid; mid, subs(:,2)];
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
355
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
356 ## 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
357 ## sub-interval and exit
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
358 if (rows (subs) > maxint)
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
359 warning ("quadgk: maximum interval count (%d) met", maxint);
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
360 q += sum (q_subs);
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
361 err += sum (q_errs);
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
362 break;
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
363 endif
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
364
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
365 ## 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
366 [q_subs, q_errs] = __quadgk_eval__ (f, subs);
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
367 endwhile
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
368
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
369 if (err > max (abstol, reltol * abs(q)))
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
370 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
371 endif
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
372 unwind_protect_cleanup
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
373 if (strcmp (warn_state.state, "on"))
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
374 warning ("on", "Octave:divide-by-zero")
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
375 endif
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
376 end_unwind_protect
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 endfunction
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
379
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
380 function [q, err] = __quadgk_eval__ (f, subs)
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
381 ## A (15,7) point pair of Guass-Konrod quadrature rules. The abscissa
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
382 ## 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
383
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
384 persistent abscissa = [-0.9914553711208126e+00, -0.9491079123427585e+00, ...
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
385 -0.8648644233597691e+00, -0.7415311855993944e+00, ...
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
386 -0.5860872354676911e+00, -0.4058451513773972e+00, ...
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
387 -0.2077849550078985e+00, 0.0000000000000000e+00, ...
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
388 0.2077849550078985e+00, 0.4058451513773972e+00, ...
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
389 0.5860872354676911e+00, 0.7415311855993944e+00, ...
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
390 0.8648644233597691e+00, 0.9491079123427585e+00, ...
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
391 0.9914553711208126e+00];
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
392
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
393 persistent weights15 = ...
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
394 diag ([0.2293532201052922e-01, 0.6309209262997855e-01, ...
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
395 0.1047900103222502e+00, 0.1406532597155259e+00, ...
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
396 0.1690047266392679e+00, 0.1903505780647854e+00, ...
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
397 0.2044329400752989e+00, 0.2094821410847278e+00, ...
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
398 0.2044329400752989e+00, 0.1903505780647854e+00, ...
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
399 0.1690047266392679e+00, 0.1406532597155259e+00, ...
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
400 0.1047900103222502e+00, 0.6309209262997855e-01, ...
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
401 0.2293532201052922e-01]);
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
402
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
403 persistent weights7 = ...
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
404 diag ([0.1294849661688697e+00, 0.2797053914892767e+00, ...
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
405 0.3818300505051889e+00, 0.4179591836734694e+00, ...
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
406 0.3818300505051889e+00, 0.2797053914892767e+00, ...
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
407 0.1294849661688697e+00]);
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
408
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
409 halfwidth = diff (subs, [], 2) ./ 2;
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
410 center = sum (subs, 2) ./ 2;;
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
411 x = bsxfun (@plus, halfwidth * abscissa, center);
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
412 y = reshape (f (x(:)), size(x));
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
413
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
414 ## 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
415 ## 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
416 ## function.
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
417 q = sum (y * weights15, 2) .* halfwidth;
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
418 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
419 endfunction
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 function t = __quadgk_finite_waypoint__ (x, a, b)
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
422 c = (-4 .* x + 2.* (b + a)) ./ (b - a);
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
423 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
424 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
425 endfunction
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
426
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
427 %error (quadgk (@sin))
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
428 %error (quadgk (@sin, -pi))
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
429 %error (quadgk (@sin, -pi, pi, 'DummyArg'))
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
430
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
431 %!assert (quadgk(@sin,-pi,pi), 0, 1e-6)
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
432 %!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
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(@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
435 %!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
436
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
437 %!assert (quadgk(@sin,-pi,0), -2, 1e-6)
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
438 %!assert (quadgk(@sin,0,pi), 2, 1e-6)
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
439 %!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
440 %!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
441 %!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
442 %!assert (quadgk (@(z) log (z), 1+1i, 1+1i, 'WayPoints', [1-1i, -1,-1i, -1+1i]), -pi * 1i, 1e-6)