|
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 |
|
|
53 @item trapz |
|
6939
|
54 Numerical integration using the trapezoidal method. |
|
6741
|
55 @end table |
|
|
56 |
|
|
57 @noindent |
|
|
58 Besides these functions Octave also allows you to perform cumulative |
|
6939
|
59 numerical integration using the trapezoidal method through the |
|
6741
|
60 @code{cumtrapz} function. |
|
|
61 |
|
3368
|
62 @DOCSTRING(quad) |
|
3294
|
63 |
|
3368
|
64 @DOCSTRING(quad_options) |
|
3294
|
65 |
|
|
66 Here is an example of using @code{quad} to integrate the function |
|
|
67 @iftex |
|
|
68 @tex |
|
|
69 $$ |
|
|
70 f(x) = x \sin (1/x) \sqrt {|1 - x|} |
|
|
71 $$ |
|
|
72 from $x = 0$ to $x = 3$. |
|
|
73 @end tex |
|
|
74 @end iftex |
|
6741
|
75 @ifnottex |
|
3294
|
76 |
|
|
77 @example |
|
|
78 @var{f}(@var{x}) = @var{x} * sin (1/@var{x}) * sqrt (abs (1 - @var{x})) |
|
|
79 @end example |
|
|
80 |
|
|
81 @noindent |
|
|
82 from @var{x} = 0 to @var{x} = 3. |
|
6741
|
83 @end ifnottex |
|
3294
|
84 |
|
|
85 This is a fairly difficult integration (plot the function over the range |
|
|
86 of integration to see why). |
|
|
87 |
|
|
88 The first step is to define the function: |
|
|
89 |
|
|
90 @example |
|
|
91 @group |
|
|
92 function y = f (x) |
|
|
93 y = x .* sin (1 ./ x) .* sqrt (abs (1 - x)); |
|
|
94 endfunction |
|
|
95 @end group |
|
|
96 @end example |
|
|
97 |
|
|
98 Note the use of the `dot' forms of the operators. This is not necessary |
|
|
99 for the call to @code{quad}, but it makes it much easier to generate a |
|
|
100 set of points for plotting (because it makes it possible to call the |
|
|
101 function with a vector argument to produce a vector result). |
|
|
102 |
|
|
103 Then we simply call quad: |
|
|
104 |
|
|
105 @example |
|
|
106 @group |
|
|
107 [v, ier, nfun, err] = quad ("f", 0, 3) |
|
|
108 @result{} 1.9819 |
|
|
109 @result{} 1 |
|
|
110 @result{} 5061 |
|
|
111 @result{} 1.1522e-07 |
|
|
112 @end group |
|
|
113 @end example |
|
|
114 |
|
|
115 Although @code{quad} returns a nonzero value for @var{ier}, the result |
|
|
116 is reasonably accurate (to see why, examine what happens to the result |
|
|
117 if you move the lower bound to 0.1, then 0.01, then 0.001, etc.). |
|
|
118 |
|
6502
|
119 @DOCSTRING(quadl) |
|
|
120 |
|
|
121 @DOCSTRING(trapz) |
|
|
122 |
|
|
123 @DOCSTRING(cumtrapz) |
|
|
124 |
|
4167
|
125 @node Orthogonal Collocation |
|
3294
|
126 @section Orthogonal Collocation |
|
|
127 |
|
3368
|
128 @DOCSTRING(colloc) |
|
3294
|
129 |
|
|
130 Here is an example of using @code{colloc} to generate weight matrices |
|
|
131 for solving the second order differential equation |
|
|
132 @iftex |
|
|
133 @tex |
|
|
134 $u^\prime - \alpha u^{\prime\prime} = 0$ with the boundary conditions |
|
|
135 $u(0) = 0$ and $u(1) = 1$. |
|
|
136 @end tex |
|
|
137 @end iftex |
|
6741
|
138 @ifnottex |
|
3294
|
139 @var{u}' - @var{alpha} * @var{u}'' = 0 with the boundary conditions |
|
|
140 @var{u}(0) = 0 and @var{u}(1) = 1. |
|
6741
|
141 @end ifnottex |
|
3294
|
142 |
|
|
143 First, we can generate the weight matrices for @var{n} points (including |
|
|
144 the endpoints of the interval), and incorporate the boundary conditions |
|
|
145 in the right hand side (for a specific value of |
|
|
146 @iftex |
|
|
147 @tex |
|
|
148 $\alpha$). |
|
|
149 @end tex |
|
|
150 @end iftex |
|
6741
|
151 @ifnottex |
|
3294
|
152 @var{alpha}). |
|
6741
|
153 @end ifnottex |
|
3294
|
154 |
|
|
155 @example |
|
|
156 @group |
|
|
157 n = 7; |
|
|
158 alpha = 0.1; |
|
|
159 [r, a, b] = colloc (n-2, "left", "right"); |
|
|
160 at = a(2:n-1,2:n-1); |
|
|
161 bt = b(2:n-1,2:n-1); |
|
|
162 rhs = alpha * b(2:n-1,n) - a(2:n-1,n); |
|
|
163 @end group |
|
|
164 @end example |
|
|
165 |
|
|
166 Then the solution at the roots @var{r} is |
|
|
167 |
|
|
168 @example |
|
|
169 u = [ 0; (at - alpha * bt) \ rhs; 1] |
|
|
170 @result{} [ 0.00; 0.004; 0.01 0.00; 0.12; 0.62; 1.00 ] |
|
|
171 @end example |
|
6741
|
172 |
|
|
173 @node Functions of Multiple Variables |
|
|
174 @section Functions of Multiple Variables |
|
|
175 |
|
|
176 Octave does not have built-in functions for computing the integral |
|
|
177 of functions of multiple variables. It is however possible to compute |
|
|
178 the integral of a function of multiple variables using the functions |
|
|
179 for one-dimensional integrals. |
|
|
180 |
|
|
181 To illustrate how the integration can be performed, we will integrate |
|
|
182 the function |
|
|
183 @iftex |
|
|
184 @tex |
|
|
185 $$ |
|
|
186 f(x, y) = \sin(\pi x y)\sqrt{x y} |
|
|
187 $$ |
|
|
188 @end tex |
|
|
189 @end iftex |
|
|
190 @ifnottex |
|
|
191 @example |
|
|
192 f(x, y) = sin(pi*x*y)*sqrt(x*y) |
|
|
193 @end example |
|
|
194 @end ifnottex |
|
|
195 for @math{x} and @math{y} between 0 and 1. |
|
|
196 |
|
|
197 The first approach creates a function that integrates @math{f} with |
|
|
198 respect to @math{x}, and then integrates that function with respect to |
|
|
199 @math{y}. Since @code{quad} is written in Fortran it cannot be called |
|
|
200 recursively. This means that @code{quad} cannot integrate a function |
|
|
201 that calls @code{quad}, and hence cannot be used to perform the double |
|
|
202 integration. It is however possible with @code{quadl}, which is what |
|
|
203 the following code does. |
|
|
204 |
|
|
205 @example |
|
|
206 function I = g(y) |
|
|
207 I = ones(1, length(y)); |
|
|
208 for i = 1:length(y) |
|
|
209 f = @@(x) sin(pi.*x.*y(i)).*sqrt(x.*y(i)); |
|
|
210 I(i) = quadl(f, 0, 1); |
|
|
211 endfor |
|
|
212 endfunction |
|
|
213 |
|
|
214 I = quadl("g", 0, 1) |
|
|
215 @result{} 0.30022 |
|
|
216 @end example |
|
|
217 |
|
|
218 The above mentioned approach works but is fairly slow, and that problem |
|
|
219 increases exponentially with the dimensionality the problem. Another |
|
|
220 possible solution is to use Orthogonal Collocation as described in the |
|
|
221 previous section. The integral of a function @math{f(x,y)} for |
|
|
222 @math{x} and @math{y} between 0 and 1 can be approximated using @math{n} |
|
|
223 points by |
|
|
224 @iftex |
|
|
225 @tex |
|
|
226 $$ |
|
|
227 \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), |
|
|
228 $$ |
|
|
229 @end tex |
|
|
230 @end iftex |
|
|
231 @ifnottex |
|
|
232 the sum over @code{i=1:n} and @code{j=1:n} of @code{q(i)*q(j)*f(r(i),r(j))}, |
|
|
233 @end ifnottex |
|
|
234 where @math{q} and @math{r} is as returned by @code{colloc(n)}. The |
|
|
235 generalisation to more than two variables is straight forward. The |
|
|
236 following code computes the studied integral using @math{n=7} points. |
|
|
237 |
|
|
238 @example |
|
|
239 f = @@(x,y) sin(pi*x*y').*sqrt(x*y'); |
|
|
240 n = 7; |
|
|
241 [t, A, B, q] = colloc(n); |
|
|
242 I = q'*f(t,t)*q; |
|
|
243 @result{} 0.30022 |
|
|
244 @end example |
|
|
245 |
|
|
246 @noindent |
|
|
247 It should be noted that the number of points determines the quality |
|
|
248 of the approximation. If the integration needs to be performed between |
|
|
249 @math{a} and @math{b} instead of 0 and 1, a change of variables is needed. |
|
|
250 |
|
|
251 |
|
|
252 |
|
|
253 |
