Work on Clenshaw-Curtis (now nested and normal)

demo/publish/show_smolyak.m

@@ -13,7 +13,7 @@
 
 %% Smolyak grid for Clenshaw-Curtis rules (2D)
 
-rule=@clenshaw_curtis_legendre_rule;
+rule=@clenshaw_curtis_nested;
 dim=2;
 
 % Smolyak grid with 4 stages
@@ -35,7 +35,7 @@
 return
 %% Smolyak grid for Clenshaw-Curtis rules (3D)
 
-rule=@clenshaw_curtis_legendre_rule;
+rule=@clenshaw_curtis_nested;
 dim=3;
 
 % Smolyak grid with 4 stages
@@ -62,7 +62,7 @@
       @gauss_legendre_rule, 'legendre';
       ; 'clenshaw-curtis' }
 
-for rule={@gauss_hermite_rule, @gauss_legendre_rule, @clenshaw_curtis_legendre_rule }
+for rule={@gauss_hermite_rule, @gauss_legendre_rule, @clenshaw_curtis_nested}
       for i=1:4
         subplot(2,2,i);
         plot(xd(1,:),xd(2,:),'*k')

demo/publish/show_smolyak_grids.m

@@ -20,7 +20,7 @@
 
 %% Smolyak grid for Clenshaw-Curtis rules (2D)
 
-rule=@clenshaw_curtis_legendre_rule;
+rule=@clenshaw_curtis_nested;
 dim=2;
 
 % Smolyak grid with 4 stages

quadrature/Contents.m

@@ -1,20 +1,42 @@
 % QUADRATURE
 %
-% Files
-%   clenshaw_curtis_legendre_rule - quad1d_cc_legendre - nodes and weights of clenshaw-curtis formula
-%   full_tensor_grid              - Return nodes and weights for full tensor product grid.
-%   gauss_hermite                 - Numerically integrate with Gauss-Hermite quadrature rule.
-%   gauss_hermite_rule            - Return the Gauss-Hermite quadrature rule with p nodes.
-%   gauss_legendre_rule           - Get Gauss points and weights for quadrature over [-1,1].
-%   integrate                     - Integrate a multivariate function.
-%   integrate_1d                  - Integrate a univariate function.
-%   smolyak_grid                  - Return nodes weights for Smolyak quadrature.
-%   tensor_mesh                   - Create D-dimensional tensor-product from 1D meshes and weights.
-%   unittest_full_tensor_grid     - Test the FULL_TENSOR_GRID function.
-%   unittest_gauss_hermite        - Test the GAUSS_HERMITE function.
-%   unittest_gauss_hermite_rule   - Test the GAUSS_HERMITE_RULE function.
-%   unittest_gauss_legendre       - Test the Gauss-Legendere quadrature methods.
-%   unittest_integrate            - Test the INTEGRATE function.
-%   unittest_integrate_1d         - Test the INTEGRATE_1D function.
-%   unittest_smolyak_grid         - Test the SMOLYAK_GRID function.
-%   unittest_tensor_mesh          - Test the TENSOR_MESH function.
+% Basic 1D quadrature rules all end in "_rule" and return the nodes as a
+% row vector and the weights as a column vector (which is often the other
+% way round, but works better with the rest of sglib). The functions ending
+% in "_nested" call the basic rule functions such that the nodes of rule
+% number N are a subset of those of rule number N+1 (i.e. they are nested
+% or embedded). The 1D rules can be used as basic quadrature formulas in
+% functions to create high-dimensional integration rules based e.g. on full
+% tensor or on Smolyak grids.
+%
+% Basic 1D Rules
+%   clenshaw_curtis_nested                - Compute the nested Clenshaw-Curtis rule.
+%   clenshaw_curtis_rule                  - Compute nodes and weights of the Clenshaw-Curtis rules.
+%   gauss_hermite_rule                    - Return the Gauss-Hermite quadrature rule with p nodes.
+%   gauss_legendre_rule                   - Get Gauss points and weights for quadrature over [-1,1].
+%   trapezoidal_nested                    - Computes the nested trapezoidal rule.
+%   trapezoidal_rule                      - Compute points and weights of the trapezoidal rule.
+%   gauss_legendre_triangle_rule          - Get Gauss points and weights for quadrature over canonical triangle.
+%
+% Grid generation for high-dimensional quadrature
+%   full_tensor_grid                      - Return nodes and weights for full tensor product grid.
+%   smolyak_grid                          - Return nodes weights for Smolyak quadrature.
+%   tensor_mesh                           - Create D-dimensional tensor-product from 1D meshes and weights.
+%
+% Easy-to-use integration methods
+%   gauss_hermite                         - Numerically integrate with Gauss-Hermite quadrature rule.
+%   integrate_1d                          - Integrate a univariate function.
+%   integrate_nd                          - Integrate a multivariate function.
+%
+% Unittests
+%   unittest_clenshaw_curtis_rule         - Test the CLENSHAW_CURTIS_RULE function.
+%   unittest_full_tensor_grid             - Test the FULL_TENSOR_GRID function.
+%   unittest_gauss_hermite                - Test the GAUSS_HERMITE function.
+%   unittest_gauss_hermite_rule           - Test the GAUSS_HERMITE_RULE function.
+%   unittest_gauss_legendre_rule          - UNITTEST_GAUSS_LEGENDRE Test the Gauss-Legendere quadrature methods.
+%   unittest_gauss_legendre_triangle_rule - Test the GAUSS_LEGENDRE_TRIANGLE_RULE function.
+%   unittest_integrate_1d                 - Test the INTEGRATE_1D function.
+%   unittest_integrate_nd                 - Test the INTEGRATE_ND function.
+%   unittest_smolyak_grid                 - Test the SMOLYAK_GRID function.
+%   unittest_tensor_mesh                  - Test the TENSOR_MESH function.
+%   unittest_trapezoidal_rule             - Test the TRAPEZOIDAL_RULE function.

