Added composite and nested trapezoidal rule

quadrature/trapezoidal_nested.m

@@ -0,0 +1,29 @@
+function [x,w] = trapezoidal_nested(n)
+% TRAPEZOIDAL_NESTED Computes the nested trapezoidal rule.
+%   [X,W] = TRAPEZOIDAL_NESTED(N) computes the trapezoidal rule on 2^(M-1)
+%   intervals. If [X1,W1] = TRAPEZOIDAL_NESTED(N+1) then X is always a
+%   proper subset of X1. This rule is mostly useful for sparse grid
+%   integration if the integrand is not very smooth.
+%
+% Example (<a href="matlab:run_example trapezoidal_nested">run</a>)
+%   clf; hold all
+%   for i = 1:5
+%     [x, w] = trapezoidal_nested(i);
+%     plot(x, i*ones(size(x)), 'x');
+%   end
+%   xlim([-1.3, 1.3]); ylim([0.5, 5.5])
+%
+% See also TRAPEZOIDAL_RULE, SMOLYAK_GRID, INTEGRATE_ND
+
+%   Elmar Zander
+%   Copyright 2013, Inst. of Scientific Computing, TU Braunschweig
+%
+%   This program is free software: you can redistribute it and/or modify it
+%   under the terms of the GNU General Public License as published by the
+%   Free Software Foundation, either version 3 of the License, or (at your
+%   option) any later version. 
+%   See the GNU General Public License for more details. You should have
+%   received a copy of the GNU General Public License along with this
+%   program.  If not, see <http://www.gnu.org/licenses/>.
+
+[x, w] = trapezoidal_rule(2^(n-1));

quadrature/trapezoidal_rule.m

@@ -0,0 +1,31 @@
+function [x,w]=trapezoidal_rule(n)
+% TRAPEZOIDAL_RULE Compute points and weights of the trapezoidal rule.
+%   [X,W]=TRAPEZOIDAL_RULE(N) computes the points and weights of the
+%   composite trapezoidal rule with N+1 integration points (i.e. on N
+%   subintervals), returning the points in the 1 x N+1 array X and the
+%   weights in the N+1 x 1 array w.
+%
+% Example (<a href="matlab:run_example trapezoidal_rule">run</a>)
+%   [x,w] = trapezoidal_rule(20);
+%   fprintf('Integral of cosine from -1 to 1:\n');
+%   fprintf('I_ex      = %g\n', sin(1) - sin(-1))
+%   fprintf('I_trap_20 = %g\n', cos(x) * w)
+%   fprintf('I_matlab  = %g\n', quad(@cos, -1, 1))
+%
+% See also INTEGRATE_1D, GAUSS_HERMITE_LEGENDRE_RULE, QUAD
+
+%   Elmar Zander
+%   Copyright 2013, Inst. of Scientific Computing, TU Braunschweig
+%
+%   This program is free software: you can redistribute it and/or modify it
+%   under the terms of the GNU General Public License as published by the
+%   Free Software Foundation, either version 3 of the License, or (at your
+%   option) any later version. 
+%   See the GNU General Public License for more details. You should have
+%   received a copy of the GNU General Public License along with this
+%   program.  If not, see <http://www.gnu.org/licenses/>.
+
+
+x = linspace(-1, 1, n+1);
+w0 = 1.0 / n * ones(n, 1);
+w = [0;w0] + [w0;0];

quadrature/unittest_trapezoidal_rule.m

@@ -0,0 +1,37 @@
+function unittest_trapezoidal_rule
+% UNITTEST_TRAPEZOIDAL_RULE Test the TRAPEZOIDAL_RULE function.
+%
+% Example (<a href="matlab:run_example unittest_trapezoidal_rule">run</a>)
+%   unittest_trapezoidal_rule
+%
+% See also TRAPEZOIDAL_RULE, TESTSUITE 
+
+%   Elmar Zander
+%   Copyright 2013, Inst. of Scientific Computing, TU Braunschweig
+%
+%   This program is free software: you can redistribute it and/or modify it
+%   under the terms of the GNU General Public License as published by the
+%   Free Software Foundation, either version 3 of the License, or (at your
+%   option) any later version. 
+%   See the GNU General Public License for more details. You should have
+%   received a copy of the GNU General Public License along with this
+%   program.  If not, see <http://www.gnu.org/licenses/>.
+
+munit_set_function( 'trapezoidal_rule' );
+
+n = 15;
+[x,w] = trapezoidal_rule(n);
+assert_equals(size(x), [1, n+1], 'size_x');
+assert_equals(size(w), [n+1, 1], 'size_w');
+assert_equals(sum(w), 2, 'sum_w');
+assert_equals(diff(diff(x)), zeros(1,n-1), 'equidist');
+assert_equals(x([1,n+1]), [-1, 1], 'endpoints');
+
+munit_set_function( 'trapezoidal_nested' );
+
+assert_equals(trapezoidal_nested(1), [-1, 1], 'nested_1');
+for i = 1:5
+    x1 = trapezoidal_nested(i);
+    x2 = trapezoidal_nested(i+1);
+    assert_true(all(ismember(x1, x2)), [], sprintf('is_nested_%d', i));
+end