1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
|
<!DOCTYPE html>
<html>
<!-- Created by GNU Texinfo 7.1.1, https://www.gnu.org/software/texinfo/ -->
<head>
<meta http-equiv="Content-Type" content="text/html; charset=utf-8">
<title>Functions of Multiple Variables (GNU Octave (version 10.3.0))</title>
<meta name="description" content="Functions of Multiple Variables (GNU Octave (version 10.3.0))">
<meta name="keywords" content="Functions of Multiple Variables (GNU Octave (version 10.3.0))">
<meta name="resource-type" content="document">
<meta name="distribution" content="global">
<meta name="Generator" content="makeinfo">
<meta name="viewport" content="width=device-width,initial-scale=1">
<link href="index.html" rel="start" title="Top">
<link href="Concept-Index.html" rel="index" title="Concept Index">
<link href="index.html#SEC_Contents" rel="contents" title="Table of Contents">
<link href="Numerical-Integration.html" rel="up" title="Numerical Integration">
<link href="Orthogonal-Collocation.html" rel="prev" title="Orthogonal Collocation">
<style type="text/css">
<!--
a.copiable-link {visibility: hidden; text-decoration: none; line-height: 0em}
div.example {margin-left: 3.2em}
span:hover a.copiable-link {visibility: visible}
strong.def-name {font-family: monospace; font-weight: bold; font-size: larger}
-->
</style>
<link rel="stylesheet" type="text/css" href="octave.css">
</head>
<body lang="en">
<div class="section-level-extent" id="Functions-of-Multiple-Variables">
<div class="nav-panel">
<p>
Previous: <a href="Orthogonal-Collocation.html" accesskey="p" rel="prev">Orthogonal Collocation</a>, Up: <a href="Numerical-Integration.html" accesskey="u" rel="up">Numerical Integration</a> [<a href="index.html#SEC_Contents" title="Table of contents" rel="contents">Contents</a>][<a href="Concept-Index.html" title="Index" rel="index">Index</a>]</p>
</div>
<hr>
<h3 class="section" id="Functions-of-Multiple-Variables-1"><span>23.3 Functions of Multiple Variables<a class="copiable-link" href="#Functions-of-Multiple-Variables-1"> ¶</a></span></h3>
<p>Octave includes several functions for computing the integral of functions of
multiple variables. This procedure can generally be performed by creating a
function that integrates <em class="math">f</em> with respect to <em class="math">x</em>, and then integrates
that function with respect to <em class="math">y</em>. This procedure can be performed
manually using the following example which integrates the function:
</p>
<div class="example">
<pre class="example-preformatted">f(x, y) = sin(pi*x*y) * sqrt(x*y)
</pre></div>
<p>for <em class="math">x</em> and <em class="math">y</em> between 0 and 1.
</p>
<p>Using <code class="code">quadgk</code> in the example below, a double integration can be
performed. (Note that any of the 1-D quadrature functions can be used in this
fashion except for <code class="code">quad</code> since it is written in Fortran and cannot be
called recursively.)
</p>
<div class="example">
<div class="group"><pre class="example-preformatted">function q = g(y)
q = ones (size (y));
for i = 1:length (y)
f = @(x) sin (pi*x.*y(i)) .* sqrt (x.*y(i));
q(i) = quadgk (f, 0, 1);
endfor
endfunction
I = quadgk ("g", 0, 1)
⇒ 0.30022
</pre></div></div>
<p>The algorithm above is implemented in the function <code class="code">dblquad</code> for integrals
over two variables. The 3-D equivalent of this process is implemented in
<code class="code">triplequad</code> for integrals over three variables. As an example, the
result above can be replicated with a call to <code class="code">dblquad</code> as shown below.
</p>
<div class="example">
<div class="group"><pre class="example-preformatted">I = dblquad (@(x, y) sin (pi*x.*y) .* sqrt (x.*y), 0, 1, 0, 1)
⇒ 0.30022
</pre></div></div>
<a class="anchor" id="XREFdblquad"></a><span style="display:block; margin-top:-4.5ex;"> </span>
<dl class="first-deftypefn">
<dt class="deftypefn" id="index-dblquad"><span><code class="def-type"><var class="var">q</var> =</code> <strong class="def-name">dblquad</strong> <code class="def-code-arguments">(<var class="var">f</var>, <var class="var">xa</var>, <var class="var">xb</var>, <var class="var">ya</var>, <var class="var">yb</var>)</code><a class="copiable-link" href="#index-dblquad"> ¶</a></span></dt>
<dt class="deftypefnx def-cmd-deftypefn" id="index-dblquad-1"><span><code class="def-type"><var class="var">q</var> =</code> <strong class="def-name">dblquad</strong> <code class="def-code-arguments">(<var class="var">f</var>, <var class="var">xa</var>, <var class="var">xb</var>, <var class="var">ya</var>, <var class="var">yb</var>, <var class="var">tol</var>)</code><a class="copiable-link" href="#index-dblquad-1"> ¶</a></span></dt>
<dt class="deftypefnx def-cmd-deftypefn" id="index-dblquad-2"><span><code class="def-type"><var class="var">q</var> =</code> <strong class="def-name">dblquad</strong> <code class="def-code-arguments">(<var class="var">f</var>, <var class="var">xa</var>, <var class="var">xb</var>, <var class="var">ya</var>, <var class="var">yb</var>, <var class="var">tol</var>, <var class="var">quadf</var>)</code><a class="copiable-link" href="#index-dblquad-2"> ¶</a></span></dt>
<dt class="deftypefnx def-cmd-deftypefn" id="index-dblquad-3"><span><code class="def-type"><var class="var">q</var> =</code> <strong class="def-name">dblquad</strong> <code class="def-code-arguments">(<var class="var">f</var>, <var class="var">xa</var>, <var class="var">xb</var>, <var class="var">ya</var>, <var class="var">yb</var>, <var class="var">tol</var>, <var class="var">quadf</var>, …)</code><a class="copiable-link" href="#index-dblquad-3"> ¶</a></span></dt>
<dd><p>Numerically evaluate the double integral of <var class="var">f</var>.
</p>
<p><var class="var">f</var> is a function handle, inline function, or string containing the name
of the function to evaluate. The function <var class="var">f</var> must have the form
<em class="math">z = f(x,y)</em> where <var class="var">x</var> is a vector and <var class="var">y</var> is a scalar. It
should return a vector of the same length and orientation as <var class="var">x</var>.
</p>
<p><var class="var">xa</var>, <var class="var">ya</var> and <var class="var">xb</var>, <var class="var">yb</var> are the lower and upper limits of
integration for x and y respectively. The underlying integrator determines
whether infinite bounds are accepted.
</p>
<p>The optional argument <var class="var">tol</var> defines the absolute tolerance used to
integrate each sub-integral. The default value is 1e-6.
</p>
<p>The optional argument <var class="var">quadf</var> specifies which underlying integrator
function to use. Any choice but <code class="code">quad</code> is available and the default
is <code class="code">quadcc</code>.
</p>
<p>Additional arguments, are passed directly to <var class="var">f</var>. To use the default
value for <var class="var">tol</var> or <var class="var">quadf</var> one may pass <code class="code">':'</code> or an empty
matrix ([]).
</p>
<p><strong class="strong">See also:</strong> <a class="ref" href="#XREFintegral2">integral2</a>, <a class="ref" href="#XREFintegral3">integral3</a>, <a class="ref" href="#XREFtriplequad">triplequad</a>, <a class="ref" href="Functions-of-One-Variable.html#XREFquad">quad</a>, <a class="ref" href="Functions-of-One-Variable.html#XREFquadv">quadv</a>, <a class="ref" href="Functions-of-One-Variable.html#XREFquadl">quadl</a>, <a class="ref" href="Functions-of-One-Variable.html#XREFquadgk">quadgk</a>, <a class="ref" href="Functions-of-One-Variable.html#XREFquadcc">quadcc</a>, <a class="ref" href="Functions-of-One-Variable.html#XREFtrapz">trapz</a>.
</p></dd></dl>
<a class="anchor" id="XREFtriplequad"></a><span style="display:block; margin-top:-4.5ex;"> </span>
<dl class="first-deftypefn">
<dt class="deftypefn" id="index-triplequad"><span><code class="def-type"><var class="var">q</var> =</code> <strong class="def-name">triplequad</strong> <code class="def-code-arguments">(<var class="var">f</var>, <var class="var">xa</var>, <var class="var">xb</var>, <var class="var">ya</var>, <var class="var">yb</var>, <var class="var">za</var>, <var class="var">zb</var>)</code><a class="copiable-link" href="#index-triplequad"> ¶</a></span></dt>
<dt class="deftypefnx def-cmd-deftypefn" id="index-triplequad-1"><span><code class="def-type"><var class="var">q</var> =</code> <strong class="def-name">triplequad</strong> <code class="def-code-arguments">(<var class="var">f</var>, <var class="var">xa</var>, <var class="var">xb</var>, <var class="var">ya</var>, <var class="var">yb</var>, <var class="var">za</var>, <var class="var">zb</var>, <var class="var">tol</var>)</code><a class="copiable-link" href="#index-triplequad-1"> ¶</a></span></dt>
<dt class="deftypefnx def-cmd-deftypefn" id="index-triplequad-2"><span><code class="def-type"><var class="var">q</var> =</code> <strong class="def-name">triplequad</strong> <code class="def-code-arguments">(<var class="var">f</var>, <var class="var">xa</var>, <var class="var">xb</var>, <var class="var">ya</var>, <var class="var">yb</var>, <var class="var">za</var>, <var class="var">zb</var>, <var class="var">tol</var>, <var class="var">quadf</var>)</code><a class="copiable-link" href="#index-triplequad-2"> ¶</a></span></dt>
<dt class="deftypefnx def-cmd-deftypefn" id="index-triplequad-3"><span><code class="def-type"><var class="var">q</var> =</code> <strong class="def-name">triplequad</strong> <code class="def-code-arguments">(<var class="var">f</var>, <var class="var">xa</var>, <var class="var">xb</var>, <var class="var">ya</var>, <var class="var">yb</var>, <var class="var">za</var>, <var class="var">zb</var>, <var class="var">tol</var>, <var class="var">quadf</var>, …)</code><a class="copiable-link" href="#index-triplequad-3"> ¶</a></span></dt>
<dd><p>Numerically evaluate the triple integral of <var class="var">f</var>.
</p>
<p><var class="var">f</var> is a function handle, inline function, or string containing the name
of the function to evaluate. The function <var class="var">f</var> must have the form
<em class="math">w = f(x,y,z)</em> where either <var class="var">x</var> or <var class="var">y</var> is a vector and the
remaining inputs are scalars. It should return a vector of the same length
and orientation as <var class="var">x</var> or <var class="var">y</var>.
</p>
<p><var class="var">xa</var>, <var class="var">ya</var>, <var class="var">za</var> and <var class="var">xb</var>, <var class="var">yb</var>, <var class="var">zb</var> are the lower
and upper limits of integration for x, y, and z respectively. The
underlying integrator determines whether infinite bounds are accepted.
</p>
<p>The optional argument <var class="var">tol</var> defines the absolute tolerance used to
integrate each sub-integral. The default value is 1e-6.
</p>
<p>The optional argument <var class="var">quadf</var> specifies which underlying integrator
function to use. Any choice but <code class="code">quad</code> is available and the default
is <code class="code">quadcc</code>.
</p>
<p>Additional arguments, are passed directly to <var class="var">f</var>. To use the default
value for <var class="var">tol</var> or <var class="var">quadf</var> one may pass <code class="code">':'</code> or an empty
matrix ([]).
</p>
<p><strong class="strong">See also:</strong> <a class="ref" href="#XREFintegral3">integral3</a>, <a class="ref" href="#XREFintegral2">integral2</a>, <a class="ref" href="#XREFdblquad">dblquad</a>, <a class="ref" href="Functions-of-One-Variable.html#XREFquad">quad</a>, <a class="ref" href="Functions-of-One-Variable.html#XREFquadv">quadv</a>, <a class="ref" href="Functions-of-One-Variable.html#XREFquadl">quadl</a>, <a class="ref" href="Functions-of-One-Variable.html#XREFquadgk">quadgk</a>, <a class="ref" href="Functions-of-One-Variable.html#XREFquadcc">quadcc</a>, <a class="ref" href="Functions-of-One-Variable.html#XREFtrapz">trapz</a>.
</p></dd></dl>
<p>The recursive algorithm for quadrature presented above is referred to as
<code class="code">"iterated"</code>. A separate 2-D integration method is implemented in the
function <code class="code">quad2d</code>. This function performs a <code class="code">"tiled"</code> integration
by subdividing the integration domain into rectangular regions and performing
separate integrations over those domains. The domains are further subdivided
in areas requiring refinement to reach the desired numerical accuracy. For
certain functions this method can be faster than the 2-D iteration used in the
other functions above.
</p>
<a class="anchor" id="XREFquad2d"></a><span style="display:block; margin-top:-4.5ex;"> </span>
<dl class="first-deftypefn">
<dt class="deftypefn" id="index-quad2d"><span><code class="def-type"><var class="var">q</var> =</code> <strong class="def-name">quad2d</strong> <code class="def-code-arguments">(<var class="var">f</var>, <var class="var">xa</var>, <var class="var">xb</var>, <var class="var">ya</var>, <var class="var">yb</var>)</code><a class="copiable-link" href="#index-quad2d"> ¶</a></span></dt>
<dt class="deftypefnx def-cmd-deftypefn" id="index-quad2d-1"><span><code class="def-type"><var class="var">q</var> =</code> <strong class="def-name">quad2d</strong> <code class="def-code-arguments">(<var class="var">f</var>, <var class="var">xa</var>, <var class="var">xb</var>, <var class="var">ya</var>, <var class="var">yb</var>, <var class="var">prop</var>, <var class="var">val</var>, …)</code><a class="copiable-link" href="#index-quad2d-1"> ¶</a></span></dt>
<dt class="deftypefnx def-cmd-deftypefn" id="index-quad2d-2"><span><code class="def-type">[<var class="var">q</var>, <var class="var">err</var>, <var class="var">iter</var>] =</code> <strong class="def-name">quad2d</strong> <code class="def-code-arguments">(…)</code><a class="copiable-link" href="#index-quad2d-2"> ¶</a></span></dt>
<dd>
<p>Numerically evaluate the two-dimensional integral of <var class="var">f</var> using adaptive
quadrature over the two-dimensional domain defined by <var class="var">xa</var>, <var class="var">xb</var>,
<var class="var">ya</var>, <var class="var">yb</var> using tiled integration. Additionally, <var class="var">ya</var> and
<var class="var">yb</var> may be scalar functions of <var class="var">x</var>, allowing for the integration
over non-rectangular domains.
</p>
<p><var class="var">f</var> is a function handle, inline function, or string containing the
name of the function to evaluate. The function <var class="var">f</var> must be of the form
<em class="math">z = f(x,y)</em>, and all operations must be vectorized such that <var class="var">x</var>
and <var class="var">y</var> accept array inputs and return array outputs of the same size.
(It can be assumed that <var class="var">x</var> and <var class="var">y</var> will either be same-size arrays
or one will be a scalar.) The underlying integrators will input arrays of
integration points into <var class="var">f</var> and/or use internal vector expansions to
speed computation that can produce unpredictable results if <var class="var">f</var> is not
restricted to elementwise operations. For integrands where this is
unavoidable, the <code class="code">("Vectorized") option described below may produce</code>
more reliable results.
</p>
<p>Additional optional parameters can be specified using
<code class="code">"<var class="var">property</var>", <var class="var">value</var></code> pairs. Valid properties are:
</p>
<dl class="table">
<dt><code class="code">AbsTol</code></dt>
<dd><p>Define the absolute error tolerance for the quadrature. The default
value is 1e-10 (1e-5 for single).
</p>
</dd>
<dt><code class="code">RelTol</code></dt>
<dd><p>Define the relative error tolerance for the quadrature. The default
value is 1e-6 (1e-4 for single).
</p>
</dd>
<dt><code class="code">MaxFunEvals</code></dt>
<dd><p>The maximum number of function calls to the vectorized function <var class="var">f</var>.
The default value is 5000.
</p>
</dd>
<dt><code class="code">Singular</code></dt>
<dd><p>Enable/disable transforms to weaken singularities on the edge of the
integration domain. The default value is <var class="var">true</var>.
</p>
</dd>
<dt><code class="code">Vectorized</code></dt>
<dd><p>Enable or disable vectorized integration. A value of <code class="code">false</code> forces
Octave to use only scalar inputs when calling the integrand, which enables
integrands <em class="math">f(x,y)</em> that have not been vectorized or only accept
scalar values of <var class="var">x</var> or <var class="var">y</var>. The default value is <code class="code">true</code>.
Note that this is achieved by wrapping <em class="math">f(x,y)</em> with the function
<code class="code">arrayfun</code>, which may significantly decrease computation speed.
</p>
</dd>
<dt><code class="code">FailurePlot</code></dt>
<dd><p>If <code class="code">quad2d</code> fails to converge to the desired error tolerance before
MaxFunEvals is reached, a plot of the areas that still need refinement
is created. The default value is <var class="var">false</var>.
</p></dd>
</dl>
<p>Adaptive quadrature is used to minimize the estimate of error until the
following is satisfied:
</p>
<div class="example">
<div class="group"><pre class="example-preformatted"> <var class="var">error</var> <= max (<var class="var">AbsTol</var>, <var class="var">RelTol</var>*|<var class="var">q</var>|)
</pre></div></div>
<p>The optional output <var class="var">err</var> is an approximate bound on the error in the
integral <code class="code">abs (<var class="var">q</var> - <var class="var">I</var>)</code>, where <var class="var">I</var> is the exact value
of the integral. The optional output <var class="var">iter</var> is the number of vectorized
function calls to the function <var class="var">f</var> that were used.
</p>
<p>Example 1 : integrate a rectangular region in x-y plane
</p>
<div class="example">
<div class="group"><pre class="example-preformatted"><var class="var">f</var> = @(<var class="var">x</var>,<var class="var">y</var>) 2*ones (size (<var class="var">x</var>));
<var class="var">q</var> = quad2d (<var class="var">f</var>, 0, 1, 0, 1)
⇒ <var class="var">q</var> = 2
</pre></div></div>
<p>The result is a volume, which for this constant-value integrand, is just
<code class="code"><var class="var">Length</var> * <var class="var">Width</var> * <var class="var">Height</var></code>.
</p>
<p>Example 2 : integrate a triangular region in x-y plane
</p>
<div class="example">
<div class="group"><pre class="example-preformatted"><var class="var">f</var> = @(<var class="var">x</var>,<var class="var">y</var>) 2*ones (size (<var class="var">x</var>));
<var class="var">ymax</var> = @(<var class="var">x</var>) 1 - <var class="var">x</var>;
<var class="var">q</var> = quad2d (<var class="var">f</var>, 0, 1, 0, <var class="var">ymax</var>)
⇒ <var class="var">q</var> = 1
</pre></div></div>
<p>The result is a volume, which for this constant-value integrand
<em class="math"><var class="var">f</var> = 2</em><!-- /@w -->, is the Triangle Area x Height or
<code class="code">1/2 * <var class="var">Base</var> * <var class="var">Width</var> * <var class="var">Height</var></code><!-- /@w -->.
</p>
<p>Example 3 : integrate a non-vectorized function over a square region
</p>
<div class="example">
<div class="group"><pre class="example-preformatted"><var class="var">f</var> = @(<var class="var">x</var>,<var class="var">y</var>) sinc (<var class="var">x</var>) * sinc (<var class="var">y</var>));
<var class="var">q</var> = quad2d (<var class="var">f</var>, -1, 1, -1, 1)
⇒ <var class="var">q</var> = 12.328 (incorrect)
<var class="var">q</var> = quad2d (<var class="var">f</var>, -1, 1, -1, 1, "Vectorized", false)
⇒ <var class="var">q</var> = 1.390 (correct)
<var class="var">f</var> = @(<var class="var">x</var>,<var class="var">y</var>) sinc (<var class="var">x</var>) .* sinc (<var class="var">y</var>);
<var class="var">q</var> = quad2d (<var class="var">f</var>, -1, 1, -1, 1)
⇒ <var class="var">q</var> = 1.390 (correct)
</pre></div></div>
<p>The first result is incorrect as the non-elementwise operator between the
sinc functions in <var class="var">f</var> create unintended matrix multiplications between
the internal integration arrays used by <code class="code">quad2d</code>. In the second
result, setting <code class="code">"Vectorized"</code> to false forces <code class="code">quad2d</code> to
perform scalar internal operations to compute the integral, resulting in
the correct numerical result at the cost of about a 20x increase in
computation time. In the third result, vectorizing the integrand <var class="var">f</var>
using the elementwise multiplication operator gets the correct result
without increasing computation time.
</p>
<p>Programming Notes: If there are singularities within the integration region
it is best to split the integral and place the singularities on the
boundary.
</p>
<p>Known <small class="sc">MATLAB</small> incompatibility: If tolerances are left unspecified, and
any integration limits are of type <code class="code">single</code>, then Octave’s integral
functions automatically reduce the default absolute and relative error
tolerances as specified above. If tighter tolerances are desired they
must be specified. <small class="sc">MATLAB</small> leaves the tighter tolerances appropriate
for <code class="code">double</code> inputs in place regardless of the class of the
integration limits.
</p>
<p>Reference: L.F. Shampine,
<cite class="cite"><small class="sc">MATLAB</small> program for quadrature in 2D</cite>, Applied Mathematics and
Computation, pp. 266–274, Vol 1, 2008.
</p>
<p><strong class="strong">See also:</strong> <a class="ref" href="#XREFintegral2">integral2</a>, <a class="ref" href="#XREFdblquad">dblquad</a>, <a class="ref" href="Functions-of-One-Variable.html#XREFintegral">integral</a>, <a class="ref" href="Functions-of-One-Variable.html#XREFquad">quad</a>, <a class="ref" href="Functions-of-One-Variable.html#XREFquadgk">quadgk</a>, <a class="ref" href="Functions-of-One-Variable.html#XREFquadv">quadv</a>, <a class="ref" href="Functions-of-One-Variable.html#XREFquadl">quadl</a>, <a class="ref" href="Functions-of-One-Variable.html#XREFquadcc">quadcc</a>, <a class="ref" href="Functions-of-One-Variable.html#XREFtrapz">trapz</a>, <a class="ref" href="#XREFintegral3">integral3</a>, <a class="ref" href="#XREFtriplequad">triplequad</a>.
</p></dd></dl>
<p>Finally, the functions <code class="code">integral2</code> and <code class="code">integral3</code> are provided
as general 2-D and 3-D integration functions. They will auto-select between
iterated and tiled integration methods and, unlike <code class="code">dblquad</code> and
<code class="code">triplequad</code>, will work with non-rectangular integration domains.
</p>
<a class="anchor" id="XREFintegral2"></a><span style="display:block; margin-top:-4.5ex;"> </span>
<dl class="first-deftypefn">
<dt class="deftypefn" id="index-integral2"><span><code class="def-type"><var class="var">q</var> =</code> <strong class="def-name">integral2</strong> <code class="def-code-arguments">(<var class="var">f</var>, <var class="var">xa</var>, <var class="var">xb</var>, <var class="var">ya</var>, <var class="var">yb</var>)</code><a class="copiable-link" href="#index-integral2"> ¶</a></span></dt>
<dt class="deftypefnx def-cmd-deftypefn" id="index-integral2-1"><span><code class="def-type"><var class="var">q</var> =</code> <strong class="def-name">integral2</strong> <code class="def-code-arguments">(<var class="var">f</var>, <var class="var">xa</var>, <var class="var">xb</var>, <var class="var">ya</var>, <var class="var">yb</var>, <var class="var">prop</var>, <var class="var">val</var>, …)</code><a class="copiable-link" href="#index-integral2-1"> ¶</a></span></dt>
<dt class="deftypefnx def-cmd-deftypefn" id="index-integral2-2"><span><code class="def-type">[<var class="var">q</var>, <var class="var">err</var>] =</code> <strong class="def-name">integral2</strong> <code class="def-code-arguments">(…)</code><a class="copiable-link" href="#index-integral2-2"> ¶</a></span></dt>
<dd>
<p>Numerically evaluate the two-dimensional integral of <var class="var">f</var> using adaptive
quadrature over the two-dimensional domain defined by <var class="var">xa</var>, <var class="var">xb</var>,
<var class="var">ya</var>, <var class="var">yb</var> (scalars may be finite or infinite). Additionally,
<var class="var">ya</var> and <var class="var">yb</var> may be scalar functions of <var class="var">x</var>, allowing for
integration over non-rectangular domains.
</p>
<p><var class="var">f</var> is a function handle, inline function, or string containing the
name of the function to evaluate. The function <var class="var">f</var> must be of the form
<em class="math">z = f(x,y)</em>, and all operations must be vectorized such that <var class="var">x</var>
and <var class="var">y</var> accept array inputs and return array outputs of the same size.
(It can be assumed that <var class="var">x</var> and <var class="var">y</var> will either be same-size arrays
or one will be a scalar.) The underlying integrators will input arrays of
integration points into <var class="var">f</var> and/or use internal vector expansions to
speed computation that can produce unpredictable results if <var class="var">f</var> is not
restricted to elementwise operations. For integrands where this is
unavoidable, the <code class="code">("Vectorized") option described below may produce</code>
more reliable results.
</p>
<p>Additional optional parameters can be specified using
<code class="code">"<var class="var">property</var>", <var class="var">value</var></code> pairs. Valid properties are:
</p>
<dl class="table">
<dt><code class="code">AbsTol</code></dt>
<dd><p>Define the absolute error tolerance for the quadrature. The default
value is 1e-10 (1e-5 for single).
</p>
</dd>
<dt><code class="code">RelTol</code></dt>
<dd><p>Define the relative error tolerance for the quadrature. The default
value is 1e-6 (1e-4 for single).
</p>
</dd>
<dt><code class="code">Method</code></dt>
<dd><p>Specify the two-dimensional integration method to be used, with valid
options being <code class="code">"auto"</code> (default), <code class="code">"tiled"</code>, or
<code class="code">"iterated"</code>. When using <code class="code">"auto"</code>, Octave will choose the
<code class="code">"tiled"</code> method unless any of the integration limits are infinite.
</p>
</dd>
<dt><code class="code">Vectorized</code></dt>
<dd><p>Enable or disable vectorized integration. A value of <code class="code">false</code> forces
Octave to use only scalar inputs when calling the integrand, which enables
integrands <em class="math">f(x,y)</em> that have not been vectorized or only accept
scalar values of <var class="var">x</var> or <var class="var">y</var>. The default value is <code class="code">true</code>.
Note that this is achieved by wrapping <em class="math">f(x,y)</em> with the function
<code class="code">arrayfun</code>, which may significantly decrease computation speed.
</p></dd>
</dl>
<p>Adaptive quadrature is used to minimize the estimate of error until the
following is satisfied:
</p>
<div class="example">
<div class="group"><pre class="example-preformatted"> <var class="var">error</var> <= max (<var class="var">AbsTol</var>, <var class="var">RelTol</var>*|<var class="var">q</var>|)
</pre></div></div>
<p><var class="var">err</var> is an approximate bound on the error in the integral
<code class="code">abs (<var class="var">q</var> - <var class="var">I</var>)</code>, where <var class="var">I</var> is the exact value of the
integral.
</p>
<p>Example 1 : integrate a rectangular region in x-y plane
</p>
<div class="example">
<div class="group"><pre class="example-preformatted"><var class="var">f</var> = @(<var class="var">x</var>,<var class="var">y</var>) 2*ones (size (<var class="var">x</var>));
<var class="var">q</var> = integral2 (<var class="var">f</var>, 0, 1, 0, 1)
⇒ <var class="var">q</var> = 2
</pre></div></div>
<p>The result is a volume, which for this constant-value integrand, is just
<code class="code"><var class="var">Length</var> * <var class="var">Width</var> * <var class="var">Height</var></code><!-- /@w -->.
</p>
<p>Example 2 : integrate a triangular region in x-y plane
</p>
<div class="example">
<div class="group"><pre class="example-preformatted"><var class="var">f</var> = @(<var class="var">x</var>,<var class="var">y</var>) 2*ones (size (<var class="var">x</var>));
<var class="var">ymax</var> = @(<var class="var">x</var>) 1 - <var class="var">x</var>;
<var class="var">q</var> = integral2 (<var class="var">f</var>, 0, 1, 0, <var class="var">ymax</var>)
⇒ <var class="var">q</var> = 1
</pre></div></div>
<p>The result is a volume, which for this constant-value integrand
<em class="math"><var class="var">f</var> = 2</em><!-- /@w -->, is the Triangle Area x Height or
<code class="code">1/2 * <var class="var">Base</var> * <var class="var">Width</var> * <var class="var">Height</var></code><!-- /@w -->.
</p>
<p>Example 3 : integrate a non-vectorized function over a square region
</p>
<div class="example">
<div class="group"><pre class="example-preformatted"><var class="var">f</var> = @(<var class="var">x</var>,<var class="var">y</var>) sinc (<var class="var">x</var>) * sinc (<var class="var">y</var>));
<var class="var">q</var> = integral2 (<var class="var">f</var>, -1, 1, -1, 1)
⇒ <var class="var">q</var> = 12.328 (incorrect)
<var class="var">q</var> = integral2 (<var class="var">f</var>, -1, 1, -1, 1, "Vectorized", false)
⇒ <var class="var">q</var> = 1.390 (correct)
<var class="var">f</var> = @(<var class="var">x</var>,<var class="var">y</var>) sinc (<var class="var">x</var>) .* sinc (<var class="var">y</var>);
<var class="var">q</var> = integral2 (<var class="var">f</var>, -1, 1, -1, 1)
⇒ <var class="var">q</var> = 1.390 (correct)
</pre></div></div>
<p>The first result is incorrect as the non-elementwise operator between the
sinc functions in <var class="var">f</var> create unintended matrix multiplications between
the internal integration arrays used by <code class="code">integral2</code>. In the second
result, setting <code class="code">"Vectorized"</code> to false forces <code class="code">integral2</code> to
perform scalar internal operations to compute the integral, resulting in
the correct numerical result at the cost of about a 20x increase in
computation time. In the third result, vectorizing the integrand <var class="var">f</var>
using the elementwise multiplication operator gets the correct result
without increasing computation time.
</p>
<p>Programming Notes: If there are singularities within the integration region
it is best to split the integral and place the singularities on the
boundary.
</p>
<p>Known <small class="sc">MATLAB</small> incompatibility: If tolerances are left unspecified, and
any integration limits are of type <code class="code">single</code>, then Octave’s integral
functions automatically reduce the default absolute and relative error
tolerances as specified above. If tighter tolerances are desired they
must be specified. <small class="sc">MATLAB</small> leaves the tighter tolerances appropriate
for <code class="code">double</code> inputs in place regardless of the class of the
integration limits.
</p>
<p>Reference: L.F. Shampine,
<cite class="cite"><small class="sc">MATLAB</small> program for quadrature in 2D</cite>, Applied Mathematics and
Computation, pp. 266–274, Vol 1, 2008.
</p>
<p><strong class="strong">See also:</strong> <a class="ref" href="#XREFquad2d">quad2d</a>, <a class="ref" href="#XREFdblquad">dblquad</a>, <a class="ref" href="Functions-of-One-Variable.html#XREFintegral">integral</a>, <a class="ref" href="Functions-of-One-Variable.html#XREFquad">quad</a>, <a class="ref" href="Functions-of-One-Variable.html#XREFquadgk">quadgk</a>, <a class="ref" href="Functions-of-One-Variable.html#XREFquadv">quadv</a>, <a class="ref" href="Functions-of-One-Variable.html#XREFquadl">quadl</a>, <a class="ref" href="Functions-of-One-Variable.html#XREFquadcc">quadcc</a>, <a class="ref" href="Functions-of-One-Variable.html#XREFtrapz">trapz</a>, <a class="ref" href="#XREFintegral3">integral3</a>, <a class="ref" href="#XREFtriplequad">triplequad</a>.
</p></dd></dl>
<a class="anchor" id="XREFintegral3"></a><span style="display:block; margin-top:-4.5ex;"> </span>
<dl class="first-deftypefn">
<dt class="deftypefn" id="index-integral3"><span><code class="def-type"><var class="var">q</var> =</code> <strong class="def-name">integral3</strong> <code class="def-code-arguments">(<var class="var">f</var>, <var class="var">xa</var>, <var class="var">xb</var>, <var class="var">ya</var>, <var class="var">yb</var>, <var class="var">za</var>, <var class="var">zb</var>)</code><a class="copiable-link" href="#index-integral3"> ¶</a></span></dt>
<dt class="deftypefnx def-cmd-deftypefn" id="index-integral3-1"><span><code class="def-type"><var class="var">q</var> =</code> <strong class="def-name">integral3</strong> <code class="def-code-arguments">(<var class="var">f</var>, <var class="var">xa</var>, <var class="var">xb</var>, <var class="var">ya</var>, <var class="var">yb</var>, <var class="var">za</var>, <var class="var">zb</var>, <var class="var">prop</var>, <var class="var">val</var>, …)</code><a class="copiable-link" href="#index-integral3-1"> ¶</a></span></dt>
<dd>
<p>Numerically evaluate the three-dimensional integral of <var class="var">f</var> using
adaptive quadrature over the three-dimensional domain defined by
<var class="var">xa</var>, <var class="var">xb</var>, <var class="var">ya</var>, <var class="var">yb</var>, <var class="var">za</var>, <var class="var">zb</var> (scalars may
be finite or infinite). Additionally, <var class="var">ya</var> and <var class="var">yb</var> may be
scalar functions of <var class="var">x</var> and <var class="var">za</var>, and <var class="var">zb</var> maybe be scalar
functions of <var class="var">x</var> and <var class="var">y</var>, allowing for integration over
non-rectangular domains.
</p>
<p><var class="var">f</var> is a function handle, inline function, or string containing the
name of the function to evaluate. The function <var class="var">f</var> must be of the form
<em class="math">z = f(x,y,z)</em>, and all operations must be vectorized such that
<var class="var">x</var>, <var class="var">y</var>, and <var class="var">z</var> accept array inputs and return array outputs
of the same size. (It can be assumed that <var class="var">x</var>, <var class="var">y</var>, and <var class="var">z</var>
will either be same-size arrays or scalars.) The underlying integrators
will input arrays of integration points into <var class="var">f</var> and/or use internal
vector expansions to speed computation that can produce unpredictable
results if <var class="var">f</var> is not restricted to elementwise operations. For
integrands where this is unavoidable, the <code class="code">("Vectorized") option</code>
described below may produce more reliable results.
</p>
<p>Additional optional parameters can be specified using
<code class="code">"<var class="var">property</var>", <var class="var">value</var></code> pairs. Valid properties are:
</p>
<dl class="table">
<dt><code class="code">AbsTol</code></dt>
<dd><p>Define the absolute error tolerance for the quadrature. The default
value is 1e-10 (1e-5 for single).
</p>
</dd>
<dt><code class="code">RelTol</code></dt>
<dd><p>Define the relative error tolerance for the quadrature. The default
value is 1e-6 (1e-4 for single).
</p>
</dd>
<dt><code class="code">Method</code></dt>
<dd><p>Specify the two-dimensional integration method to be used, with valid
options being <code class="code">"auto"</code> (default), <code class="code">"tiled"</code>, or
<code class="code">"iterated"</code>. When using <code class="code">"auto"</code>, Octave will choose the
<code class="code">"tiled"</code> method unless any of the integration limits are infinite.
</p>
</dd>
<dt><code class="code">Vectorized</code></dt>
<dd><p>Enable or disable vectorized integration. A value of <code class="code">false</code> forces
Octave to use only scalar inputs when calling the integrand, which enables
integrands <em class="math">f(x,y,z)</em> that have not been vectorized or only accept
scalar values of <var class="var">x</var>, <var class="var">y</var>, or <var class="var">z</var>. The default value is
<code class="code">true</code>. Note that this is achieved by wrapping <em class="math">f(x,y,z)</em> with
the function <code class="code">arrayfun</code>, which may significantly decrease computation
speed.
</p></dd>
</dl>
<p>Adaptive quadrature is used to minimize the estimate of error until the
following is satisfied:
</p>
<div class="example">
<div class="group"><pre class="example-preformatted"> <var class="var">error</var> <= max (<var class="var">AbsTol</var>, <var class="var">RelTol</var>*|<var class="var">q</var>|)
</pre></div></div>
<p><var class="var">err</var> is an approximate bound on the error in the integral
<code class="code">abs (<var class="var">q</var> - <var class="var">I</var>)</code>, where <var class="var">I</var> is the exact value of the
integral.
</p>
<p>Example 1 : integrate over a rectangular volume
</p>
<div class="example">
<div class="group"><pre class="example-preformatted"><var class="var">f</var> = @(<var class="var">x</var>,<var class="var">y</var>,<var class="var">z</var>) ones (size (<var class="var">x</var>));
<var class="var">q</var> = integral3 (<var class="var">f</var>, 0, 1, 0, 1, 0, 1)
⇒ <var class="var">q</var> = 1.00000
</pre></div></div>
<p>For this constant-value integrand, the result is a volume which is just
<code class="code"><var class="var">Length</var> * <var class="var">Width</var> * <var class="var">Height</var></code>.
</p>
<p>Example 2 : integrate over a spherical volume
</p>
<div class="example">
<div class="group"><pre class="example-preformatted"><var class="var">f</var> = @(<var class="var">x</var>,<var class="var">y</var>) ones (size (<var class="var">x</var>));
<var class="var">ymax</var> = @(<var class="var">x</var>) sqrt (1 - <var class="var">x</var>.^2);
<var class="var">zmax</var> = @(<var class="var">x</var>,<var class="var">y</var>) sqrt (1 - <var class="var">x</var>.^2 - <var class="var">y</var>.^2);
<var class="var">q</var> = integral3 (<var class="var">f</var>, 0, 1, 0, <var class="var">ymax</var>, 0, <var class="var">zmax</var>)
⇒ <var class="var">q</var> = 0.52360
</pre></div></div>
<p>For this constant-value integrand, the result is a volume which is 1/8th
of a unit sphere or <code class="code">1/8 * 4/3 * pi</code>.
</p>
<p>Example 3 : integrate a non-vectorized function over a cubic volume
</p>
<div class="example">
<div class="group"><pre class="example-preformatted"><var class="var">f</var> = @(<var class="var">x</var>,<var class="var">y</var>) sinc (<var class="var">x</var>) * sinc (<var class="var">y</var>), * sinc (<var class="var">z</var>);
<var class="var">q</var> = integral3 (<var class="var">f</var>, -1, 1, -1, 1, -1, 1)
⇒ <var class="var">q</var> = 14.535 (incorrect)
<var class="var">q</var> = integral3 (<var class="var">f</var>, -1, 1, -1, 1, -1, 1, "Vectorized", false)
⇒ <var class="var">q</var> = 1.6388 (correct)
<var class="var">f</var> = @(<var class="var">x</var>,<var class="var">y</var>,<var class="var">z</var>) sinc (<var class="var">x</var>) .* sinc (<var class="var">y</var>), .* sinc (<var class="var">z</var>);
<var class="var">q</var> = integral3 (<var class="var">f</var>, -1, 1, -1, 1, -1, 1)
⇒ <var class="var">q</var> = 1.6388 (correct)
</pre></div></div>
<p>The first result is incorrect as the non-elementwise operator between the
sinc functions in <var class="var">f</var> create unintended matrix multiplications between
the internal integration arrays used by <code class="code">integral3</code>. In the second
result, setting <code class="code">"Vectorized"</code> to false forces <code class="code">integral3</code> to
perform scalar internal operations to compute the integral, resulting in
the correct numerical result at the cost of about a 30x increase in
computation time. In the third result, vectorizing the integrand <var class="var">f</var>
using the elementwise multiplication operator gets the correct result
without increasing computation time.
</p>
<p>Programming Notes: If there are singularities within the integration region
it is best to split the integral and place the singularities on the
boundary.
</p>
<p>Known <small class="sc">MATLAB</small> incompatibility: If tolerances are left unspecified, and
any integration limits are of type <code class="code">single</code>, then Octave’s integral
functions automatically reduce the default absolute and relative error
tolerances as specified above. If tighter tolerances are desired they
must be specified. <small class="sc">MATLAB</small> leaves the tighter tolerances appropriate
for <code class="code">double</code> inputs in place regardless of the class of the
integration limits.
</p>
<p>Reference: L.F. Shampine,
<cite class="cite"><small class="sc">MATLAB</small> program for quadrature in 2D</cite>, Applied Mathematics and
Computation, pp. 266–274, Vol 1, 2008.
</p>
<p><strong class="strong">See also:</strong> <a class="ref" href="#XREFtriplequad">triplequad</a>, <a class="ref" href="Functions-of-One-Variable.html#XREFintegral">integral</a>, <a class="ref" href="Functions-of-One-Variable.html#XREFquad">quad</a>, <a class="ref" href="Functions-of-One-Variable.html#XREFquadgk">quadgk</a>, <a class="ref" href="Functions-of-One-Variable.html#XREFquadv">quadv</a>, <a class="ref" href="Functions-of-One-Variable.html#XREFquadl">quadl</a>, <a class="ref" href="Functions-of-One-Variable.html#XREFquadcc">quadcc</a>, <a class="ref" href="Functions-of-One-Variable.html#XREFtrapz">trapz</a>, <a class="ref" href="#XREFintegral2">integral2</a>, <a class="ref" href="#XREFquad2d">quad2d</a>, <a class="ref" href="#XREFdblquad">dblquad</a>.
</p></dd></dl>
<p>The above integrations can be fairly slow, and that problem increases
exponentially with the dimensionality of the integral. Another possible
solution for 2-D integration is to use Orthogonal Collocation as described in
the previous section (see <a class="pxref" href="Orthogonal-Collocation.html">Orthogonal Collocation</a>). The integral of a
function <em class="math">f(x,y)</em> for <em class="math">x</em> and <em class="math">y</em> between 0 and 1 can be
approximated using <em class="math">n</em> points by
the sum over <code class="code">i=1:n</code> and <code class="code">j=1:n</code> of <code class="code">q(i)*q(j)*f(r(i),r(j))</code>,
where <em class="math">q</em> and <em class="math">r</em> is as returned by <code class="code">colloc (n)</code>. The
generalization to more than two variables is straight forward. The
following code computes the studied integral using <em class="math">n=8</em> points.
</p>
<div class="example">
<div class="group"><pre class="example-preformatted">f = @(x,y) sin (pi*x*y') .* sqrt (x*y');
n = 8;
[t, ~, ~, q] = colloc (n);
I = q'*f(t,t)*q;
⇒ 0.30022
</pre></div></div>
<p>It should be noted that the number of points determines the quality
of the approximation. If the integration needs to be performed between
<em class="math">a</em> and <em class="math">b</em>, instead of 0 and 1, then a change of variables is needed.
</p>
</div>
<hr>
<div class="nav-panel">
<p>
Previous: <a href="Orthogonal-Collocation.html">Orthogonal Collocation</a>, Up: <a href="Numerical-Integration.html">Numerical Integration</a> [<a href="index.html#SEC_Contents" title="Table of contents" rel="contents">Contents</a>][<a href="Concept-Index.html" title="Index" rel="index">Index</a>]</p>
</div>
</body>
</html>
|
