Mercurial > octave-nkf

annotate doc/interpreter/quad.txi @ 12612:16cca721117b stable

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
doc: Update all documentation for chapter on Numerical Integration * cumtrapz.m, dblquad.m, quadgk.m, quadl.m, quadv.m, trapz.m, triplequad.m, quad.cc, quadcc.cc: Improve docstrings. * Quad-opts.in: Keep code sample together on a single line. * mk-opts.pl: Update quad-options function description * octave.texi: Update order of detailmenu to match order in quad.texi. * quad.txi: Add language about when to use each quad function, add examples of using trapz. * aspell-octave.en.pws: Add new spelling words from quad.texi chapter
author Rik <octave@nomad.inbox5.com>
date Sun, 17 Apr 2011 19:57:07 -0700
parents 1811f4f8b3b5
children 72c96de7a403
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
11523
fd0a3ac60b0e update copyright notices
John W. Eaton <jwe@octave.org>
parents: 10828
diff changeset
1 @c Copyright (C) 1996-2011 John W. Eaton
7018
fd42779a8428 [project @ 2007-10-13 00:52:12 by jwe]
jwe
parents: 6939
diff changeset
2 @c
fd42779a8428 [project @ 2007-10-13 00:52:12 by jwe]
jwe
parents: 6939
diff changeset
3 @c This file is part of Octave.
fd42779a8428 [project @ 2007-10-13 00:52:12 by jwe]
jwe
parents: 6939
diff changeset
4 @c
fd42779a8428 [project @ 2007-10-13 00:52:12 by jwe]
jwe
parents: 6939
diff changeset
5 @c Octave is free software; you can redistribute it and/or modify it
fd42779a8428 [project @ 2007-10-13 00:52:12 by jwe]
jwe
parents: 6939
diff changeset
6 @c under the terms of the GNU General Public License as published by the
fd42779a8428 [project @ 2007-10-13 00:52:12 by jwe]
jwe
parents: 6939
diff changeset
7 @c Free Software Foundation; either version 3 of the License, or (at
fd42779a8428 [project @ 2007-10-13 00:52:12 by jwe]
jwe
parents: 6939
diff changeset
8 @c your option) any later version.
fd42779a8428 [project @ 2007-10-13 00:52:12 by jwe]
jwe
parents: 6939
diff changeset
9 @c
fd42779a8428 [project @ 2007-10-13 00:52:12 by jwe]
jwe
parents: 6939
diff changeset
10 @c Octave is distributed in the hope that it will be useful, but WITHOUT
fd42779a8428 [project @ 2007-10-13 00:52:12 by jwe]
jwe
parents: 6939
diff changeset
11 @c ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
fd42779a8428 [project @ 2007-10-13 00:52:12 by jwe]
jwe
parents: 6939
diff changeset
12 @c FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
fd42779a8428 [project @ 2007-10-13 00:52:12 by jwe]
jwe
parents: 6939
diff changeset
13 @c for more details.
fd42779a8428 [project @ 2007-10-13 00:52:12 by jwe]
jwe
parents: 6939
diff changeset
14 @c
fd42779a8428 [project @ 2007-10-13 00:52:12 by jwe]
jwe
parents: 6939
diff changeset
15 @c You should have received a copy of the GNU General Public License
fd42779a8428 [project @ 2007-10-13 00:52:12 by jwe]
jwe
parents: 6939
diff changeset
16 @c along with Octave; see the file COPYING. If not, see
fd42779a8428 [project @ 2007-10-13 00:52:12 by jwe]
jwe
parents: 6939
diff changeset
17 @c <http://www.gnu.org/licenses/>.
3294
bfe1573bd2ae [project @ 1999-10-19 10:06:07 by jwe]
jwe
parents:
diff changeset
18
6741
00116015904d [project @ 2007-06-18 16:07:14 by jwe]
jwe
parents: 6502
diff changeset
19 @node Numerical Integration
00116015904d [project @ 2007-06-18 16:07:14 by jwe]
jwe
parents: 6502
diff changeset
20 @chapter Numerical Integration
00116015904d [project @ 2007-06-18 16:07:14 by jwe]
jwe
parents: 6502
diff changeset
21
00116015904d [project @ 2007-06-18 16:07:14 by jwe]
jwe
parents: 6502
diff changeset
22 Octave comes with several built-in functions for computing the integral
12612
16cca721117b doc: Update all documentation for chapter on Numerical Integration
Rik <octave@nomad.inbox5.com>
parents: 11562
diff changeset
23 of a function numerically (termed quadrature). These functions all solve
16cca721117b doc: Update all documentation for chapter on Numerical Integration
Rik <octave@nomad.inbox5.com>
parents: 11562
diff changeset
24 1-dimensional integration problems.
3294
bfe1573bd2ae [project @ 1999-10-19 10:06:07 by jwe]
jwe
parents:
diff changeset
25
bfe1573bd2ae [project @ 1999-10-19 10:06:07 by jwe]
jwe
parents:
diff changeset
26 @menu
6741
00116015904d [project @ 2007-06-18 16:07:14 by jwe]
jwe
parents: 6502
diff changeset
27 * Functions of One Variable::
12612
16cca721117b doc: Update all documentation for chapter on Numerical Integration
Rik <octave@nomad.inbox5.com>
parents: 11562
diff changeset
28 * Orthogonal Collocation::
6741
00116015904d [project @ 2007-06-18 16:07:14 by jwe]
jwe
parents: 6502
diff changeset
29 * Functions of Multiple Variables::
3294
bfe1573bd2ae [project @ 1999-10-19 10:06:07 by jwe]
jwe
parents:
diff changeset
30 @end menu
bfe1573bd2ae [project @ 1999-10-19 10:06:07 by jwe]
jwe
parents:
diff changeset
31
4167
aae05d51353c [project @ 2002-11-12 02:52:50 by jwe]
jwe
parents: 3368
diff changeset
32 @node Functions of One Variable
3294
bfe1573bd2ae [project @ 1999-10-19 10:06:07 by jwe]
jwe
parents:
diff changeset
33 @section Functions of One Variable
bfe1573bd2ae [project @ 1999-10-19 10:06:07 by jwe]
jwe
parents:
diff changeset
34
12612
16cca721117b doc: Update all documentation for chapter on Numerical Integration
Rik <octave@nomad.inbox5.com>
parents: 11562
diff changeset
35 Octave supports five different algorithms for computing the integral
6741
00116015904d [project @ 2007-06-18 16:07:14 by jwe]
jwe
parents: 6502
diff changeset
36 @tex
00116015904d [project @ 2007-06-18 16:07:14 by jwe]
jwe
parents: 6502
diff changeset
37 $$
00116015904d [project @ 2007-06-18 16:07:14 by jwe]
jwe
parents: 6502
diff changeset
38 \int_a^b f(x) d x
00116015904d [project @ 2007-06-18 16:07:14 by jwe]
jwe
parents: 6502
diff changeset
39 $$
00116015904d [project @ 2007-06-18 16:07:14 by jwe]
jwe
parents: 6502
diff changeset
40 @end tex
00116015904d [project @ 2007-06-18 16:07:14 by jwe]
jwe
parents: 6502
diff changeset
41 of a function @math{f} over the interval from @math{a} to @math{b}.
00116015904d [project @ 2007-06-18 16:07:14 by jwe]
jwe
parents: 6502
diff changeset
42 These are
00116015904d [project @ 2007-06-18 16:07:14 by jwe]
jwe
parents: 6502
diff changeset
43
00116015904d [project @ 2007-06-18 16:07:14 by jwe]
jwe
parents: 6502
diff changeset
44 @table @code
00116015904d [project @ 2007-06-18 16:07:14 by jwe]
jwe
parents: 6502
diff changeset
45 @item quad
00116015904d [project @ 2007-06-18 16:07:14 by jwe]
jwe
parents: 6502
diff changeset
46 Numerical integration based on Gaussian quadrature.
00116015904d [project @ 2007-06-18 16:07:14 by jwe]
jwe
parents: 6502
diff changeset
47
12612
16cca721117b doc: Update all documentation for chapter on Numerical Integration
Rik <octave@nomad.inbox5.com>
parents: 11562
diff changeset
48 @item quadv
16cca721117b doc: Update all documentation for chapter on Numerical Integration
Rik <octave@nomad.inbox5.com>
parents: 11562
diff changeset
49 Numerical integration using an adaptive vectorized Simpson's rule.
16cca721117b doc: Update all documentation for chapter on Numerical Integration
Rik <octave@nomad.inbox5.com>
parents: 11562
diff changeset
50
6741
00116015904d [project @ 2007-06-18 16:07:14 by jwe]
jwe
parents: 6502
diff changeset
51 @item quadl
00116015904d [project @ 2007-06-18 16:07:14 by jwe]
jwe
parents: 6502
diff changeset
52 Numerical integration using an adaptive Lobatto rule.
00116015904d [project @ 2007-06-18 16:07:14 by jwe]
jwe
parents: 6502
diff changeset
53
7984
bbaa5d7d0143 Some documentation updates
David Bateman <dbateman@free.fr>
parents: 7018
diff changeset
54 @item quadgk
9067
8970b4b10e9f Cleanup documentation for quad.texi and diffeq.texi
Rik <rdrider0-list@yahoo.com>
parents: 8920
diff changeset
55 Numerical integration using an adaptive Gauss-Konrod rule.
7984
bbaa5d7d0143 Some documentation updates
David Bateman <dbateman@free.fr>
parents: 7018
diff changeset
56
11562
1811f4f8b3b5 Rename cquad to quadcc and add to documentation.
Rik <octave@nomad.inbox5.com>
parents: 11523
diff changeset
57 @item quadcc
1811f4f8b3b5 Rename cquad to quadcc and add to documentation.
Rik <octave@nomad.inbox5.com>
parents: 11523
diff changeset
58 Numerical integration using adaptive Clenshaw-Curtis rules.
1811f4f8b3b5 Rename cquad to quadcc and add to documentation.
Rik <octave@nomad.inbox5.com>
parents: 11523
diff changeset
59
12612
16cca721117b doc: Update all documentation for chapter on Numerical Integration
Rik <octave@nomad.inbox5.com>
parents: 11562
diff changeset
60 @item trapz, cumtrapz
16cca721117b doc: Update all documentation for chapter on Numerical Integration
Rik <octave@nomad.inbox5.com>
parents: 11562
diff changeset
61 Numerical integration of data using the trapezoidal method.
6741
00116015904d [project @ 2007-06-18 16:07:14 by jwe]
jwe
parents: 6502
diff changeset
62 @end table
00116015904d [project @ 2007-06-18 16:07:14 by jwe]
jwe
parents: 6502
diff changeset
63
00116015904d [project @ 2007-06-18 16:07:14 by jwe]
jwe
parents: 6502
diff changeset
64 @noindent
12612
16cca721117b doc: Update all documentation for chapter on Numerical Integration
Rik <octave@nomad.inbox5.com>
parents: 11562
diff changeset
65 The best quadrature algorithm to use depends on the integrand. If you have
16cca721117b doc: Update all documentation for chapter on Numerical Integration
Rik <octave@nomad.inbox5.com>
parents: 11562
diff changeset
66 empirical data, rather than a function, the choice is @code{trapz} or
16cca721117b doc: Update all documentation for chapter on Numerical Integration
Rik <octave@nomad.inbox5.com>
parents: 11562
diff changeset
67 @code{cumtrapz}. If you are uncertain about the characteristics of the
16cca721117b doc: Update all documentation for chapter on Numerical Integration
Rik <octave@nomad.inbox5.com>
parents: 11562
diff changeset
68 integrand, @code{quadcc} will be the most robust as it can handle
16cca721117b doc: Update all documentation for chapter on Numerical Integration
Rik <octave@nomad.inbox5.com>
parents: 11562
diff changeset
69 discontinuities, singularities, oscillatory functions, and infinite intervals.
16cca721117b doc: Update all documentation for chapter on Numerical Integration
Rik <octave@nomad.inbox5.com>
parents: 11562
diff changeset
70 When the integrand is smooth @code{quadgk} may be the fastest of the
16cca721117b doc: Update all documentation for chapter on Numerical Integration
Rik <octave@nomad.inbox5.com>
parents: 11562
diff changeset
71 algorithms.
6741
00116015904d [project @ 2007-06-18 16:07:14 by jwe]
jwe
parents: 6502
diff changeset
72
12612
16cca721117b doc: Update all documentation for chapter on Numerical Integration
Rik <octave@nomad.inbox5.com>
parents: 11562
diff changeset
73 @multitable @columnfractions 0.05 0.15 .80
16cca721117b doc: Update all documentation for chapter on Numerical Integration
Rik <octave@nomad.inbox5.com>
parents: 11562
diff changeset
74 @headitem @tab Function @tab Characteristics
16cca721117b doc: Update all documentation for chapter on Numerical Integration
Rik <octave@nomad.inbox5.com>
parents: 11562
diff changeset
75 @item @tab quad @tab Low accuracy with nonsmooth integrands
16cca721117b doc: Update all documentation for chapter on Numerical Integration
Rik <octave@nomad.inbox5.com>
parents: 11562
diff changeset
76 @item @tab quadv @tab Medium accuracy with smooth integrands
16cca721117b doc: Update all documentation for chapter on Numerical Integration
Rik <octave@nomad.inbox5.com>
parents: 11562
diff changeset
77 @item @tab quadl @tab Medium accuracy with smooth integrands. Slower than quadgk.
16cca721117b doc: Update all documentation for chapter on Numerical Integration
Rik <octave@nomad.inbox5.com>
parents: 11562
diff changeset
78 @item @tab quadgk @tab Medium accuracy (@math{1e^{-6}}--@math{1e^{-9}}) with smooth integrands.
16cca721117b doc: Update all documentation for chapter on Numerical Integration
Rik <octave@nomad.inbox5.com>
parents: 11562
diff changeset
79 @item @tab @tab Handles oscillatory functions and infinite bounds
16cca721117b doc: Update all documentation for chapter on Numerical Integration
Rik <octave@nomad.inbox5.com>
parents: 11562
diff changeset
80 @item @tab quadcc @tab Low to High accuracy with nonsmooth/smooth integrands
16cca721117b doc: Update all documentation for chapter on Numerical Integration
Rik <octave@nomad.inbox5.com>
parents: 11562
diff changeset
81 @item @tab @tab Handles oscillatory functions, singularities, and infinite bounds
16cca721117b doc: Update all documentation for chapter on Numerical Integration
Rik <octave@nomad.inbox5.com>
parents: 11562
diff changeset
82 @end multitable
3294
bfe1573bd2ae [project @ 1999-10-19 10:06:07 by jwe]
jwe
parents:
diff changeset
83
bfe1573bd2ae [project @ 1999-10-19 10:06:07 by jwe]
jwe
parents:
diff changeset
84
bfe1573bd2ae [project @ 1999-10-19 10:06:07 by jwe]
jwe
parents:
diff changeset
85 Here is an example of using @code{quad} to integrate the function
bfe1573bd2ae [project @ 1999-10-19 10:06:07 by jwe]
jwe
parents:
diff changeset
86 @tex
bfe1573bd2ae [project @ 1999-10-19 10:06:07 by jwe]
jwe
parents:
diff changeset
87 $$
bfe1573bd2ae [project @ 1999-10-19 10:06:07 by jwe]
jwe
parents:
diff changeset
88 f(x) = x \sin (1/x) \sqrt {|1 - x|}
bfe1573bd2ae [project @ 1999-10-19 10:06:07 by jwe]
jwe
parents:
diff changeset
89 $$
bfe1573bd2ae [project @ 1999-10-19 10:06:07 by jwe]
jwe
parents:
diff changeset
90 from $x = 0$ to $x = 3$.
bfe1573bd2ae [project @ 1999-10-19 10:06:07 by jwe]
jwe
parents:
diff changeset
91 @end tex
6741
00116015904d [project @ 2007-06-18 16:07:14 by jwe]
jwe
parents: 6502
diff changeset
92 @ifnottex
3294
bfe1573bd2ae [project @ 1999-10-19 10:06:07 by jwe]
jwe
parents:
diff changeset
93
bfe1573bd2ae [project @ 1999-10-19 10:06:07 by jwe]
jwe
parents:
diff changeset
94 @example
bfe1573bd2ae [project @ 1999-10-19 10:06:07 by jwe]
jwe
parents:
diff changeset
95 @var{f}(@var{x}) = @var{x} * sin (1/@var{x}) * sqrt (abs (1 - @var{x}))
bfe1573bd2ae [project @ 1999-10-19 10:06:07 by jwe]
jwe
parents:
diff changeset
96 @end example
bfe1573bd2ae [project @ 1999-10-19 10:06:07 by jwe]
jwe
parents:
diff changeset
97
bfe1573bd2ae [project @ 1999-10-19 10:06:07 by jwe]
jwe
parents:
diff changeset
98 @noindent
bfe1573bd2ae [project @ 1999-10-19 10:06:07 by jwe]
jwe
parents:
diff changeset
99 from @var{x} = 0 to @var{x} = 3.
6741
00116015904d [project @ 2007-06-18 16:07:14 by jwe]
jwe
parents: 6502
diff changeset
100 @end ifnottex
3294
bfe1573bd2ae [project @ 1999-10-19 10:06:07 by jwe]
jwe
parents:
diff changeset
101
bfe1573bd2ae [project @ 1999-10-19 10:06:07 by jwe]
jwe
parents:
diff changeset
102 This is a fairly difficult integration (plot the function over the range
bfe1573bd2ae [project @ 1999-10-19 10:06:07 by jwe]
jwe
parents:
diff changeset
103 of integration to see why).
bfe1573bd2ae [project @ 1999-10-19 10:06:07 by jwe]
jwe
parents:
diff changeset
104
bfe1573bd2ae [project @ 1999-10-19 10:06:07 by jwe]
jwe
parents:
diff changeset
105 The first step is to define the function:
bfe1573bd2ae [project @ 1999-10-19 10:06:07 by jwe]
jwe
parents:
diff changeset
106
bfe1573bd2ae [project @ 1999-10-19 10:06:07 by jwe]
jwe
parents:
diff changeset
107 @example
bfe1573bd2ae [project @ 1999-10-19 10:06:07 by jwe]
jwe
parents:
diff changeset
108 @group
bfe1573bd2ae [project @ 1999-10-19 10:06:07 by jwe]
jwe
parents:
diff changeset
109 function y = f (x)
12612
16cca721117b doc: Update all documentation for chapter on Numerical Integration
Rik <octave@nomad.inbox5.com>
parents: 11562
diff changeset
110 y = x .* sin (1./x) .* sqrt (abs (1 - x));
3294
bfe1573bd2ae [project @ 1999-10-19 10:06:07 by jwe]
jwe
parents:
diff changeset
111 endfunction
bfe1573bd2ae [project @ 1999-10-19 10:06:07 by jwe]
jwe
parents:
diff changeset
112 @end group
bfe1573bd2ae [project @ 1999-10-19 10:06:07 by jwe]
jwe
parents:
diff changeset
113 @end example
bfe1573bd2ae [project @ 1999-10-19 10:06:07 by jwe]
jwe
parents:
diff changeset
114
12612
16cca721117b doc: Update all documentation for chapter on Numerical Integration
Rik <octave@nomad.inbox5.com>
parents: 11562
diff changeset
115 Note the use of the `dot' forms of the operators. This is not necessary for
16cca721117b doc: Update all documentation for chapter on Numerical Integration
Rik <octave@nomad.inbox5.com>
parents: 11562
diff changeset
116 the @code{quad} integrator, but is required by the other integrators. In any
16cca721117b doc: Update all documentation for chapter on Numerical Integration
Rik <octave@nomad.inbox5.com>
parents: 11562
diff changeset
117 case, it makes it much easier to generate a set of points for plotting because
16cca721117b doc: Update all documentation for chapter on Numerical Integration
Rik <octave@nomad.inbox5.com>
parents: 11562
diff changeset
118 it is possible to call the function with a vector argument to produce a vector
16cca721117b doc: Update all documentation for chapter on Numerical Integration
Rik <octave@nomad.inbox5.com>
parents: 11562
diff changeset
119 result.
3294
bfe1573bd2ae [project @ 1999-10-19 10:06:07 by jwe]
jwe
parents:
diff changeset
120
12612
16cca721117b doc: Update all documentation for chapter on Numerical Integration
Rik <octave@nomad.inbox5.com>
parents: 11562
diff changeset
121 The second step is to call quad with the limits of integration:
3294
bfe1573bd2ae [project @ 1999-10-19 10:06:07 by jwe]
jwe
parents:
diff changeset
122
bfe1573bd2ae [project @ 1999-10-19 10:06:07 by jwe]
jwe
parents:
diff changeset
123 @example
bfe1573bd2ae [project @ 1999-10-19 10:06:07 by jwe]
jwe
parents:
diff changeset
124 @group
12612
16cca721117b doc: Update all documentation for chapter on Numerical Integration
Rik <octave@nomad.inbox5.com>
parents: 11562
diff changeset
125 [q, ier, nfun, err] = quad ("f", 0, 3)
3294
bfe1573bd2ae [project @ 1999-10-19 10:06:07 by jwe]
jwe
parents:
diff changeset
126 @result{} 1.9819
bfe1573bd2ae [project @ 1999-10-19 10:06:07 by jwe]
jwe
parents:
diff changeset
127 @result{} 1
bfe1573bd2ae [project @ 1999-10-19 10:06:07 by jwe]
jwe
parents:
diff changeset
128 @result{} 5061
bfe1573bd2ae [project @ 1999-10-19 10:06:07 by jwe]
jwe
parents:
diff changeset
129 @result{} 1.1522e-07
bfe1573bd2ae [project @ 1999-10-19 10:06:07 by jwe]
jwe
parents:
diff changeset
130 @end group
bfe1573bd2ae [project @ 1999-10-19 10:06:07 by jwe]
jwe
parents:
diff changeset
131 @end example
bfe1573bd2ae [project @ 1999-10-19 10:06:07 by jwe]
jwe
parents:
diff changeset
132
bfe1573bd2ae [project @ 1999-10-19 10:06:07 by jwe]
jwe
parents:
diff changeset
133 Although @code{quad} returns a nonzero value for @var{ier}, the result
bfe1573bd2ae [project @ 1999-10-19 10:06:07 by jwe]
jwe
parents:
diff changeset
134 is reasonably accurate (to see why, examine what happens to the result
bfe1573bd2ae [project @ 1999-10-19 10:06:07 by jwe]
jwe
parents:
diff changeset
135 if you move the lower bound to 0.1, then 0.01, then 0.001, etc.).
bfe1573bd2ae [project @ 1999-10-19 10:06:07 by jwe]
jwe
parents:
diff changeset
136
12612
16cca721117b doc: Update all documentation for chapter on Numerical Integration
Rik <octave@nomad.inbox5.com>
parents: 11562
diff changeset
137 The function "f" can be the string name of a function, a function handle, or
16cca721117b doc: Update all documentation for chapter on Numerical Integration
Rik <octave@nomad.inbox5.com>
parents: 11562
diff changeset
138 an inline function. These options make it quite easy to do integration
16cca721117b doc: Update all documentation for chapter on Numerical Integration
Rik <octave@nomad.inbox5.com>
parents: 11562
diff changeset
139 without having to fully define a function in an m-file. For example:
16cca721117b doc: Update all documentation for chapter on Numerical Integration
Rik <octave@nomad.inbox5.com>
parents: 11562
diff changeset
140
16cca721117b doc: Update all documentation for chapter on Numerical Integration
Rik <octave@nomad.inbox5.com>
parents: 11562
diff changeset
141 @example
16cca721117b doc: Update all documentation for chapter on Numerical Integration
Rik <octave@nomad.inbox5.com>
parents: 11562
diff changeset
142 @group
16cca721117b doc: Update all documentation for chapter on Numerical Integration
Rik <octave@nomad.inbox5.com>
parents: 11562
diff changeset
143 # Verify integral (x^3) = x^4/4
16cca721117b doc: Update all documentation for chapter on Numerical Integration
Rik <octave@nomad.inbox5.com>
parents: 11562
diff changeset
144 f = inline ("x.^3");
16cca721117b doc: Update all documentation for chapter on Numerical Integration
Rik <octave@nomad.inbox5.com>
parents: 11562
diff changeset
145 quadgk (f, 0, 1)
16cca721117b doc: Update all documentation for chapter on Numerical Integration
Rik <octave@nomad.inbox5.com>
parents: 11562
diff changeset
146 @result{} 0.25000
16cca721117b doc: Update all documentation for chapter on Numerical Integration
Rik <octave@nomad.inbox5.com>
parents: 11562
diff changeset
147
16cca721117b doc: Update all documentation for chapter on Numerical Integration
Rik <octave@nomad.inbox5.com>
parents: 11562
diff changeset
148 # Verify gamma function = (n-1)! for n = 4
16cca721117b doc: Update all documentation for chapter on Numerical Integration
Rik <octave@nomad.inbox5.com>
parents: 11562
diff changeset
149 f = @@(x) x.^3 .* exp (-x);
16cca721117b doc: Update all documentation for chapter on Numerical Integration
Rik <octave@nomad.inbox5.com>
parents: 11562
diff changeset
150 quadcc (f, 0, Inf)
16cca721117b doc: Update all documentation for chapter on Numerical Integration
Rik <octave@nomad.inbox5.com>
parents: 11562
diff changeset
151 @result{} 6.0000
16cca721117b doc: Update all documentation for chapter on Numerical Integration
Rik <octave@nomad.inbox5.com>
parents: 11562
diff changeset
152 @end group
16cca721117b doc: Update all documentation for chapter on Numerical Integration
Rik <octave@nomad.inbox5.com>
parents: 11562
diff changeset
153 @end example
16cca721117b doc: Update all documentation for chapter on Numerical Integration
Rik <octave@nomad.inbox5.com>
parents: 11562
diff changeset
154
16cca721117b doc: Update all documentation for chapter on Numerical Integration
Rik <octave@nomad.inbox5.com>
parents: 11562
diff changeset
155 @DOCSTRING(quad)
16cca721117b doc: Update all documentation for chapter on Numerical Integration
Rik <octave@nomad.inbox5.com>
parents: 11562
diff changeset
156
16cca721117b doc: Update all documentation for chapter on Numerical Integration
Rik <octave@nomad.inbox5.com>
parents: 11562
diff changeset
157 @DOCSTRING(quad_options)
16cca721117b doc: Update all documentation for chapter on Numerical Integration
Rik <octave@nomad.inbox5.com>
parents: 11562
diff changeset
158
16cca721117b doc: Update all documentation for chapter on Numerical Integration
Rik <octave@nomad.inbox5.com>
parents: 11562
diff changeset
159 @DOCSTRING(quadv)
16cca721117b doc: Update all documentation for chapter on Numerical Integration
Rik <octave@nomad.inbox5.com>
parents: 11562
diff changeset
160
6502
6ab0a8767780 [project @ 2007-04-06 03:32:06 by jwe]
jwe
parents: 4167
diff changeset
161 @DOCSTRING(quadl)
6ab0a8767780 [project @ 2007-04-06 03:32:06 by jwe]
jwe
parents: 4167
diff changeset
162
7984
bbaa5d7d0143 Some documentation updates
David Bateman <dbateman@free.fr>
parents: 7018
diff changeset
163 @DOCSTRING(quadgk)
bbaa5d7d0143 Some documentation updates
David Bateman <dbateman@free.fr>
parents: 7018
diff changeset
164
12612
16cca721117b doc: Update all documentation for chapter on Numerical Integration
Rik <octave@nomad.inbox5.com>
parents: 11562
diff changeset
165 @DOCSTRING(quadcc)
16cca721117b doc: Update all documentation for chapter on Numerical Integration
Rik <octave@nomad.inbox5.com>
parents: 11562
diff changeset
166
16cca721117b doc: Update all documentation for chapter on Numerical Integration
Rik <octave@nomad.inbox5.com>
parents: 11562
diff changeset
167 Sometimes one does not have the function, but only the raw (x, y) points from
16cca721117b doc: Update all documentation for chapter on Numerical Integration
Rik <octave@nomad.inbox5.com>
parents: 11562
diff changeset
168 which to perform an integration. This can occur when collecting data in an
16cca721117b doc: Update all documentation for chapter on Numerical Integration
Rik <octave@nomad.inbox5.com>
parents: 11562
diff changeset
169 experiment. The @code{trapz} function can integrate these values as shown in
16cca721117b doc: Update all documentation for chapter on Numerical Integration
Rik <octave@nomad.inbox5.com>
parents: 11562
diff changeset
170 the following example where "data" has been collected on the cosine function
16cca721117b doc: Update all documentation for chapter on Numerical Integration
Rik <octave@nomad.inbox5.com>
parents: 11562
diff changeset
171 over the range [0, pi/2).
16cca721117b doc: Update all documentation for chapter on Numerical Integration
Rik <octave@nomad.inbox5.com>
parents: 11562
diff changeset
172
16cca721117b doc: Update all documentation for chapter on Numerical Integration
Rik <octave@nomad.inbox5.com>
parents: 11562
diff changeset
173 @example
16cca721117b doc: Update all documentation for chapter on Numerical Integration
Rik <octave@nomad.inbox5.com>
parents: 11562
diff changeset
174 @group
16cca721117b doc: Update all documentation for chapter on Numerical Integration
Rik <octave@nomad.inbox5.com>
parents: 11562
diff changeset
175 x = 0:0.1:pi/2; # Uniformly spaced points
16cca721117b doc: Update all documentation for chapter on Numerical Integration
Rik <octave@nomad.inbox5.com>
parents: 11562
diff changeset
176 y = cos (x);
16cca721117b doc: Update all documentation for chapter on Numerical Integration
Rik <octave@nomad.inbox5.com>
parents: 11562
diff changeset
177 trapz (x, y)
16cca721117b doc: Update all documentation for chapter on Numerical Integration
Rik <octave@nomad.inbox5.com>
parents: 11562
diff changeset
178 @result{} 0.99666
16cca721117b doc: Update all documentation for chapter on Numerical Integration
Rik <octave@nomad.inbox5.com>
parents: 11562
diff changeset
179 @end group
16cca721117b doc: Update all documentation for chapter on Numerical Integration
Rik <octave@nomad.inbox5.com>
parents: 11562
diff changeset
180 @end example
16cca721117b doc: Update all documentation for chapter on Numerical Integration
Rik <octave@nomad.inbox5.com>
parents: 11562
diff changeset
181
16cca721117b doc: Update all documentation for chapter on Numerical Integration
Rik <octave@nomad.inbox5.com>
parents: 11562
diff changeset
182 The answer is reasonably close to the exact value of 1. Ordinary quadrature
16cca721117b doc: Update all documentation for chapter on Numerical Integration
Rik <octave@nomad.inbox5.com>
parents: 11562
diff changeset
183 is sensitive to the characteristics of the integrand. Empirical integration
16cca721117b doc: Update all documentation for chapter on Numerical Integration
Rik <octave@nomad.inbox5.com>
parents: 11562
diff changeset
184 depends not just on the integrand, but also on the particular points chosen to
16cca721117b doc: Update all documentation for chapter on Numerical Integration
Rik <octave@nomad.inbox5.com>
parents: 11562
diff changeset
185 represent the function. Repeating the example above with the sine function
16cca721117b doc: Update all documentation for chapter on Numerical Integration
Rik <octave@nomad.inbox5.com>
parents: 11562
diff changeset
186 over the range [0, pi/2) produces far inferior results.
7984
bbaa5d7d0143 Some documentation updates
David Bateman <dbateman@free.fr>
parents: 7018
diff changeset
187
12612
16cca721117b doc: Update all documentation for chapter on Numerical Integration
Rik <octave@nomad.inbox5.com>
parents: 11562
diff changeset
188 @example
16cca721117b doc: Update all documentation for chapter on Numerical Integration
Rik <octave@nomad.inbox5.com>
parents: 11562
diff changeset
189 @group
16cca721117b doc: Update all documentation for chapter on Numerical Integration
Rik <octave@nomad.inbox5.com>
parents: 11562
diff changeset
190 x = 0:0.1:pi/2; # Uniformly spaced points
16cca721117b doc: Update all documentation for chapter on Numerical Integration
Rik <octave@nomad.inbox5.com>
parents: 11562
diff changeset
191 y = sin (x);
16cca721117b doc: Update all documentation for chapter on Numerical Integration
Rik <octave@nomad.inbox5.com>
parents: 11562
diff changeset
192 trapz (x, y)
16cca721117b doc: Update all documentation for chapter on Numerical Integration
Rik <octave@nomad.inbox5.com>
parents: 11562
diff changeset
193 @result{} 0.92849
16cca721117b doc: Update all documentation for chapter on Numerical Integration
Rik <octave@nomad.inbox5.com>
parents: 11562
diff changeset
194 @end group
16cca721117b doc: Update all documentation for chapter on Numerical Integration
Rik <octave@nomad.inbox5.com>
parents: 11562
diff changeset
195 @end example
16cca721117b doc: Update all documentation for chapter on Numerical Integration
Rik <octave@nomad.inbox5.com>
parents: 11562
diff changeset
196
16cca721117b doc: Update all documentation for chapter on Numerical Integration
Rik <octave@nomad.inbox5.com>
parents: 11562
diff changeset
197 However, a slightly different choice of data points can change the result
16cca721117b doc: Update all documentation for chapter on Numerical Integration
Rik <octave@nomad.inbox5.com>
parents: 11562
diff changeset
198 significantly. The same integration, with the same number of points, but
16cca721117b doc: Update all documentation for chapter on Numerical Integration
Rik <octave@nomad.inbox5.com>
parents: 11562
diff changeset
199 spaced differently produces a more accurate answer.
16cca721117b doc: Update all documentation for chapter on Numerical Integration
Rik <octave@nomad.inbox5.com>
parents: 11562
diff changeset
200
16cca721117b doc: Update all documentation for chapter on Numerical Integration
Rik <octave@nomad.inbox5.com>
parents: 11562
diff changeset
201 @example
16cca721117b doc: Update all documentation for chapter on Numerical Integration
Rik <octave@nomad.inbox5.com>
parents: 11562
diff changeset
202 @group
16cca721117b doc: Update all documentation for chapter on Numerical Integration
Rik <octave@nomad.inbox5.com>
parents: 11562
diff changeset
203 x = linspace (0, pi/2, 16); # Uniformly spaced, but including endpoint
16cca721117b doc: Update all documentation for chapter on Numerical Integration
Rik <octave@nomad.inbox5.com>
parents: 11562
diff changeset
204 y = sin (x);
16cca721117b doc: Update all documentation for chapter on Numerical Integration
Rik <octave@nomad.inbox5.com>
parents: 11562
diff changeset
205 trapz (x, y)
16cca721117b doc: Update all documentation for chapter on Numerical Integration
Rik <octave@nomad.inbox5.com>
parents: 11562
diff changeset
206 @result{} 0.99909
16cca721117b doc: Update all documentation for chapter on Numerical Integration
Rik <octave@nomad.inbox5.com>
parents: 11562
diff changeset
207 @end group
16cca721117b doc: Update all documentation for chapter on Numerical Integration
Rik <octave@nomad.inbox5.com>
parents: 11562
diff changeset
208 @end example
16cca721117b doc: Update all documentation for chapter on Numerical Integration
Rik <octave@nomad.inbox5.com>
parents: 11562
diff changeset
209
16cca721117b doc: Update all documentation for chapter on Numerical Integration
Rik <octave@nomad.inbox5.com>
parents: 11562
diff changeset
210 In general there may be no way of knowing the best distribution of points ahead
16cca721117b doc: Update all documentation for chapter on Numerical Integration
Rik <octave@nomad.inbox5.com>
parents: 11562
diff changeset
211 of time. Or the points may come from an experiment where there is no freedom to
16cca721117b doc: Update all documentation for chapter on Numerical Integration
Rik <octave@nomad.inbox5.com>
parents: 11562
diff changeset
212 select the best distribution. In any case, one must remain aware of this
16cca721117b doc: Update all documentation for chapter on Numerical Integration
Rik <octave@nomad.inbox5.com>
parents: 11562
diff changeset
213 issue when using @code{trapz}.
11562
1811f4f8b3b5 Rename cquad to quadcc and add to documentation.
Rik <octave@nomad.inbox5.com>
parents: 11523
diff changeset
214
6502
6ab0a8767780 [project @ 2007-04-06 03:32:06 by jwe]
jwe
parents: 4167
diff changeset
215 @DOCSTRING(trapz)
6ab0a8767780 [project @ 2007-04-06 03:32:06 by jwe]
jwe
parents: 4167
diff changeset
216
6ab0a8767780 [project @ 2007-04-06 03:32:06 by jwe]
jwe
parents: 4167
diff changeset
217 @DOCSTRING(cumtrapz)
6ab0a8767780 [project @ 2007-04-06 03:32:06 by jwe]
jwe
parents: 4167
diff changeset
218
4167
aae05d51353c [project @ 2002-11-12 02:52:50 by jwe]
jwe
parents: 3368
diff changeset
219 @node Orthogonal Collocation
3294
bfe1573bd2ae [project @ 1999-10-19 10:06:07 by jwe]
jwe
parents:
diff changeset
220 @section Orthogonal Collocation
bfe1573bd2ae [project @ 1999-10-19 10:06:07 by jwe]
jwe
parents:
diff changeset
221
3368
a4cd1e9d9962 [project @ 1999-11-20 17:22:48 by jwe]
jwe
parents: 3294
diff changeset
222 @DOCSTRING(colloc)
3294
bfe1573bd2ae [project @ 1999-10-19 10:06:07 by jwe]
jwe
parents:
diff changeset
223
bfe1573bd2ae [project @ 1999-10-19 10:06:07 by jwe]
jwe
parents:
diff changeset
224 Here is an example of using @code{colloc} to generate weight matrices
bfe1573bd2ae [project @ 1999-10-19 10:06:07 by jwe]
jwe
parents:
diff changeset
225 for solving the second order differential equation
bfe1573bd2ae [project @ 1999-10-19 10:06:07 by jwe]
jwe
parents:
diff changeset
226 @tex
bfe1573bd2ae [project @ 1999-10-19 10:06:07 by jwe]
jwe
parents:
diff changeset
227 $u^\prime - \alpha u^{\prime\prime} = 0$ with the boundary conditions
bfe1573bd2ae [project @ 1999-10-19 10:06:07 by jwe]
jwe
parents:
diff changeset
228 $u(0) = 0$ and $u(1) = 1$.
bfe1573bd2ae [project @ 1999-10-19 10:06:07 by jwe]
jwe
parents:
diff changeset
229 @end tex
6741
00116015904d [project @ 2007-06-18 16:07:14 by jwe]
jwe
parents: 6502
diff changeset
230 @ifnottex
3294
bfe1573bd2ae [project @ 1999-10-19 10:06:07 by jwe]
jwe
parents:
diff changeset
231 @var{u}' - @var{alpha} * @var{u}'' = 0 with the boundary conditions
bfe1573bd2ae [project @ 1999-10-19 10:06:07 by jwe]
jwe
parents:
diff changeset
232 @var{u}(0) = 0 and @var{u}(1) = 1.
6741
00116015904d [project @ 2007-06-18 16:07:14 by jwe]
jwe
parents: 6502
diff changeset
233 @end ifnottex
3294
bfe1573bd2ae [project @ 1999-10-19 10:06:07 by jwe]
jwe
parents:
diff changeset
234
bfe1573bd2ae [project @ 1999-10-19 10:06:07 by jwe]
jwe
parents:
diff changeset
235 First, we can generate the weight matrices for @var{n} points (including
bfe1573bd2ae [project @ 1999-10-19 10:06:07 by jwe]
jwe
parents:
diff changeset
236 the endpoints of the interval), and incorporate the boundary conditions
bfe1573bd2ae [project @ 1999-10-19 10:06:07 by jwe]
jwe
parents:
diff changeset
237 in the right hand side (for a specific value of
bfe1573bd2ae [project @ 1999-10-19 10:06:07 by jwe]
jwe
parents:
diff changeset
238 @tex
bfe1573bd2ae [project @ 1999-10-19 10:06:07 by jwe]
jwe
parents:
diff changeset
239 $\alpha$).
bfe1573bd2ae [project @ 1999-10-19 10:06:07 by jwe]
jwe
parents:
diff changeset
240 @end tex
6741
00116015904d [project @ 2007-06-18 16:07:14 by jwe]
jwe
parents: 6502
diff changeset
241 @ifnottex
3294
bfe1573bd2ae [project @ 1999-10-19 10:06:07 by jwe]
jwe
parents:
diff changeset
242 @var{alpha}).
6741
00116015904d [project @ 2007-06-18 16:07:14 by jwe]
jwe
parents: 6502
diff changeset
243 @end ifnottex
3294
bfe1573bd2ae [project @ 1999-10-19 10:06:07 by jwe]
jwe
parents:
diff changeset
244
bfe1573bd2ae [project @ 1999-10-19 10:06:07 by jwe]
jwe
parents:
diff changeset
245 @example
bfe1573bd2ae [project @ 1999-10-19 10:06:07 by jwe]
jwe
parents:
diff changeset
246 @group
bfe1573bd2ae [project @ 1999-10-19 10:06:07 by jwe]
jwe
parents:
diff changeset
247 n = 7;
bfe1573bd2ae [project @ 1999-10-19 10:06:07 by jwe]
jwe
parents:
diff changeset
248 alpha = 0.1;
bfe1573bd2ae [project @ 1999-10-19 10:06:07 by jwe]
jwe
parents:
diff changeset
249 [r, a, b] = colloc (n-2, "left", "right");
bfe1573bd2ae [project @ 1999-10-19 10:06:07 by jwe]
jwe
parents:
diff changeset
250 at = a(2:n-1,2:n-1);
bfe1573bd2ae [project @ 1999-10-19 10:06:07 by jwe]
jwe
parents:
diff changeset
251 bt = b(2:n-1,2:n-1);
bfe1573bd2ae [project @ 1999-10-19 10:06:07 by jwe]
jwe
parents:
diff changeset
252 rhs = alpha * b(2:n-1,n) - a(2:n-1,n);
bfe1573bd2ae [project @ 1999-10-19 10:06:07 by jwe]
jwe
parents:
diff changeset
253 @end group
bfe1573bd2ae [project @ 1999-10-19 10:06:07 by jwe]
jwe
parents:
diff changeset
254 @end example
bfe1573bd2ae [project @ 1999-10-19 10:06:07 by jwe]
jwe
parents:
diff changeset
255
bfe1573bd2ae [project @ 1999-10-19 10:06:07 by jwe]
jwe
parents:
diff changeset
256 Then the solution at the roots @var{r} is
bfe1573bd2ae [project @ 1999-10-19 10:06:07 by jwe]
jwe
parents:
diff changeset
257
bfe1573bd2ae [project @ 1999-10-19 10:06:07 by jwe]
jwe
parents:
diff changeset
258 @example
9067
8970b4b10e9f Cleanup documentation for quad.texi and diffeq.texi
Rik <rdrider0-list@yahoo.com>
parents: 8920
diff changeset
259 @group
3294
bfe1573bd2ae [project @ 1999-10-19 10:06:07 by jwe]
jwe
parents:
diff changeset
260 u = [ 0; (at - alpha * bt) \ rhs; 1]
bfe1573bd2ae [project @ 1999-10-19 10:06:07 by jwe]
jwe
parents:
diff changeset
261 @result{} [ 0.00; 0.004; 0.01 0.00; 0.12; 0.62; 1.00 ]
9067
8970b4b10e9f Cleanup documentation for quad.texi and diffeq.texi
Rik <rdrider0-list@yahoo.com>
parents: 8920
diff changeset
262 @end group
3294
bfe1573bd2ae [project @ 1999-10-19 10:06:07 by jwe]
jwe
parents:
diff changeset
263 @end example
6741
00116015904d [project @ 2007-06-18 16:07:14 by jwe]
jwe
parents: 6502
diff changeset
264
00116015904d [project @ 2007-06-18 16:07:14 by jwe]
jwe
parents: 6502
diff changeset
265 @node Functions of Multiple Variables
00116015904d [project @ 2007-06-18 16:07:14 by jwe]
jwe
parents: 6502
diff changeset
266 @section Functions of Multiple Variables
00116015904d [project @ 2007-06-18 16:07:14 by jwe]
jwe
parents: 6502
diff changeset
267
7984
bbaa5d7d0143 Some documentation updates
David Bateman <dbateman@free.fr>
parents: 7018
diff changeset
268 Octave does not have built-in functions for computing the integral of
12612
16cca721117b doc: Update all documentation for chapter on Numerical Integration
Rik <octave@nomad.inbox5.com>
parents: 11562
diff changeset
269 functions of multiple variables directly. It is possible, however, to
7984
bbaa5d7d0143 Some documentation updates
David Bateman <dbateman@free.fr>
parents: 7018
diff changeset
270 compute the integral of a function of multiple variables using the
12612
16cca721117b doc: Update all documentation for chapter on Numerical Integration
Rik <octave@nomad.inbox5.com>
parents: 11562
diff changeset
271 existing functions for one-dimensional integrals.
6741
00116015904d [project @ 2007-06-18 16:07:14 by jwe]
jwe
parents: 6502
diff changeset
272
00116015904d [project @ 2007-06-18 16:07:14 by jwe]
jwe
parents: 6502
diff changeset
273 To illustrate how the integration can be performed, we will integrate
00116015904d [project @ 2007-06-18 16:07:14 by jwe]
jwe
parents: 6502
diff changeset
274 the function
00116015904d [project @ 2007-06-18 16:07:14 by jwe]
jwe
parents: 6502
diff changeset
275 @tex
00116015904d [project @ 2007-06-18 16:07:14 by jwe]
jwe
parents: 6502
diff changeset
276 $$
00116015904d [project @ 2007-06-18 16:07:14 by jwe]
jwe
parents: 6502
diff changeset
277 f(x, y) = \sin(\pi x y)\sqrt{x y}
00116015904d [project @ 2007-06-18 16:07:14 by jwe]
jwe
parents: 6502
diff changeset
278 $$
00116015904d [project @ 2007-06-18 16:07:14 by jwe]
jwe
parents: 6502
diff changeset
279 @end tex
00116015904d [project @ 2007-06-18 16:07:14 by jwe]
jwe
parents: 6502
diff changeset
280 @ifnottex
10828
322f43e0e170 Grammarcheck .txi documentation files.
Rik <octave@nomad.inbox5.com>
parents: 9245
diff changeset
281
6741
00116015904d [project @ 2007-06-18 16:07:14 by jwe]
jwe
parents: 6502
diff changeset
282 @example
00116015904d [project @ 2007-06-18 16:07:14 by jwe]
jwe
parents: 6502
diff changeset
283 f(x, y) = sin(pi*x*y)*sqrt(x*y)
00116015904d [project @ 2007-06-18 16:07:14 by jwe]
jwe
parents: 6502
diff changeset
284 @end example
10828
322f43e0e170 Grammarcheck .txi documentation files.
Rik <octave@nomad.inbox5.com>
parents: 9245
diff changeset
285
6741
00116015904d [project @ 2007-06-18 16:07:14 by jwe]
jwe
parents: 6502
diff changeset
286 @end ifnottex
00116015904d [project @ 2007-06-18 16:07:14 by jwe]
jwe
parents: 6502
diff changeset
287 for @math{x} and @math{y} between 0 and 1.
00116015904d [project @ 2007-06-18 16:07:14 by jwe]
jwe
parents: 6502
diff changeset
288
00116015904d [project @ 2007-06-18 16:07:14 by jwe]
jwe
parents: 6502
diff changeset
289 The first approach creates a function that integrates @math{f} with
00116015904d [project @ 2007-06-18 16:07:14 by jwe]
jwe
parents: 6502
diff changeset
290 respect to @math{x}, and then integrates that function with respect to
12612
16cca721117b doc: Update all documentation for chapter on Numerical Integration
Rik <octave@nomad.inbox5.com>
parents: 11562
diff changeset
291 @math{y}. Because @code{quad} is written in Fortran it cannot be called
6741
00116015904d [project @ 2007-06-18 16:07:14 by jwe]
jwe
parents: 6502
diff changeset
292 recursively. This means that @code{quad} cannot integrate a function
00116015904d [project @ 2007-06-18 16:07:14 by jwe]
jwe
parents: 6502
diff changeset
293 that calls @code{quad}, and hence cannot be used to perform the double
12612
16cca721117b doc: Update all documentation for chapter on Numerical Integration
Rik <octave@nomad.inbox5.com>
parents: 11562
diff changeset
294 integration. Any of the other integrators, however, can be used which is
16cca721117b doc: Update all documentation for chapter on Numerical Integration
Rik <octave@nomad.inbox5.com>
parents: 11562
diff changeset
295 what the following code demonstrates.
6741
00116015904d [project @ 2007-06-18 16:07:14 by jwe]
jwe
parents: 6502
diff changeset
296
00116015904d [project @ 2007-06-18 16:07:14 by jwe]
jwe
parents: 6502
diff changeset
297 @example
9067
8970b4b10e9f Cleanup documentation for quad.texi and diffeq.texi
Rik <rdrider0-list@yahoo.com>
parents: 8920
diff changeset
298 @group
12612
16cca721117b doc: Update all documentation for chapter on Numerical Integration
Rik <octave@nomad.inbox5.com>
parents: 11562
diff changeset
299 function q = g(y)
16cca721117b doc: Update all documentation for chapter on Numerical Integration
Rik <octave@nomad.inbox5.com>
parents: 11562
diff changeset
300 q = ones (size (y));
16cca721117b doc: Update all documentation for chapter on Numerical Integration
Rik <octave@nomad.inbox5.com>
parents: 11562
diff changeset
301 for i = 1:length (y)
16cca721117b doc: Update all documentation for chapter on Numerical Integration
Rik <octave@nomad.inbox5.com>
parents: 11562
diff changeset
302 f = @@(x) sin (pi*x.*y(i)) .* sqrt (x.*y(i));
16cca721117b doc: Update all documentation for chapter on Numerical Integration
Rik <octave@nomad.inbox5.com>
parents: 11562
diff changeset
303 q(i) = quadgk (f, 0, 1);
6741
00116015904d [project @ 2007-06-18 16:07:14 by jwe]
jwe
parents: 6502
diff changeset
304 endfor
00116015904d [project @ 2007-06-18 16:07:14 by jwe]
jwe
parents: 6502
diff changeset
305 endfunction
00116015904d [project @ 2007-06-18 16:07:14 by jwe]
jwe
parents: 6502
diff changeset
306
12612
16cca721117b doc: Update all documentation for chapter on Numerical Integration
Rik <octave@nomad.inbox5.com>
parents: 11562
diff changeset
307 I = quadgk ("g", 0, 1)
6741
00116015904d [project @ 2007-06-18 16:07:14 by jwe]
jwe
parents: 6502
diff changeset
308 @result{} 0.30022
9067
8970b4b10e9f Cleanup documentation for quad.texi and diffeq.texi
Rik <rdrider0-list@yahoo.com>
parents: 8920
diff changeset
309 @end group
6741
00116015904d [project @ 2007-06-18 16:07:14 by jwe]
jwe
parents: 6502
diff changeset
310 @end example
00116015904d [project @ 2007-06-18 16:07:14 by jwe]
jwe
parents: 6502
diff changeset
311
7984
bbaa5d7d0143 Some documentation updates
David Bateman <dbateman@free.fr>
parents: 7018
diff changeset
312 The above process can be simplified with the @code{dblquad} and
bbaa5d7d0143 Some documentation updates
David Bateman <dbateman@free.fr>
parents: 7018
diff changeset
313 @code{triplequad} functions for integrals over two and three
10828
322f43e0e170 Grammarcheck .txi documentation files.
Rik <octave@nomad.inbox5.com>
parents: 9245
diff changeset
314 variables. For example:
7984
bbaa5d7d0143 Some documentation updates
David Bateman <dbateman@free.fr>
parents: 7018
diff changeset
315
bbaa5d7d0143 Some documentation updates
David Bateman <dbateman@free.fr>
parents: 7018
diff changeset
316 @example
9067
8970b4b10e9f Cleanup documentation for quad.texi and diffeq.texi
Rik <rdrider0-list@yahoo.com>
parents: 8920
diff changeset
317 @group
12612
16cca721117b doc: Update all documentation for chapter on Numerical Integration
Rik <octave@nomad.inbox5.com>
parents: 11562
diff changeset
318 I = dblquad (@@(x, y) sin (pi*x.*y) .* sqrt (x.*y), 0, 1, 0, 1)
7984
bbaa5d7d0143 Some documentation updates
David Bateman <dbateman@free.fr>
parents: 7018
diff changeset
319 @result{} 0.30022
9067
8970b4b10e9f Cleanup documentation for quad.texi and diffeq.texi
Rik <rdrider0-list@yahoo.com>
parents: 8920
diff changeset
320 @end group
7984
bbaa5d7d0143 Some documentation updates
David Bateman <dbateman@free.fr>
parents: 7018
diff changeset
321 @end example
bbaa5d7d0143 Some documentation updates
David Bateman <dbateman@free.fr>
parents: 7018
diff changeset
322
bbaa5d7d0143 Some documentation updates
David Bateman <dbateman@free.fr>
parents: 7018
diff changeset
323 @DOCSTRING(dblquad)
bbaa5d7d0143 Some documentation updates
David Bateman <dbateman@free.fr>
parents: 7018
diff changeset
324
bbaa5d7d0143 Some documentation updates
David Bateman <dbateman@free.fr>
parents: 7018
diff changeset
325 @DOCSTRING(triplequad)
bbaa5d7d0143 Some documentation updates
David Bateman <dbateman@free.fr>
parents: 7018
diff changeset
326
12612
16cca721117b doc: Update all documentation for chapter on Numerical Integration
Rik <octave@nomad.inbox5.com>
parents: 11562
diff changeset
327 The above mentioned approach works, but is fairly slow, and that problem
16cca721117b doc: Update all documentation for chapter on Numerical Integration
Rik <octave@nomad.inbox5.com>
parents: 11562
diff changeset
328 increases exponentially with the dimensionality of the integral. Another
6741
00116015904d [project @ 2007-06-18 16:07:14 by jwe]
jwe
parents: 6502
diff changeset
329 possible solution is to use Orthogonal Collocation as described in the
12612
16cca721117b doc: Update all documentation for chapter on Numerical Integration
Rik <octave@nomad.inbox5.com>
parents: 11562
diff changeset
330 previous section (@pxref{Orthogonal Collocation}). The integral of a function
16cca721117b doc: Update all documentation for chapter on Numerical Integration
Rik <octave@nomad.inbox5.com>
parents: 11562
diff changeset
331 @math{f(x,y)} for @math{x} and @math{y} between 0 and 1 can be approximated
16cca721117b doc: Update all documentation for chapter on Numerical Integration
Rik <octave@nomad.inbox5.com>
parents: 11562
diff changeset
332 using @math{n} points by
6741
00116015904d [project @ 2007-06-18 16:07:14 by jwe]
jwe
parents: 6502
diff changeset
333 @tex
00116015904d [project @ 2007-06-18 16:07:14 by jwe]
jwe
parents: 6502
diff changeset
334 $$
00116015904d [project @ 2007-06-18 16:07:14 by jwe]
jwe
parents: 6502
diff changeset
335 \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),
00116015904d [project @ 2007-06-18 16:07:14 by jwe]
jwe
parents: 6502
diff changeset
336 $$
00116015904d [project @ 2007-06-18 16:07:14 by jwe]
jwe
parents: 6502
diff changeset
337 @end tex
00116015904d [project @ 2007-06-18 16:07:14 by jwe]
jwe
parents: 6502
diff changeset
338 @ifnottex
00116015904d [project @ 2007-06-18 16:07:14 by jwe]
jwe
parents: 6502
diff changeset
339 the sum over @code{i=1:n} and @code{j=1:n} of @code{q(i)*q(j)*f(r(i),r(j))},
00116015904d [project @ 2007-06-18 16:07:14 by jwe]
jwe
parents: 6502
diff changeset
340 @end ifnottex
00116015904d [project @ 2007-06-18 16:07:14 by jwe]
jwe
parents: 6502
diff changeset
341 where @math{q} and @math{r} is as returned by @code{colloc(n)}. The
9067
8970b4b10e9f Cleanup documentation for quad.texi and diffeq.texi
Rik <rdrider0-list@yahoo.com>
parents: 8920
diff changeset
342 generalization to more than two variables is straight forward. The
12612
16cca721117b doc: Update all documentation for chapter on Numerical Integration
Rik <octave@nomad.inbox5.com>
parents: 11562
diff changeset
343 following code computes the studied integral using @math{n=8} points.
6741
00116015904d [project @ 2007-06-18 16:07:14 by jwe]
jwe
parents: 6502
diff changeset
344
00116015904d [project @ 2007-06-18 16:07:14 by jwe]
jwe
parents: 6502
diff changeset
345 @example
9067
8970b4b10e9f Cleanup documentation for quad.texi and diffeq.texi
Rik <rdrider0-list@yahoo.com>
parents: 8920
diff changeset
346 @group
12612
16cca721117b doc: Update all documentation for chapter on Numerical Integration
Rik <octave@nomad.inbox5.com>
parents: 11562
diff changeset
347 f = @@(x,y) sin (pi*x*y') .* sqrt (x*y');
16cca721117b doc: Update all documentation for chapter on Numerical Integration
Rik <octave@nomad.inbox5.com>
parents: 11562
diff changeset
348 n = 8;
16cca721117b doc: Update all documentation for chapter on Numerical Integration
Rik <octave@nomad.inbox5.com>
parents: 11562
diff changeset
349 [t, ~, ~, q] = colloc (n);
6741
00116015904d [project @ 2007-06-18 16:07:14 by jwe]
jwe
parents: 6502
diff changeset
350 I = q'*f(t,t)*q;
00116015904d [project @ 2007-06-18 16:07:14 by jwe]
jwe
parents: 6502
diff changeset
351 @result{} 0.30022
9067
8970b4b10e9f Cleanup documentation for quad.texi and diffeq.texi
Rik <rdrider0-list@yahoo.com>
parents: 8920
diff changeset
352 @end group
6741
00116015904d [project @ 2007-06-18 16:07:14 by jwe]
jwe
parents: 6502
diff changeset
353 @end example
00116015904d [project @ 2007-06-18 16:07:14 by jwe]
jwe
parents: 6502
diff changeset
354
00116015904d [project @ 2007-06-18 16:07:14 by jwe]
jwe
parents: 6502
diff changeset
355 @noindent
00116015904d [project @ 2007-06-18 16:07:14 by jwe]
jwe
parents: 6502
diff changeset
356 It should be noted that the number of points determines the quality
00116015904d [project @ 2007-06-18 16:07:14 by jwe]
jwe
parents: 6502
diff changeset
357 of the approximation. If the integration needs to be performed between
12612
16cca721117b doc: Update all documentation for chapter on Numerical Integration
Rik <octave@nomad.inbox5.com>
parents: 11562
diff changeset
358 @math{a} and @math{b}, instead of 0 and 1, then a change of variables is needed.
6741
00116015904d [project @ 2007-06-18 16:07:14 by jwe]
jwe
parents: 6502
diff changeset
359