|
2670
|
1 @c Copyright (C) 1996, 1997 John W. Eaton |
|
2333
|
2 @c This is part of the Octave manual. |
|
|
3 @c For copying conditions, see the file gpl.texi. |
|
|
4 |
|
2653
|
5 @node Quadrature, Differential Equations, Nonlinear Equations, Top |
|
2333
|
6 @chapter Quadrature |
|
|
7 |
|
|
8 @menu |
|
2670
|
9 * Functions of One Variable:: |
|
2333
|
10 * Orthogonal Collocation:: |
|
|
11 @end menu |
|
|
12 |
|
2670
|
13 @node Functions of One Variable, Orthogonal Collocation, Quadrature, Quadrature |
|
|
14 @section Functions of One Variable |
|
2333
|
15 |
|
2670
|
16 @deftypefn {Loadable Function} {[@var{v}, @var{ier}, @var{nfun}, @var{err}] =} quad (@var{f}, @var{a}, @var{b}, @var{tol}, @var{sing}) |
|
2333
|
17 Integrate a nonlinear function of one variable using Quadpack. |
|
2449
|
18 The first argument is the name of the function to call to compute the |
|
|
19 value of the integrand. It must have the form |
|
2333
|
20 |
|
|
21 @example |
|
|
22 y = f (x) |
|
|
23 @end example |
|
|
24 |
|
|
25 @noindent |
|
2449
|
26 where @var{y} and @var{x} are scalars. |
|
2333
|
27 |
|
|
28 The second and third arguments are limits of integration. Either or |
|
|
29 both may be infinite. |
|
|
30 |
|
2449
|
31 The optional argument @var{tol} is a vector that specifies the desired |
|
2333
|
32 accuracy of the result. The first element of the vector is the desired |
|
|
33 absolute tolerance, and the second element is the desired relative |
|
|
34 tolerance. To choose a relative test only, set the absolute |
|
|
35 tolerance to zero. To choose an absolute test only, set the relative |
|
|
36 tolerance to zero. |
|
|
37 |
|
|
38 The optional argument @var{sing} is a vector of values at which the |
|
2449
|
39 integrand is known to be singular. |
|
2670
|
40 |
|
|
41 The result of the integration is returned in @var{v} and @var{ier} |
|
|
42 contains an integer error code (0 indicates a successful integration). |
|
|
43 The value of @var{nfun} indicates how many function evaluations were |
|
|
44 required, and @var{err} contains an estimate of the error in the |
|
|
45 solution. |
|
2449
|
46 @end deftypefn |
|
2333
|
47 |
|
2465
|
48 @deftypefn {Loadable Function} {} quad_options (@var{opt}, @var{val}) |
|
2449
|
49 When called with two arguments, this function allows you set options |
|
|
50 parameters for the function @code{quad}. Given one argument, |
|
|
51 @code{quad_options} returns the value of the corresponding option. If |
|
|
52 no arguments are supplied, the names of all the available options and |
|
|
53 their current values are displayed. |
|
|
54 @end deftypefn |
|
2333
|
55 |
|
2670
|
56 Here is an example of using @code{quad} to integrate the function |
|
|
57 @iftex |
|
|
58 @tex |
|
|
59 $$ |
|
|
60 f(x) = x \sin (1/x) \sqrt {|1 - x|} |
|
|
61 $$ |
|
|
62 from $x = 0$ to $x = 3$. |
|
|
63 @end tex |
|
|
64 @end iftex |
|
|
65 @ifinfo |
|
2689
|
66 |
|
2670
|
67 @example |
|
|
68 @var{f}(@var{x}) = @var{x} * sin (1/@var{x}) * sqrt (abs (1 - @var{x})) |
|
|
69 @end example |
|
|
70 |
|
|
71 @noindent |
|
|
72 from @var{x} = 0 to @var{x} = 3. |
|
|
73 @end ifinfo |
|
|
74 |
|
|
75 This is a fairly difficult integration (plot the function over the range |
|
|
76 of integration to see why). |
|
|
77 |
|
|
78 The first step is to define the function: |
|
|
79 |
|
|
80 @example |
|
|
81 @group |
|
|
82 function y = f (x) |
|
|
83 y = x .* sin (1 ./ x) .* sqrt (abs (1 - x)); |
|
|
84 endfunction |
|
|
85 @end group |
|
|
86 @end example |
|
|
87 |
|
|
88 Note the use of the `dot' forms of the operators. This is not necessary |
|
|
89 for the call to @code{quad}, but it makes it much easier to generate a |
|
2689
|
90 set of points for plotting (because it makes it possible to call the |
|
|
91 function with a vector argument to produce a vector result). |
|
2670
|
92 |
|
|
93 Then we simply call quad: |
|
|
94 |
|
|
95 @example |
|
|
96 @group |
|
|
97 [v, ier, nfun, err] = quad ("f", 0, 3) |
|
|
98 @result{} 1.9819 |
|
|
99 @result{} 1 |
|
|
100 @result{} 5061 |
|
|
101 @result{} 1.1522e-07 |
|
|
102 @end group |
|
|
103 @end example |
|
|
104 |
|
|
105 Although @code{quad} returns a nonzero value for @var{ier}, the result |
|
2689
|
106 is reasonably accurate (to see why, examine what happens to the result |
|
|
107 if you move the lower bound to 0.1, then 0.01, then 0.001, etc.). |
|
2670
|
108 |
|
|
109 @node Orthogonal Collocation, , Functions of One Variable, Quadrature |
|
2333
|
110 @section Orthogonal Collocation |
|
|
111 |
|
2465
|
112 @deftypefn {Loadable Function} {[@var{r}, @var{A}, @var{B}, @var{q}] =} colloc (@var{n}, "left", "right") |
|
2333
|
113 Compute derivative and integral weight matrices for orthogonal |
|
|
114 collocation using the subroutines given in J. Villadsen and |
|
|
115 M. L. Michelsen, @cite{Solution of Differential Equation Models by |
|
|
116 Polynomial Approximation}. |
|
2449
|
117 @end deftypefn |
|
2689
|
118 |
|
|
119 Here is an example of using @code{colloc} to generate weight matrices |
|
|
120 for solving the second order differential equation |
|
|
121 @iftex |
|
|
122 @tex |
|
|
123 $u^\prime - \alpha u^{\prime\prime} = 0$ with the boundary conditions |
|
|
124 $u(0) = 0$ and $u(1) = 1$. |
|
|
125 @end tex |
|
|
126 @end iftex |
|
|
127 @ifinfo |
|
|
128 @var{u}' - @var{alpha} * @var{u}'' = 0 with the boundary conditions |
|
|
129 @var{u}(0) = 0 and @var{u}(1) = 1. |
|
|
130 @end ifinfo |
|
|
131 |
|
|
132 First, we can generate the weight matrices for @var{n} points (including |
|
|
133 the endpoints of the interval), and incorporate the boundary conditions |
|
|
134 in the right hand side (for a specific value of |
|
|
135 @iftex |
|
|
136 @tex |
|
|
137 $\alpha$). |
|
|
138 @end tex |
|
|
139 @end iftex |
|
|
140 @ifinfo |
|
|
141 @var{alpha}). |
|
|
142 @end ifinfo |
|
|
143 |
|
|
144 @example |
|
|
145 @group |
|
|
146 n = 7; |
|
|
147 alpha = 0.1; |
|
|
148 [r, a, b] = colloc (n-2, "left", "right"); |
|
|
149 at = a(2:n-1,2:n-1); |
|
|
150 bt = b(2:n-1,2:n-1); |
|
|
151 rhs = alpha * b(2:n-1,n) - a(2:n-1,n); |
|
|
152 @end group |
|
|
153 @end example |
|
|
154 |
|
|
155 Then the solution at the roots @var{r} is |
|
|
156 |
|
|
157 @example |
|
|
158 u = [ 0; (at - alpha * bt) \ rhs; 1] |
|
|
159 @result{} [ 0.00; 0.004; 0.01 0.00; 0.12; 0.62; 1.00 ] |
|
|
160 @end example |
