Mercurial > hg > octave-lyh

annotate scripts/general/quadgk.m @ 7795:df9519e9990c

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Handle single precision eps values
author David Bateman <dbateman@free.fr>
date Mon, 12 May 2008 22:57:11 +0200
parents 680631e787aa
children 23c248d415b5
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 ##
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
35 ## The bounds of the quadrature @code{[@var{a}, @var{b}] can be finite
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
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
77 ## function to integrate, then these can be flagged with teh 'WayPoints'
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
78 ## property. This forces the ends of a sub-interval to fall on the
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
79 ## breakpoints of the function and can result in significantly improved
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
80 ## estimated of the error in the integral, faster computation or both.
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
81 ## For example
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
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
101 ## quadgk (@@ (z) log (z), 1+1i, 1+1i, 'WayPoints', [1-1i, -1,-1i, -1+1i])
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
102 ## @end example
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
103 ##
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
104 ## @noindent
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
105 ## 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
106 ## 1-1i, -1-1i, -1+1i]}
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
107 ##
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
108 ## 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
109 ## 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
110 ## @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
111 ##
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
112 ## @seealso{triplequad, dblquad, quad, quadl, quadv, trapz}
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
113 ## @end deftypefn
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 function [q, err] = quadgk (f, a, b, varargin)
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
116 if (nargin < 3)
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
117 print_usage ();
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
118 endif
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
119
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
120 if (b < a)
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
121 [q, err] = quadgk (f, b, a, varargin{:});
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
122 q = -q;
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
123 else
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
124 abstol = 1e-10;
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
125 reltol = 1e-5;
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
126 waypoints = [];
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
127 maxint = 650;
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
128 trace = false;
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
129
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
130 if (nargin > 3)
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
131 if (! ischar (varargin{1}))
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
132 if (!isempty (varargin{1}))
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
133 abstol = varargin{1};
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
134 reltol = 0;
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
135 endif
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
136 if (nargin > 4)
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
137 trace = varargin{2};
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 > 5)
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
140 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
141 endif
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
142 else
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
143 idx = 1;
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
144 while (idx < nargin - 3)
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
145 if (ischar (varargin{idx}))
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
146 str = varargin{idx++};
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
147 if (strcmpi (str, "reltol"))
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
148 reltol = varargin{idx++};
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
149 elseif (strcmpi (str, "abstol"))
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
150 abstol = varargin{idx++};
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
151 elseif (strcmpi (str, "waypoints"))
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
152 waypoints = varargin{idx++} (:);
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
153 if (isreal(waypoints))
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
154 waypoints (waypoints < a | waypoints > b) = [];
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
155 endif
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
156 elseif (strcmpi (str, "maxintervalcount"))
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
157 maxint = varargin{idx++};
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
158 elseif (strcmpi (str, "trace"))
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
159 trace = varargin{idx++};
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
160 else
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
161 error ("quadgk: unknown property %s", str);
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
162 endif
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: expecting property to be a string");
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 endwhile
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
167 if (idx != nargin - 2)
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
168 error ("quadgk: expecting properties in pairs")
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
169 endif
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
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
173 ## 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
174 if (ischar (f))
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
175 f = @(x) feval (f, x);
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
176 endif
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
177
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
178 ## 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
179 ## 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
180 ## contour integrals
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
181 if (iscomplex (a) || iscomplex (b) || iscomplex(waypoints))
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
182 ## contour integral, no transform
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
183 subs = [a; waypoints; b];
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
184 h = b - a;
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
185 trans = @(t) t;
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
186 elseif (isinf (a) && isinf(b))
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
187 ## Standard Infinite to finite integral transformation.
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
188 ## \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
189 ## where
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
190 ## g(t) = t / (1 - t^2)
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
191 ## g'(t) = 1 / (1 + t^2) ^ 2
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
192 ## waypoint transform is then
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
193 ## 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
194 if (!isempty (waypoints))
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
195 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
196 subs = [-1; trans(waypoints); 1];
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
197 else
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
198 subs = linspace (-1, 1, 11)';
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
199 endif
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
200 h = 2;
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
201 trans = @(t) t ./ (1 - t.^2);
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
202 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
203 elseif (isinf(a))
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
204 ## 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
205 ## 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
206 ## a finite interval. The singularity weakening transform is
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
207 ## \int_{-\infinity}^b f(x) dx =
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
208 ## - \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
209 ## (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
210 ## \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
211 ## \int_{-1}^0 f (b - g(s) ^ 2 ) 2 g(s) g'(s) ds
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
212 ## where
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
213 ## g(s) = s / (1 + s)
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
214 ## g'(s) = 1 / (1 + s) ^ 2
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
215 ## waypoint transform is then
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
216 ## t = sqrt (b - x)
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
217 ## s = - t / (t + 1)
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
218 if (!isempty (waypoints))
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
219 tmp = sqrt (b - waypoints);
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
220 trans = @(x) - x ./ (x + 1);
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
221 subs = [0; trans(tmp); 1];
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
222 else
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
223 subs = linspace (0, 1, 11)';
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
224 endif
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
225 h = 1;
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
226 trans = @(t) b - (t ./ (1 + t)).^2;
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
227 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
228 elseif (isinf(b))
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
229 ## 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
230 ## 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
231 ## a finite interval. The singularity weakening transform is
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
232 ## \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
233 ## and the finite interval transform is
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
234 ## \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
235 ## \int_0^1 f (a + g(s) ^ 2 ) 2 g(s) g'(s) ds
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
236 ## where
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
237 ## g(s) = s / (1 - s)
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
238 ## g'(s) = 1 / (1 - s) ^ 2
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
239 ## waypoint transform is then
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
240 ## t = sqrt (x - a)
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
241 ## s = t / (t + 1)
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
242 if (!isempty (waypoints))
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
243 tmp = sqrt (waypoints - a);
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
244 trans = @(x) x ./ (x + 1);
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
245 subs = [0; trans(tmp); 1];
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
246 else
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
247 subs = linspace (0, 1, 11)';
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
248 endif
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
249 h = 1;
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
250 trans = @(t) a + (t ./ (1 - t)).^2;
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
251 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
252 else
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
253 ## Davis, Rabinowitz, "Methods of Numerical Integration" p441 2ed.
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
254 ## 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
255 ## 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
256 ## g'(t) = ((b-a)/4) * (3 - 3t^2);
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
257 ## 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
258 ## 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
259 ## 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
260 ## ignored. The third is
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
261 ## c = (-4*x + 2*(b+a)) / (b-a);
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
262 ## 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
263 ## 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
264 if (! isempty (waypoints))
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
265 trans = @__quadgk_finite_waypoint__;
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
266 subs = [-1; trans(waypoints, a, b); 1];
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
267 else
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
268 subs = linspace(-1, 1, 11)';
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
269 endif
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
270 h = 2;
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
271 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
272 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
273 3 .* (b - a) ./ 4 .* (1 - t.^2);
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
274 endif
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
275
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
276 ## 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
277 ## 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
278 while (length (subs) < 11)
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
279 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
280 endwhile
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
281 subs = [subs(1:end-1), subs(2:end)];
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
282
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
283 warn_state = warning ("query", "Octave:divide-by-zero");
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
284
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
285 unwind_protect
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
286 ## Singularity will cause divide by zero warnings
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
287 warning ("off", "Octave:divide-by-zero")
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
288
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
289 ## 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
290 [q_subs, q_errs] = __quadgk_eval__ (f, subs);
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
291 q0 = sum (q_subs);
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
292 err0 = sum (q_errs);
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
293
7795
df9519e9990c Handle single precision eps values
David Bateman <dbateman@free.fr>
parents: 7771
diff changeset
294 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
295 myeps = eps ("single");
df9519e9990c Handle single precision eps values
David Bateman <dbateman@free.fr>
parents: 7771
diff changeset
296 else
df9519e9990c Handle single precision eps values
David Bateman <dbateman@free.fr>
parents: 7771
diff changeset
297 myeps = eps;
df9519e9990c Handle single precision eps values
David Bateman <dbateman@free.fr>
parents: 7771
diff changeset
298 endif
df9519e9990c Handle single precision eps values
David Bateman <dbateman@free.fr>
parents: 7771
diff changeset
299
7771
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
300 first = true;
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
301 while (true)
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
302 ## 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
303 ## 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
304 ## value for this test. Shampine suggests 100*eps
7795
df9519e9990c Handle single precision eps values
David Bateman <dbateman@free.fr>
parents: 7771
diff changeset
305 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
306 q = q0;
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
307 err = err0;
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
308 break;
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
309 endif
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
310
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
311 ## 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
312 if (any (! isfinite (q_subs)))
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
313 warning ("quadgk: non finite integrand encountered");
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
314 q = q0;
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
315 err = err0;
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
316 break;
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
317 endif
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
318
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
319 tol = max (abstol, reltol .* abs (q0));
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
320
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
321 ## If the global error estimate is meet exit
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
322 if (err0 < tol)
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
323 q = q0;
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
324 err = err0;
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
325 break;
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
326 endif
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
327
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
328 ## 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
329 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
330 if (first)
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
331 q = sum (q_subs (idx));
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
332 err = sum (q_errs(idx));
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
333 first = false;
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
334 else
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
335 q0 = q + sum (q_subs);
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
336 err0 = err + sum (q_errs);
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
337 q += sum (q_subs (idx));
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
338 err += sum (q_errs(idx));
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
339 endif
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
340 subs(idx,:) = [];
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
341
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
342 ## If no remaining sub-intervals exit
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
343 if (rows (subs) == 0)
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
344 break;
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
345 endif
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
346
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
347 if (trace)
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
348 disp([fcnt, rows(subs), err, q0]);
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
349 endif
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
350
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
351 ## Split remaining sub-intervals in two
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
352 mid = (subs(:,2) + subs(:,1)) ./ 2;
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
353 subs = [subs(:,1), mid; mid, subs(:,2)];
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
354
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
355 ## 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
356 ## sub-interval and exit
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
357 if (rows (subs) > maxint)
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
358 warning ("quadgk: maximum interval count (%d) met", maxint);
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
359 q += sum (q_subs);
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
360 err += sum (q_errs);
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
361 break;
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
362 endif
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
363
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
364 ## 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
365 [q_subs, q_errs] = __quadgk_eval__ (f, subs);
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
366 endwhile
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
367
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
368 if (err > max (abstol, reltol * abs(q)))
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
369 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
370 endif
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
371 unwind_protect_cleanup
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
372 if (strcmp (warn_state.state, "on"))
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
373 warning ("on", "Octave:divide-by-zero")
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
374 endif
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
375 end_unwind_protect
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
376 endif
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
377 endfunction
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
378
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
379 function [q, err] = __quadgk_eval__ (f, subs)
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
380 ## 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
381 ## 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
382
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
383 persistent abscissa = [-0.9914553711208126e+00, -0.9491079123427585e+00, ...
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
384 -0.8648644233597691e+00, -0.7415311855993944e+00, ...
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
385 -0.5860872354676911e+00, -0.4058451513773972e+00, ...
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
386 -0.2077849550078985e+00, 0.0000000000000000e+00, ...
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
387 0.2077849550078985e+00, 0.4058451513773972e+00, ...
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
388 0.5860872354676911e+00, 0.7415311855993944e+00, ...
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
389 0.8648644233597691e+00, 0.9491079123427585e+00, ...
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
390 0.9914553711208126e+00];
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
391
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
392 persistent weights15 = ...
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
393 diag ([0.2293532201052922e-01, 0.6309209262997855e-01, ...
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
394 0.1047900103222502e+00, 0.1406532597155259e+00, ...
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
395 0.1690047266392679e+00, 0.1903505780647854e+00, ...
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
396 0.2044329400752989e+00, 0.2094821410847278e+00, ...
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
397 0.2044329400752989e+00, 0.1903505780647854e+00, ...
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
398 0.1690047266392679e+00, 0.1406532597155259e+00, ...
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
399 0.1047900103222502e+00, 0.6309209262997855e-01, ...
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
400 0.2293532201052922e-01]);
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
401
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
402 persistent weights7 = ...
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
403 diag ([0.1294849661688697e+00, 0.2797053914892767e+00, ...
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
404 0.3818300505051889e+00, 0.4179591836734694e+00, ...
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
405 0.3818300505051889e+00, 0.2797053914892767e+00, ...
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
406 0.1294849661688697e+00]);
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
407
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
408 halfwidth = diff (subs, [], 2) ./ 2;
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
409 center = sum (subs, 2) ./ 2;;
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
410 x = bsxfun (@plus, halfwidth * abscissa, center);
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
411 y = reshape (f (x(:)), size(x));
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
412
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
413 ## 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
414 ## 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
415 ## function.
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
416 q = sum (y * weights15, 2) .* halfwidth;
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
417 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
418 endfunction
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
419
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
420 function t = __quadgk_finite_waypoint__ (x, a, b)
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
421 c = (-4 .* x + 2.* (b + a)) ./ (b - a);
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
422 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
423 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
424 endfunction
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
425
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
426 %error (quadgk (@sin))
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
427 %error (quadgk (@sin, -pi))
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
428 %error (quadgk (@sin, -pi, pi, 'DummyArg'))
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
429
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
430 %!assert (quadgk(@sin,-pi,pi), 0, 1e-6)
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
431 %!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
432 %!assert (quadgk('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,'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
434 %!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
435
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
436 %!assert (quadgk(@sin,-pi,0), -2, 1e-6)
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
437 %!assert (quadgk(@sin,0,pi), 2, 1e-6)
680631e787aa Add quadv, quadgk, dblquad and triplequad functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
438 %!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
439 %!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
440 %!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
441 %!assert (quadgk (@(z) log (z), 1+1i, 1+1i, 'WayPoints', [1-1i, -1,-1i, -1+1i]), -pi * 1i, 1e-6)