Mercurial > hg > octave-nkf

comparison doc/interpreter/quad.txi @ 6741:00116015904d

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
[project @ 2007-06-18 16:07:14 by jwe]
author jwe
date Mon, 18 Jun 2007 16:07:14 +0000
parents 6ab0a8767780
children 083721ae3dfa
comparison
equal deleted inserted replaced
6740:605ea655366d 6741:00116015904d
1 @c Copyright (C) 1996, 1997 John W. Eaton 1 @c Copyright (C) 1996, 1997 John W. Eaton
2 @c This is part of the Octave manual. 2 @c This is part of the Octave manual.
3 @c For copying conditions, see the file gpl.texi. 3 @c For copying conditions, see the file gpl.texi.
4 4
5 @node Quadrature 5 @node Numerical Integration
6 @chapter Quadrature 6 @chapter Numerical Integration
7
8 Octave comes with several built-in functions for computing the integral
9 of a function numerically. These functions all solve 1-dimensional
10 integration problems.
7 11
8 @menu 12 @menu
9 * Functions of One Variable:: 13 * Functions of One Variable::
14 * Functions of Multiple Variables::
10 * Orthogonal Collocation:: 15 * Orthogonal Collocation::
11 @end menu 16 @end menu
12 17
13 @node Functions of One Variable 18 @node Functions of One Variable
14 @section Functions of One Variable 19 @section Functions of One Variable
15 20
21 Octave supports three different algorithms for computing the integral
22 @iftex
23 @tex
24 $$
25 \int_a^b f(x) d x
26 $$
27 @end tex
28 @end iftex
29 of a function @math{f} over the interval from @math{a} to @math{b}.
30 These are
31
32 @table @code
33 @item quad
34 Numerical integration based on Gaussian quadrature.
35
36 @item quadl
37 Numerical integration using an adaptive Lobatto rule.
38
39 @item trapz
40 Numerical integration using the trapezodial method.
41 @end table
42
43 @noindent
44 Besides these functions Octave also allows you to perform cumulative
45 numerical integration using the trapezodial method through the
46 @code{cumtrapz} function.
47
16 @DOCSTRING(quad) 48 @DOCSTRING(quad)
17 49
18 @DOCSTRING(quad_options) 50 @DOCSTRING(quad_options)
19 51
20 Here is an example of using @code{quad} to integrate the function 52 Here is an example of using @code{quad} to integrate the function
24 f(x) = x \sin (1/x) \sqrt {|1 - x|} 56 f(x) = x \sin (1/x) \sqrt {|1 - x|}
25 $$ 57 $$
26 from $x = 0$ to $x = 3$. 58 from $x = 0$ to $x = 3$.
27 @end tex 59 @end tex
28 @end iftex 60 @end iftex
29 @ifinfo 61 @ifnottex
30 62
31 @example 63 @example
32 @var{f}(@var{x}) = @var{x} * sin (1/@var{x}) * sqrt (abs (1 - @var{x})) 64 @var{f}(@var{x}) = @var{x} * sin (1/@var{x}) * sqrt (abs (1 - @var{x}))
33 @end example 65 @end example
34 66
35 @noindent 67 @noindent
36 from @var{x} = 0 to @var{x} = 3. 68 from @var{x} = 0 to @var{x} = 3.
37 @end ifinfo 69 @end ifnottex
38 70
39 This is a fairly difficult integration (plot the function over the range 71 This is a fairly difficult integration (plot the function over the range
40 of integration to see why). 72 of integration to see why).
41 73
42 The first step is to define the function: 74 The first step is to define the function:
87 @tex 119 @tex
88 $u^\prime - \alpha u^{\prime\prime} = 0$ with the boundary conditions 120 $u^\prime - \alpha u^{\prime\prime} = 0$ with the boundary conditions
89 $u(0) = 0$ and $u(1) = 1$. 121 $u(0) = 0$ and $u(1) = 1$.
90 @end tex 122 @end tex
91 @end iftex 123 @end iftex
92 @ifinfo 124 @ifnottex
93 @var{u}' - @var{alpha} * @var{u}'' = 0 with the boundary conditions 125 @var{u}' - @var{alpha} * @var{u}'' = 0 with the boundary conditions
94 @var{u}(0) = 0 and @var{u}(1) = 1. 126 @var{u}(0) = 0 and @var{u}(1) = 1.
95 @end ifinfo 127 @end ifnottex
96 128
97 First, we can generate the weight matrices for @var{n} points (including 129 First, we can generate the weight matrices for @var{n} points (including
98 the endpoints of the interval), and incorporate the boundary conditions 130 the endpoints of the interval), and incorporate the boundary conditions
99 in the right hand side (for a specific value of 131 in the right hand side (for a specific value of
100 @iftex 132 @iftex
101 @tex 133 @tex
102 $\alpha$). 134 $\alpha$).
103 @end tex 135 @end tex
104 @end iftex 136 @end iftex
105 @ifinfo 137 @ifnottex
106 @var{alpha}). 138 @var{alpha}).
107 @end ifinfo 139 @end ifnottex
108 140
109 @example 141 @example
110 @group 142 @group
111 n = 7; 143 n = 7;
112 alpha = 0.1; 144 alpha = 0.1;
121 153
122 @example 154 @example
123 u = [ 0; (at - alpha * bt) \ rhs; 1] 155 u = [ 0; (at - alpha * bt) \ rhs; 1]
124 @result{} [ 0.00; 0.004; 0.01 0.00; 0.12; 0.62; 1.00 ] 156 @result{} [ 0.00; 0.004; 0.01 0.00; 0.12; 0.62; 1.00 ]
125 @end example 157 @end example
158
159 @node Functions of Multiple Variables
160 @section Functions of Multiple Variables
161
162 Octave does not have built-in functions for computing the integral
163 of functions of multiple variables. It is however possible to compute
164 the integral of a function of multiple variables using the functions
165 for one-dimensional integrals.
166
167 To illustrate how the integration can be performed, we will integrate
168 the function
169 @iftex
170 @tex
171 $$
172 f(x, y) = \sin(\pi x y)\sqrt{x y}
173 $$
174 @end tex
175 @end iftex
176 @ifnottex
177 @example
178 f(x, y) = sin(pi*x*y)*sqrt(x*y)
179 @end example
180 @end ifnottex
181 for @math{x} and @math{y} between 0 and 1.
182
183 The first approach creates a function that integrates @math{f} with
184 respect to @math{x}, and then integrates that function with respect to
185 @math{y}. Since @code{quad} is written in Fortran it cannot be called
186 recursively. This means that @code{quad} cannot integrate a function
187 that calls @code{quad}, and hence cannot be used to perform the double
188 integration. It is however possible with @code{quadl}, which is what
189 the following code does.
190
191 @example
192 function I = g(y)
193 I = ones(1, length(y));
194 for i = 1:length(y)
195 f = @@(x) sin(pi.*x.*y(i)).*sqrt(x.*y(i));
196 I(i) = quadl(f, 0, 1);
197 endfor
198 endfunction
199
200 I = quadl("g", 0, 1)
201 @result{} 0.30022
202 @end example
203
204 The above mentioned approach works but is fairly slow, and that problem
205 increases exponentially with the dimensionality the problem. Another
206 possible solution is to use Orthogonal Collocation as described in the
207 previous section. The integral of a function @math{f(x,y)} for
208 @math{x} and @math{y} between 0 and 1 can be approximated using @math{n}
209 points by
210 @iftex
211 @tex
212 $$
213 \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),
214 $$
215 @end tex
216 @end iftex
217 @ifnottex
218 the sum over @code{i=1:n} and @code{j=1:n} of @code{q(i)*q(j)*f(r(i),r(j))},
219 @end ifnottex
220 where @math{q} and @math{r} is as returned by @code{colloc(n)}. The
221 generalisation to more than two variables is straight forward. The
222 following code computes the studied integral using @math{n=7} points.
223
224 @example
225 f = @@(x,y) sin(pi*x*y').*sqrt(x*y');
226 n = 7;
227 [t, A, B, q] = colloc(n);
228 I = q'*f(t,t)*q;
229 @result{} 0.30022
230 @end example
231
232 @noindent
233 It should be noted that the number of points determines the quality
234 of the approximation. If the integration needs to be performed between
235 @math{a} and @math{b} instead of 0 and 1, a change of variables is needed.
236
237
238
239