quadrature/clenshaw_curtis_nested.m

@@ -0,0 +1,29 @@
+function [x,w]=clenshaw_curtis_nested(n)
+% CLENSHAW_CURTIS_NESTED Compute the nested Clenshaw-Curtis rule.
+%   [X,W]=CLENSHAW_CURTIS_NESTED(N) computes the Clenshaw-Curtis rule of
+%   order M=2^(N-1)+1. This is mainly for sparse grid integration where
+%   nested rules are advantageous.
+%
+% Example (<a href="matlab:run_example clenshaw_curtis_nested">run</a>)
+%   clf; hold all
+%   for i = 1:5
+%     [x, w] = clenshaw_curtis_nested(i);
+%     plot(x, i*ones(size(x)), 'x');
+%   end
+%   xlim([-1.3, 1.3]); ylim([0.5, 5.5])
+%
+% See also SMOLYAK_GRID, INTEGRATE_ND, CLENSHAW_CURTIS_RULE
+
+%   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/>.
+
+m = 2^(n-1)+1;
+[x,w]=clenshaw_curtis_rule(m);

quadrature/clenshaw_curtis_rule.m

@@ -0,0 +1,83 @@
+function [x,w]=clenshaw_curtis_rule(n)
+% CLENSHAW_CURTIS_RULE Compute nodes and weights of the Clenshaw-Curtis rules.
+%   [X,W]=CLENSHAW_CURTIS_RULE(N) 
+%
+% References
+%   [1] J. Waldvogel: Fast Construction of the Fejer and Clenshaw-Curtis
+%       Quadrature Rules, BIT (46) 2006, pp 195-202
+%       doi: 10.1007/s10543-006-0045-4
+%
+% Example (<a href="matlab:run_example clenshaw_curtis_rule">run</a>)
+%
+% See also
+
+%   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/>.
+
+if n == 1
+    x = 0;
+    w = 2;
+elseif n == 2
+    x = [-1, 1];
+    w = [1; 1];
+else
+    x = sin(pi * linspace(-0.5, 0.5, n));
+    w = cc_fejer_weights(n, 3);
+end
+
+function [x,w]=clencurt(n)
+% Computes the Clenshaw Curtis nodes and weights
+% Adapted a code by G. von Winckel
+% http://www.scientificpython.net/1/post/2012/04/clenshaw-curtis-quadrature.html
+% Problem for n==2
+if n == 1
+    x = 0;
+    w = 2;
+else
+    C = zeros(n,2);
+    k = 2*(1:floor((n-1)/2));
+    C(1:2:end,1) = 2./[1, 1-k.*k];
+    C(2,2) = -(n-1);
+    V = [C; flipud(C(2:n-1,:))];
+    F = real(ifft(V)); %, n=None, axis=0))
+    x = F(1:n,2)';
+    w = [F(1,1); 2*F(2:n-1,1); F(n,1)];
+end
+
+
+function w = cc_fejer_weights(n, mode)
+[wf1,wf2,wcc] = fejer(n-1);
+switch mode
+    case 1
+        w = wf1;
+    case 2
+        w = wf2(2:end);
+    case 3
+        w = [wcc; wcc(1)];
+    otherwise
+        error('foo');
+end
+
+
+function [wf1,wf2,wcc] = fejer(n)
+% FEJER Computes Fejer and CC weights (taken from [1])
+% Weights of the Fejer2, Clenshaw-Curtis and Fejer1 quadratures
+% by DFTs. Nodes: x_k = cos(k*pi/n), n>1
+N=(1:2:n-1)'; l=length(N); m=n-l; K=(0:m-1)';
+% Fejer2 nodes: k=0,1,...,n; weights: wf2, wf2_n=wf2_0=0
+v0=[2./N./(N-2); 1/N(end); zeros(m,1)];
+v2=-v0(1:end-1)-v0(end:-1:2); wf2=ifft(v2);
+%Clenshaw-Curtis nodes: k=0,1,...,n; weights: wcc, wcc_n=wcc_0
+g0=-ones(n,1); g0(1+l)=g0(1+l)+n; g0(1+m)=g0(1+m)+n;
+g=g0/(n^2-1+mod(n,2)); wcc=ifft(v2+g);
+% Fejer1 nodes: k=1/2,3/2,...,n-1/2; vector of weights: wf1
+v0=[2*exp(1i*pi*K/n)./(1-4*K.^2); zeros(l+1,1)];
+v1=v0(1:end-1)+conj(v0(end:-1:2)); wf1=ifft(v1);

