1 | function anss = gaussquad(f,a,b,tol) |
---|
2 | % GAUSSQUAD(f,a,b,tol) Definite integration using Gauss-Legendre |
---|
3 | % quadrature. Finds the definite integral of function f from a to b. |
---|
4 | % User can define the relative convergence tolerence, tol; the default is |
---|
5 | % 10^-14. If f is an inline function, it needs to accept vector args. |
---|
6 | % Function f can also be ananymous or a function handle. |
---|
7 | % |
---|
8 | % Example: |
---|
9 | % |
---|
10 | % f = @(x) sin(x)./x+exp(-x); % Vectorized function handle. |
---|
11 | % gaussquad(f,-1,1) |
---|
12 | % |
---|
13 | % See also quad, quadl, quadv and rombquad (on the FEX) |
---|
14 | % |
---|
15 | % Author: Matt Fig |
---|
16 | % Date: revised 1/31/2006 |
---|
17 | % Contact: popkenai@yahoo.com |
---|
18 | |
---|
19 | |
---|
20 | % Check for infinite boundaries. |
---|
21 | if nargin<3 |
---|
22 | error ('Not enough arguments. See help.') |
---|
23 | elseif nargin < 4 |
---|
24 | tol = 1e-14; % Default. |
---|
25 | end |
---|
26 | |
---|
27 | if isinf(f(a)) || isinf(f(b)) |
---|
28 | error('Infinite limits not allowed.') |
---|
29 | end |
---|
30 | |
---|
31 | |
---|
32 | % Initially use a 40 point Gauss-Legendre quadrature in 15 subintervals. |
---|
33 | gp = 40; |
---|
34 | cnt = 1; |
---|
35 | ints = 15; |
---|
36 | maxcount = 35; |
---|
37 | anss1 = 1; |
---|
38 | anss2 = 0; % Initialization for while. |
---|
39 | |
---|
40 | % Main loop. Each pass through, the number of intervals and Gauss points |
---|
41 | % is increased until either the tol is met or the maxcount is exceeded. |
---|
42 | while abs(anss2-anss1) >= tol && cnt < maxcount |
---|
43 | anss1 = core(f,a,b,gp,ints); |
---|
44 | anss2 = core(f,a,b,gp+8,ints+3); |
---|
45 | gp = gp+16; |
---|
46 | ints = ints+6; |
---|
47 | cnt = cnt+1; |
---|
48 | end |
---|
49 | |
---|
50 | % Give a warning if max iterations is reached and give an error estimate. |
---|
51 | if cnt>=maxcount |
---|
52 | str = sprintf(['Maximum iterations reached, results may be\n'... |
---|
53 | ' inaccurate. Estimated error: %.2d'],... |
---|
54 | abs(anss2-anss1)); |
---|
55 | warning(str); %#ok |
---|
56 | end |
---|
57 | |
---|
58 | anss = anss2; |
---|
59 | |
---|
60 | |
---|
61 | function anss = core(f,a,b,gp,ints) |
---|
62 | % Core subfunction that performs the integrations. |
---|
63 | [abs1, wgt1] = Gauss(gp); % Get Gauss points. |
---|
64 | bb(1) = a; % Get subintervals. |
---|
65 | bb(2:ints) = (2:ints)*(b-a)/ints+a; |
---|
66 | dif = diff(bb)/2; |
---|
67 | % Evaluate f at scaled intervals, need to map each interval to [-1 1]. |
---|
68 | an = f((abs1+1)*dif+repmat(bb(1:end-1),gp,1)); % Function evaluations. |
---|
69 | new = dif*an'*wgt1; % Multiply by the weights and intervals. |
---|
70 | anss = sum(new(:)); |
---|
71 | |
---|
72 | |
---|
73 | function [x, w] = Gauss(n) |
---|
74 | % Generates the abscissa and weights for a Gauss-Legendre quadrature. |
---|
75 | % Reference: Numerical Recipes in Fortran 77, Cornell press. |
---|
76 | x = zeros(n,1); % Preallocations. |
---|
77 | w = x; |
---|
78 | m = (n+1)/2; |
---|
79 | for ii=1:m |
---|
80 | z = cos(pi*(ii-.25)/(n+.5)); % Initial estimate. |
---|
81 | z1 = z+1; |
---|
82 | while abs(z-z1)>eps |
---|
83 | p1 = 1; |
---|
84 | p2 = 0; |
---|
85 | for jj = 1:n |
---|
86 | p3 = p2; |
---|
87 | p2 = p1; |
---|
88 | p1 = ((2*jj-1)*z*p2-(jj-1)*p3)/jj; % The Legendre polynomial. |
---|
89 | end |
---|
90 | pp = n*(z*p1-p2)/(z^2-1); % The L.P. derivative. |
---|
91 | z1 = z; |
---|
92 | z = z1-p1/pp; |
---|
93 | end |
---|
94 | x(ii) = -z; % Build up the abscissas. |
---|
95 | x(n+1-ii) = z; |
---|
96 | w(ii) = 2/((1-z^2)*(pp^2)); % Build up the weights. |
---|
97 | w(n+1-ii) = w(ii); |
---|
98 | end |
---|