Mercurial > octave-nkf
annotate doc/interpreter/quad.txi @ 7984:bbaa5d7d0143
Some documentation updates
| author | David Bateman <dbateman@free.fr> |
|---|---|
| date | Mon, 28 Jul 2008 15:47:40 +0200 |
| parents | fd42779a8428 |
| children | 23c248d415b5 |
| rev | line source |
|---|---|
| 7018 | 1 @c Copyright (C) 1996, 1997, 1999, 2002, 2007 John W. Eaton |
| 2 @c | |
| 3 @c This file is part of Octave. | |
| 4 @c | |
| 5 @c Octave is free software; you can redistribute it and/or modify it | |
| 6 @c under the terms of the GNU General Public License as published by the | |
| 7 @c Free Software Foundation; either version 3 of the License, or (at | |
| 8 @c your option) any later version. | |
| 9 @c | |
| 10 @c Octave is distributed in the hope that it will be useful, but WITHOUT | |
| 11 @c ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or | |
| 12 @c FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License | |
| 13 @c for more details. | |
| 14 @c | |
| 15 @c You should have received a copy of the GNU General Public License | |
| 16 @c along with Octave; see the file COPYING. If not, see | |
| 17 @c <http://www.gnu.org/licenses/>. | |
| 3294 | 18 |
| 6741 | 19 @node Numerical Integration |
| 20 @chapter Numerical Integration | |
| 21 | |
| 22 Octave comes with several built-in functions for computing the integral | |
| 23 of a function numerically. These functions all solve 1-dimensional | |
| 24 integration problems. | |
| 3294 | 25 |
| 26 @menu | |
| 6741 | 27 * Functions of One Variable:: |
| 28 * Functions of Multiple Variables:: | |
| 3294 | 29 * Orthogonal Collocation:: |
| 30 @end menu | |
| 31 | |
| 4167 | 32 @node Functions of One Variable |
| 3294 | 33 @section Functions of One Variable |
| 34 | |
| 6741 | 35 Octave supports three different algorithms for computing the integral |
| 36 @iftex | |
| 37 @tex | |
| 38 $$ | |
| 39 \int_a^b f(x) d x | |
| 40 $$ | |
| 41 @end tex | |
| 42 @end iftex | |
| 43 of a function @math{f} over the interval from @math{a} to @math{b}. | |
| 44 These are | |
| 45 | |
| 46 @table @code | |
| 47 @item quad | |
| 48 Numerical integration based on Gaussian quadrature. | |
| 49 | |
| 50 @item quadl | |
| 51 Numerical integration using an adaptive Lobatto rule. | |
| 52 | |
|
7984
bbaa5d7d0143
Some documentation updates
David Bateman <dbateman@free.fr>
parents:
7018
diff
changeset
|
53 @item quadgk |
|
bbaa5d7d0143
Some documentation updates
David Bateman <dbateman@free.fr>
parents:
7018
diff
changeset
|
54 Numerical integration using an adaptive Guass-Konrod rule. |
|
bbaa5d7d0143
Some documentation updates
David Bateman <dbateman@free.fr>
parents:
7018
diff
changeset
|
55 |
|
bbaa5d7d0143
Some documentation updates
David Bateman <dbateman@free.fr>
parents:
7018
diff
changeset
|
56 @item quadv |
|
bbaa5d7d0143
Some documentation updates
David Bateman <dbateman@free.fr>
parents:
7018
diff
changeset
|
57 Numerical integration using an adaptive vectorized Simpson's rule. |
|
bbaa5d7d0143
Some documentation updates
David Bateman <dbateman@free.fr>
parents:
7018
diff
changeset
|
58 |
| 6741 | 59 @item trapz |
| 6939 | 60 Numerical integration using the trapezoidal method. |
| 6741 | 61 @end table |
| 62 | |
| 63 @noindent | |
| 64 Besides these functions Octave also allows you to perform cumulative | |
| 6939 | 65 numerical integration using the trapezoidal method through the |
| 6741 | 66 @code{cumtrapz} function. |
| 67 | |
| 3368 | 68 @DOCSTRING(quad) |
| 3294 | 69 |
| 3368 | 70 @DOCSTRING(quad_options) |
| 3294 | 71 |
| 72 Here is an example of using @code{quad} to integrate the function | |
| 73 @iftex | |
| 74 @tex | |
| 75 $$ | |
| 76 f(x) = x \sin (1/x) \sqrt {|1 - x|} | |
| 77 $$ | |
| 78 from $x = 0$ to $x = 3$. | |
| 79 @end tex | |
| 80 @end iftex | |
| 6741 | 81 @ifnottex |
| 3294 | 82 |
| 83 @example | |
| 84 @var{f}(@var{x}) = @var{x} * sin (1/@var{x}) * sqrt (abs (1 - @var{x})) | |
| 85 @end example | |
| 86 | |
| 87 @noindent | |
| 88 from @var{x} = 0 to @var{x} = 3. | |
| 6741 | 89 @end ifnottex |
| 3294 | 90 |
| 91 This is a fairly difficult integration (plot the function over the range | |
| 92 of integration to see why). | |
| 93 | |
| 94 The first step is to define the function: | |
| 95 | |
| 96 @example | |
| 97 @group | |
| 98 function y = f (x) | |
| 99 y = x .* sin (1 ./ x) .* sqrt (abs (1 - x)); | |
| 100 endfunction | |
| 101 @end group | |
| 102 @end example | |
| 103 | |
| 104 Note the use of the `dot' forms of the operators. This is not necessary | |
| 105 for the call to @code{quad}, but it makes it much easier to generate a | |
| 106 set of points for plotting (because it makes it possible to call the | |
| 107 function with a vector argument to produce a vector result). | |
| 108 | |
| 109 Then we simply call quad: | |
| 110 | |
| 111 @example | |
| 112 @group | |
| 113 [v, ier, nfun, err] = quad ("f", 0, 3) | |
| 114 @result{} 1.9819 | |
| 115 @result{} 1 | |
| 116 @result{} 5061 | |
| 117 @result{} 1.1522e-07 | |
| 118 @end group | |
| 119 @end example | |
| 120 | |
| 121 Although @code{quad} returns a nonzero value for @var{ier}, the result | |
| 122 is reasonably accurate (to see why, examine what happens to the result | |
| 123 if you move the lower bound to 0.1, then 0.01, then 0.001, etc.). | |
| 124 | |
| 6502 | 125 @DOCSTRING(quadl) |
| 126 | |
|
7984
bbaa5d7d0143
Some documentation updates
David Bateman <dbateman@free.fr>
parents:
7018
diff
changeset
|
127 @DOCSTRING(quadgk) |
|
bbaa5d7d0143
Some documentation updates
David Bateman <dbateman@free.fr>
parents:
7018
diff
changeset
|
128 |
|
bbaa5d7d0143
Some documentation updates
David Bateman <dbateman@free.fr>
parents:
7018
diff
changeset
|
129 @DOCSTRING(quadv) |
|
bbaa5d7d0143
Some documentation updates
David Bateman <dbateman@free.fr>
parents:
7018
diff
changeset
|
130 |
| 6502 | 131 @DOCSTRING(trapz) |
| 132 | |
| 133 @DOCSTRING(cumtrapz) | |
| 134 | |
| 4167 | 135 @node Orthogonal Collocation |
| 3294 | 136 @section Orthogonal Collocation |
| 137 | |
| 3368 | 138 @DOCSTRING(colloc) |
| 3294 | 139 |
| 140 Here is an example of using @code{colloc} to generate weight matrices | |
| 141 for solving the second order differential equation | |
| 142 @iftex | |
| 143 @tex | |
| 144 $u^\prime - \alpha u^{\prime\prime} = 0$ with the boundary conditions | |
| 145 $u(0) = 0$ and $u(1) = 1$. | |
| 146 @end tex | |
| 147 @end iftex | |
| 6741 | 148 @ifnottex |
| 3294 | 149 @var{u}' - @var{alpha} * @var{u}'' = 0 with the boundary conditions |
| 150 @var{u}(0) = 0 and @var{u}(1) = 1. | |
| 6741 | 151 @end ifnottex |
| 3294 | 152 |
| 153 First, we can generate the weight matrices for @var{n} points (including | |
| 154 the endpoints of the interval), and incorporate the boundary conditions | |
| 155 in the right hand side (for a specific value of | |
| 156 @iftex | |
| 157 @tex | |
| 158 $\alpha$). | |
| 159 @end tex | |
| 160 @end iftex | |
| 6741 | 161 @ifnottex |
| 3294 | 162 @var{alpha}). |
| 6741 | 163 @end ifnottex |
| 3294 | 164 |
| 165 @example | |
| 166 @group | |
| 167 n = 7; | |
| 168 alpha = 0.1; | |
| 169 [r, a, b] = colloc (n-2, "left", "right"); | |
| 170 at = a(2:n-1,2:n-1); | |
| 171 bt = b(2:n-1,2:n-1); | |
| 172 rhs = alpha * b(2:n-1,n) - a(2:n-1,n); | |
| 173 @end group | |
| 174 @end example | |
| 175 | |
| 176 Then the solution at the roots @var{r} is | |
| 177 | |
| 178 @example | |
| 179 u = [ 0; (at - alpha * bt) \ rhs; 1] | |
| 180 @result{} [ 0.00; 0.004; 0.01 0.00; 0.12; 0.62; 1.00 ] | |
| 181 @end example | |
| 6741 | 182 |
| 183 @node Functions of Multiple Variables | |
| 184 @section Functions of Multiple Variables | |
| 185 | |
|
7984
bbaa5d7d0143
Some documentation updates
David Bateman <dbateman@free.fr>
parents:
7018
diff
changeset
|
186 Octave does not have built-in functions for computing the integral of |
|
bbaa5d7d0143
Some documentation updates
David Bateman <dbateman@free.fr>
parents:
7018
diff
changeset
|
187 functions of multiple variables directly. It is however possible to |
|
bbaa5d7d0143
Some documentation updates
David Bateman <dbateman@free.fr>
parents:
7018
diff
changeset
|
188 compute the integral of a function of multiple variables using the |
|
bbaa5d7d0143
Some documentation updates
David Bateman <dbateman@free.fr>
parents:
7018
diff
changeset
|
189 functions for one-dimensional integrals. |
| 6741 | 190 |
| 191 To illustrate how the integration can be performed, we will integrate | |
| 192 the function | |
| 193 @iftex | |
| 194 @tex | |
| 195 $$ | |
| 196 f(x, y) = \sin(\pi x y)\sqrt{x y} | |
| 197 $$ | |
| 198 @end tex | |
| 199 @end iftex | |
| 200 @ifnottex | |
| 201 @example | |
| 202 f(x, y) = sin(pi*x*y)*sqrt(x*y) | |
| 203 @end example | |
| 204 @end ifnottex | |
| 205 for @math{x} and @math{y} between 0 and 1. | |
| 206 | |
| 207 The first approach creates a function that integrates @math{f} with | |
| 208 respect to @math{x}, and then integrates that function with respect to | |
| 209 @math{y}. Since @code{quad} is written in Fortran it cannot be called | |
| 210 recursively. This means that @code{quad} cannot integrate a function | |
| 211 that calls @code{quad}, and hence cannot be used to perform the double | |
| 212 integration. It is however possible with @code{quadl}, which is what | |
| 213 the following code does. | |
| 214 | |
| 215 @example | |
| 216 function I = g(y) | |
| 217 I = ones(1, length(y)); | |
| 218 for i = 1:length(y) | |
| 219 f = @@(x) sin(pi.*x.*y(i)).*sqrt(x.*y(i)); | |
| 220 I(i) = quadl(f, 0, 1); | |
| 221 endfor | |
| 222 endfunction | |
| 223 | |
| 224 I = quadl("g", 0, 1) | |
| 225 @result{} 0.30022 | |
| 226 @end example | |
| 227 | |
|
7984
bbaa5d7d0143
Some documentation updates
David Bateman <dbateman@free.fr>
parents:
7018
diff
changeset
|
228 The above process can be simplified with the @code{dblquad} and |
|
bbaa5d7d0143
Some documentation updates
David Bateman <dbateman@free.fr>
parents:
7018
diff
changeset
|
229 @code{triplequad} functions for integrals over two and three |
|
bbaa5d7d0143
Some documentation updates
David Bateman <dbateman@free.fr>
parents:
7018
diff
changeset
|
230 variables. For example |
|
bbaa5d7d0143
Some documentation updates
David Bateman <dbateman@free.fr>
parents:
7018
diff
changeset
|
231 |
|
bbaa5d7d0143
Some documentation updates
David Bateman <dbateman@free.fr>
parents:
7018
diff
changeset
|
232 @example |
|
bbaa5d7d0143
Some documentation updates
David Bateman <dbateman@free.fr>
parents:
7018
diff
changeset
|
233 I = dblquad (@(x, y) sin(pi.*x.*y).*sqrt(x.*y), 0, 1, 0, 1) |
|
bbaa5d7d0143
Some documentation updates
David Bateman <dbateman@free.fr>
parents:
7018
diff
changeset
|
234 @result{} 0.30022 |
|
bbaa5d7d0143
Some documentation updates
David Bateman <dbateman@free.fr>
parents:
7018
diff
changeset
|
235 @end example |
|
bbaa5d7d0143
Some documentation updates
David Bateman <dbateman@free.fr>
parents:
7018
diff
changeset
|
236 |
|
bbaa5d7d0143
Some documentation updates
David Bateman <dbateman@free.fr>
parents:
7018
diff
changeset
|
237 @DOCSTRING(dblquad) |
|
bbaa5d7d0143
Some documentation updates
David Bateman <dbateman@free.fr>
parents:
7018
diff
changeset
|
238 |
|
bbaa5d7d0143
Some documentation updates
David Bateman <dbateman@free.fr>
parents:
7018
diff
changeset
|
239 @DOCSTRING(triplequad) |
|
bbaa5d7d0143
Some documentation updates
David Bateman <dbateman@free.fr>
parents:
7018
diff
changeset
|
240 |
| 6741 | 241 The above mentioned approach works but is fairly slow, and that problem |
| 242 increases exponentially with the dimensionality the problem. Another | |
| 243 possible solution is to use Orthogonal Collocation as described in the | |
| 244 previous section. The integral of a function @math{f(x,y)} for | |
| 245 @math{x} and @math{y} between 0 and 1 can be approximated using @math{n} | |
| 246 points by | |
| 247 @iftex | |
| 248 @tex | |
| 249 $$ | |
| 250 \int_0^1 \int_0^1 f(x,y) d x d y \approx \sum_{i=1}^n \sum_{j=1}^n q_i q_j f(r_i, r_j), | |
| 251 $$ | |
| 252 @end tex | |
| 253 @end iftex | |
| 254 @ifnottex | |
| 255 the sum over @code{i=1:n} and @code{j=1:n} of @code{q(i)*q(j)*f(r(i),r(j))}, | |
| 256 @end ifnottex | |
| 257 where @math{q} and @math{r} is as returned by @code{colloc(n)}. The | |
| 258 generalisation to more than two variables is straight forward. The | |
| 259 following code computes the studied integral using @math{n=7} points. | |
| 260 | |
| 261 @example | |
| 262 f = @@(x,y) sin(pi*x*y').*sqrt(x*y'); | |
| 263 n = 7; | |
| 264 [t, A, B, q] = colloc(n); | |
| 265 I = q'*f(t,t)*q; | |
| 266 @result{} 0.30022 | |
| 267 @end example | |
| 268 | |
| 269 @noindent | |
| 270 It should be noted that the number of points determines the quality | |
| 271 of the approximation. If the integration needs to be performed between | |
| 272 @math{a} and @math{b} instead of 0 and 1, a change of variables is needed. | |
| 273 | |
| 274 | |
| 275 |