quadrature/full_tensor_grid.m

@@ -18,7 +18,7 @@
 %   subplot(2,2,2); plot(xd(1,:),xd(2,:),'*k')
 %   [xd,wd]=full_tensor_grid( [], 5, {@gauss_hermite_rule; @gauss_legendre_rule}  );
 %   subplot(2,2,3); plot(xd(1,:),xd(2,:),'*k')
-%   [xd,wd]=full_tensor_grid( [], [6 6 3], {@gauss_hermite_rule; @gauss_legendre_rule; @clenshaw_curtis_legendre_rule}  );
+%   [xd,wd]=full_tensor_grid( [], [6 6 3], {@gauss_hermite_rule; @gauss_legendre_rule; @clenshaw_curtis_nested}  );
 %   subplot(2,2,4); plot3(xd(1,:),xd(2,:),xd(3,:),'*k')
 %
 % See also SMOLYAK_GRID, TENSOR_MESH

quadrature/smolyak_grid.m

@@ -22,7 +22,7 @@
 %      option for compatibility with old A. K. versions.)
 %
 % Example (<a href="matlab:run_example smolyak_grid">run</a>)
-%   for rule={@gauss_hermite_rule, @gauss_legendre_rule, @clenshaw_curtis_legendre_rule }
+%   for rule={@gauss_hermite_rule, @gauss_legendre_rule, @clenshaw_curtis_nested}
 %     for i=1:4
 %       [xd,wd]=smolyak_grid( 2, 3+i, rule );
 %       subplot(2,2,i);

quadrature/unittest_clenshaw_curtis_rule.m

@@ -0,0 +1,60 @@
+function unittest_clenshaw_curtis_rule
+% UNITTEST_CLENSHAW_CURTIS_RULE Test the CLENSHAW_CURTIS_RULE function.
+%
+% Example (<a href="matlab:run_example unittest_clenshaw_curtis_rule">run</a>)
+%   unittest_clenshaw_curtis_rule
+%
+% See also CLENSHAW_CURTIS_RULE, TESTSUITE 
+
+%   Elmar Zander
+%   Copyright 2010, Inst. of Scientific Computing, TU Braunschweig
+%   $Id$ 
+%
+%   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( 'clenshaw_curtis_rule' );
+
+
+[x,w]=clenshaw_curtis_rule(1);
+assert_equals( size(x), [1, 1], 'x1_row' );
+assert_equals( size(w), [1, 1], 'w1_column' );
+assert_equals( sum(w), 2, 'w1_sum' );
+
+[x,w]=clenshaw_curtis_rule(2);
+assert_equals( size(x), [1, 2], 'x2_row2' );
+assert_equals( size(w), [2, 1], 'w2_column' );
+assert_equals( sum(w), 2, 'w2_sum' );
+
+[x,w]=clenshaw_curtis_rule(3);
+assert_equals( x, [-1, 0, 1], 'x3' );
+assert_equals( w, [1; 4; 1]/3, 'w3' );
+
+n = 17;
+[x,w]=clenshaw_curtis_rule(n);
+assert_equals( size(x), [1, 17], 'xn_row' );
+assert_equals( size(w), [17, 1], 'wn_column' );
+assert_equals( sum(w), 2, 'wn_sum' );
+
+assert_equals( sum(x*w), 0, 'int_ord1_ex' );
+assert_equals( sum((x.^2)*w), 2/3, 'int_ord2_ex' );
+assert_equals( sum((x.^3)*w), 0, 'int_ord3_ex' );
+assert_equals( sum((x.^4)*w), 2/5, 'int_ord4_ex' );
+assert_equals( sum((x.^14)*w), 2/15, 'int_ord14_ex' );
+assert_equals( sum((x.^16)*w), 2/17, 'int_ord16_ex' );
+
+
+
+munit_set_function( 'clenshaw_curtis_nested' );
+
+assert_equals(clenshaw_curtis_nested(1), [-1, 1], 'nested_1');
+for i = 1:5
+    x1 = clenshaw_curtis_nested(i);
+    x2 = clenshaw_curtis_nested(i+1);
+    assert_true(all(ismember(x1, x2)), [], sprintf('is_nested_%d', i));
+end

sglib_startup.m

@@ -59,10 +59,10 @@
 
 % show greeting if user wants that
 if appdata.settings.show_greeting
-    [versionstr, msg] = sglib_version('as_string', true);
+    [versionstr, msgs] = sglib_version('as_string', true);
     fprintf( '\nSGLIB v%s\n', versionstr );
-    if ~isempty(msg)
-        fprintf( '%s\n', msg );
+    for msg = msgs
+        fprintf( '%s\n', msg{1} );
     end
     fprintf( '\nChecking toolboxes:\n' );
 end