diff --git a/demo/publish/show_smolyak.m b/demo/publish/show_smolyak.m index 0f3ac584..9bbf505a 100644 --- a/demo/publish/show_smolyak.m +++ b/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') diff --git a/demo/publish/show_smolyak_grids.m b/demo/publish/show_smolyak_grids.m index 8397097b..6801b392 100644 --- a/demo/publish/show_smolyak_grids.m +++ b/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 diff --git a/quadrature/Contents.m b/quadrature/Contents.m index 6d2c34ed..b3cbfb80 100644 --- a/quadrature/Contents.m +++ b/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. diff --git a/quadrature/clenshaw_curtis_legendre_rule.m b/quadrature/clenshaw_curtis_legendre_rule.m deleted file mode 100644 index 8e72b39c..00000000 --- a/quadrature/clenshaw_curtis_legendre_rule.m +++ /dev/null @@ -1,30 +0,0 @@ -function [y,w]=clenshaw_curtis_legendre_rule(k) -% quad1d_cc_legendre - nodes and weights of clenshaw-curtis formula -% for Legendre measure - -%n = floor(order/2)+1; -%assert(k>0); - -if k==1 - y=0; - w=1; - return -end - -% order: -m = 2^(k-1)+1; - -% -J = 1:m; -y = -cos(pi * (J-1)/(m-1)); - -w = 1 - cos(pi * (J-1))/(m*(m-2)); -for k=1:(m-3)/2 - w = w - 2*1/(4*k^2-1)*cos(2 * pi * k * (J-1)/(m-1)); -end - -w = w * 2/(m-1); -w=w'; - -w(1)=1/(m*(m-2)); -w(m)=w(1); diff --git a/quadrature/clenshaw_curtis_nested.m b/quadrature/clenshaw_curtis_nested.m new file mode 100644 index 00000000..967677e2 --- /dev/null +++ b/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 (run) +% 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 . + +m = 2^(n-1)+1; +[x,w]=clenshaw_curtis_rule(m); diff --git a/quadrature/clenshaw_curtis_rule.m b/quadrature/clenshaw_curtis_rule.m new file mode 100644 index 00000000..0bde290f --- /dev/null +++ b/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 (run) +% +% 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 . + +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); diff --git a/quadrature/full_tensor_grid.m b/quadrature/full_tensor_grid.m index 83be9521..7847009e 100644 --- a/quadrature/full_tensor_grid.m +++ b/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 diff --git a/quadrature/smolyak_grid.m b/quadrature/smolyak_grid.m index 146b7e95..6c725217 100644 --- a/quadrature/smolyak_grid.m +++ b/quadrature/smolyak_grid.m @@ -22,7 +22,7 @@ % option for compatibility with old A. K. versions.) % % Example (run) -% 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); diff --git a/quadrature/unittest_clenshaw_curtis_legendre_rule.m b/quadrature/unittest_clenshaw_curtis_legendre_rule.m deleted file mode 100644 index 54c33932..00000000 --- a/quadrature/unittest_clenshaw_curtis_legendre_rule.m +++ /dev/null @@ -1,33 +0,0 @@ -function unittest_clenshaw_curtis_legendre_rule -% UNITTEST_CLENSHAW_CURTIS_LEGENDRE_RULE Test the CLENSHAW_CURTIS_LEGENDRE_RULE function. -% -% Example (run) -% unittest_clenshaw_curtis_legendre_rule -% -% See also CLENSHAW_CURTIS_LEGENDRE_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 . - -munit_set_function( 'clenshaw_curtis_legendre_rule' ); - -k=5; -[x,w]=clenshaw_curtis_legendre_rule(k); -assert_equals( size(x,1), 1, 'x_row' ); -assert_equals( size(w,2), 1, 'w_column' ); -assert_equals( sum(w), max(x)-min(x), 'w_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' ); diff --git a/quadrature/unittest_clenshaw_curtis_rule.m b/quadrature/unittest_clenshaw_curtis_rule.m new file mode 100644 index 00000000..aabfa3a1 --- /dev/null +++ b/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 (run) +% 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 . + +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 diff --git a/sglib_startup.m b/sglib_startup.m index c71d6868..e5cfa206 100644 --- a/sglib_startup.m +++ b/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 diff --git a/sglib_version.m b/sglib_version.m index cf411386..9d2dac31 100644 --- a/sglib_version.m +++ b/sglib_version.m @@ -1,8 +1,16 @@ -function [version, msg]=sglib_version(varargin) +function [version, msgs]=sglib_version(varargin) % SGLIB_VERSION Returns version information for sglib. -% SGLIB_VERSION Returns version information for sglib either as array or -% in string format (if the option 'as_string' is specified). +% VERSION=SGLIB_VERSION Returns version information for sglib either as +% array or in string format (if the option 'as_string' is specified). % +% VERSION=SGLIB_VERSION('as_string', true) Returns version information +% for sglib in string format. +% +% [VERSION, MSGS]=SGLIB_VERSION Returns an addtional cell array +% containing string messages with important information about possibly +% incompatible changes. +% +% Note: % The code of SGLIB_VERSION will also contain a log of the major % improvements of sglib (some kind of 'whats_new' file. % @@ -12,7 +20,7 @@ % fprintf('Version as string: '); % disp(sglib_version('as_string', true)) % -% See also +% See also VERSION, SGLIB_STARTUP % Elmar Zander % Copyright 2013, Inst. of Scientific Computing, TU Braunschweig @@ -29,6 +37,8 @@ [as_string, options] = get_option(options, 'as_string', false); check_unsupported_options(options, mfilename); +msgs = {}; + % Version information follows: % If the change is minor add to the last digit, large changes the middle % digit, and if the change is large and may be incompatible to previous @@ -45,7 +55,7 @@ % * Added option to make ordering of multiindices compatible with UQToolkit % * Incompatible change to 'multiindex' interface when used with more than % two arguments (removed optional 'combine' parameter) -msg = 'Attention: incompatible change in ''multiindex'' when called with more than two parameters (see help).'; +msgs{end+1} = 'Attention: incompatible change in ''multiindex'' when called with more than two parameters (see help).'; % Version 0.9.4 % * Extended SQRSPACE to cope with negative values and different exponents. @@ -56,11 +66,15 @@ % * Added gpc_integrate for integration over GPC spaces % * Added gpc_basis_eval for quick evaluation of the GPC basis functions % * Added monomial evaluation to gpc_evaluate -version = [0, 9, 5]; -%msg = ''; -% Version 0.9.6 (upcoming) -% * Cleanup of linalg +% Version 0.9.6 +% * Cleanup of linalg (esp. subspace_distance and subspace_angles) +% * Added composite and nested trapezoidal rule. +% * Improved Clenshaw-Curtis rule +msgs{end+1} = 'Attention: clenshaw_curtis_legendre_rule removed. Use clenshaw_curtis_nested instead.'; +version = [0, 9, 6]; + +% Version 0.9.7 (upcoming) % If Version information is requested as string, convert the arrary diff --git a/testing/test_fejer.m b/testing/test_fejer.m index 15fc19e1..d620890e 100644 --- a/testing/test_fejer.m +++ b/testing/test_fejer.m @@ -1,17 +1,17 @@ function test_fejer +clc +format compact +format short g + +n=1 +cc_fejer(n, 1) +cc_fejer(n, 2) +cc_fejer(n, 3) +sum(cc_fejer(n, 1)) +sum(cc_fejer(n, 2)) +sum(cc_fejer(n, 3)) -cc_fejer(5, 1) -cc_fejer(5, 2) -cc_fejer(5, 3) -sum(cc_fejer(5, 1)) -sum(cc_fejer(5, 2)) -sum(cc_fejer(5, 3)) - -[w1,w2,wc]=fejer(16); -wc(end+1)=wc(1); -sum(wc) -length(wc) function w = cc_fejer(n, mode) diff --git a/testing/test_grid.m b/testing/test_grid.m index d6cb305e..0f01db2e 100644 --- a/testing/test_grid.m +++ b/testing/test_grid.m @@ -2,7 +2,7 @@ format short g clf -f = @clenshaw_curtis_legendre_rule; +f = @clenshaw_curtis_nested; x = smolyak_grid(2, 5, f); size(x) plot(x(1,:), x(2,:), 'k.'); diff --git a/testing/test_nd_exact.m b/testing/test_nd_exact.m index 0ab39af3..31f040d2 100644 --- a/testing/test_nd_exact.m +++ b/testing/test_nd_exact.m @@ -7,7 +7,7 @@ subplot(3,4,p+4) test(13, p, @gauss_legendre_rule, 'full_tensor') subplot(3,4,p+8) - test(13, p, @clenshaw_curtis_legendre_rule, 'smolyak') + test(13, p, @clenshaw_curtis_nested, 'smolyak') end function test(pt, p, rule_func, grid) diff --git a/testing/testing_sparse2.m b/testing/testing_sparse2.m index 01def6e1..e4d9ab54 100644 --- a/testing/testing_sparse2.m +++ b/testing/testing_sparse2.m @@ -2,13 +2,13 @@ l = 4; -[x, w] = full_tensor_grid(d, l, @clenshaw_curtis_legendre_rule); +[x, w] = full_tensor_grid(d, l, @clenshaw_curtis_nested); subplot(2,2,1) plot3(x(1,:), x(2,:), x(3,:), '.') axis square -[x, w] = smolyak_grid(d, l, @clenshaw_curtis_legendre_rule); +[x, w] = smolyak_grid(d, l, @clenshaw_curtis_nested); subplot(2,2,2) plot3(x(1,:), x(2,:), x(3,:), '.') axis square