From 13ea4211371744a68a8fa30384e1748191eba65c Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?St=C3=A9phane=20Adjemian=20=28Charybdis=29?= Date: Wed, 28 Nov 2012 16:15:35 +0100 Subject: [PATCH] Added routine computing Gaussian cubature weights and nodes (implementation of algorithms described in Stroud 1971). --- matlab/cubature_with_gaussian_weight.m | 314 +++++++++++++++++++++++++ 1 file changed, 314 insertions(+) create mode 100644 matlab/cubature_with_gaussian_weight.m diff --git a/matlab/cubature_with_gaussian_weight.m b/matlab/cubature_with_gaussian_weight.m new file mode 100644 index 000000000..c0f396144 --- /dev/null +++ b/matlab/cubature_with_gaussian_weight.m @@ -0,0 +1,314 @@ +function [nodes, weights] = cubature_with_gaussian_weight(d,n,method) + +%@info: +%! @deftypefn {Function File} {@var{nodes}, @var{weights} =} cubature_with_gaussian_weight (@var{d}, @var{n}) +%! @anchor{cubature_with_gaussian_weight} +%! @sp 1 +%! Computes nodes and weights for a n-order cubature with gaussian weight. +%! @sp 2 +%! @strong{Inputs} +%! @sp 1 +%! @table @ @var +%! @item d +%! Scalar integer, dimension of the region of integration. +%! @item n +%! Scalar integer equal to 3 or 5, approximation order. +%! @end table +%! @sp 2 +%! @strong{Outputs} +%! @sp 1 +%! @table @ @var +%! @item nodes +%! n*m matrix of doubles, the m nodes where the integrated function has to be evaluated. The number of nodes, m, is equal to 2*@var{d} is @var{n}==3 or 2*@var{d}^2+1 if @var{n}==5. +%! @item weights +%! m*1 vector of doubles, weights associated to the nodes. +%! @end table +%! @sp 2 +%! @strong{Remarks} +%! @sp 1 +%! The routine returns nodes and associated weights to compute a multivariate integral of the form: +%! +%! \int_D f(x)*\exp(-) dx +%! +%! +%! @end deftypefn +%@eod: + +% Copyright (C) 2012 Dynare Team +% +% This file is part of Dynare. +% +% Dynare 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. +% +% Dynare is distributed in the hope that it will be useful, +% but WITHOUT ANY WARRANTY; without even the implied warranty of +% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +% GNU General Public License for more details. +% +% You should have received a copy of the GNU General Public License +% along with Dynare. If not, see . + +% AUTHOR(S) stephane DOT adjemian AT univ DASH lemans DOT fr + +% Set default. +if nargin<3 || isempty(method) + method = 'Stroud'; +end + +if strcmp(method,'Stroud') && isequal(n,3) + %V = sqrt(pi)^d; + r = sqrt(d/2); + nodes = r*[eye(d),-eye(d)]; + weights = ones(2*d,1)/(2*d); + return +end + +if strcmp(method,'Stroud') && isequal(n,5) + r = sqrt((d+2)/2); + s = sqrt((d+2)/4); + m = 2*d^2+1; + A = 2/(n+2); + B = (4-d)/(2*(n+2)^2); + C = 1/(n+2)^2; + % Initialize the outputs + nodes = zeros(d,m); + weights = zeros(m,1); + % Set the weight for the first node (0) + weights(1) = A; + skip = 1; + % Set the remaining nodes and associated weights. + nodes(:,skip+(1:d)) = r*eye(d); + weights(skip+(1:d)) = B; + skip = skip+d; + nodes(:,skip+(1:d)) = -r*eye(d); + weights(skip+(1:d)) = B; + skip = skip+d; + for i=1:d-1 + for j = i+1:d + nodes(:,skip+(1:4)) = s*ee(d,i,j); + weights(skip+(1:4)) = C; + skip = skip+4; + end + end + return +end + +if strcmp(method,'Stroud') + error(['cubature_with_gaussian_weight:: Cubature (Stroud tables) is not yet implemented with n = ' int2str(n) '!']) +end + + + + +function v = e(n,i) + v = zeros(n,1); + v(i) = 1; + +function m = ee(n,i,j) + m = zeros(n,4); + m(:,1) = e(n,i)+e(n,j); + m(:,2) = e(n,i)-e(n,j); + m(:,3) = -m(:,2); + m(:,4) = -m(:,1); + +%@test:1 +%$ % Set problem +%$ d = 4; +%$ +%$ t = zeros(5,1); +%$ +%$ % Call the tested routine +%$ try +%$ [nodes,weights] = cubature_with_gaussian_weight(d,3); +%$ t(1) = 1; +%$ catch exception +%$ t = t(1); +%$ T = all(t); +%$ LOG = getReport(exception,'extended'); +%$ return +%$ end +%$ +%$ % Check the results. +%$ nodes = sqrt(2)*nodes; +%$ +%$ % Compute (approximated) first order moments. +%$ m1 = nodes*weights; +%$ +%$ % Compute (approximated) second order moments. +%$ m2 = nodes.^2*weights; +%$ +%$ % Compute (approximated) third order moments. +%$ m3 = nodes.^3*weights; +%$ +%$ % Compute (approximated) fourth order moments. +%$ m4 = nodes.^4*weights; +%$ +%$ t(2) = dyn_assert(m1,zeros(d,1),1e-12); +%$ t(3) = dyn_assert(m2,ones(d,1),1e-12); +%$ t(4) = dyn_assert(m3,zeros(d,1),1e-12); +%$ t(5) = dyn_assert(m4,d*ones(d,1),1e-10); +%$ T = all(t); +%@eof:1 + +%@test:2 +%$ % Set problem +%$ d = 4; +%$ Sigma = diag(1:d); +%$ Omega = diag(sqrt(1:d)); +%$ +%$ t = zeros(5,1); +%$ +%$ % Call the tested routine +%$ try +%$ [nodes,weights] = cubature_with_gaussian_weight(d,3); +%$ t(1) = 1; +%$ catch exception +%$ t = t(1); +%$ T = all(t); +%$ LOG = getReport(exception,'extended'); +%$ return +%$ end +%$ +%$ % Check the results. +%$ nodes = sqrt(2)*Omega*nodes; +%$ +%$ % Compute (approximated) first order moments. +%$ m1 = nodes*weights; +%$ +%$ % Compute (approximated) second order moments. +%$ m2 = nodes.^2*weights; +%$ +%$ % Compute (approximated) third order moments. +%$ m3 = nodes.^3*weights; +%$ +%$ % Compute (approximated) fourth order moments. +%$ m4 = nodes.^4*weights; +%$ +%$ t(2) = dyn_assert(m1,zeros(d,1),1e-12); +%$ t(3) = dyn_assert(m2,transpose(1:d),1e-12); +%$ t(4) = dyn_assert(m3,zeros(d,1),1e-12); +%$ t(5) = dyn_assert(m4,d*transpose(1:d).^2,1e-10); +%$ T = all(t); +%@eof:2 + +%@test:3 +%$ % Set problem +%$ d = 4; +%$ Sigma = diag(1:d); +%$ Omega = diag(sqrt(1:d)); +%$ +%$ t = zeros(4,1); +%$ +%$ % Call the tested routine +%$ try +%$ [nodes,weights] = cubature_with_gaussian_weight(d,3); +%$ t(1) = 1; +%$ catch exception +%$ t = t(1); +%$ T = all(t); +%$ LOG = getReport(exception,'extended'); +%$ return +%$ end +%$ +%$ % Check the results. +%$ nodes = sqrt(2)*Omega*nodes; +%$ +%$ % Compute (approximated) first order moments. +%$ m1 = nodes*weights; +%$ +%$ % Compute (approximated) second order moments. +%$ m2 = bsxfun(@times,nodes,transpose(weights))*transpose(nodes); +%$ +%$ t(2) = dyn_assert(m1,zeros(d,1),1e-12); +%$ t(3) = dyn_assert(diag(m2),transpose(1:d),1e-12); +%$ t(4) = dyn_assert(m2(:),vec(diag(diag(m2))),1e-12); +%$ T = all(t); +%@eof:3 + +%@test:4 +%$ % Set problem +%$ d = 10; +%$ a = randn(d,2*d); +%$ Sigma = a*a'; +%$ Omega = chol(Sigma,'lower'); +%$ +%$ t = zeros(4,1); +%$ +%$ % Call the tested routine +%$ try +%$ [nodes,weights] = cubature_with_gaussian_weight(d,3); +%$ t(1) = 1; +%$ catch exception +%$ t = t(1); +%$ T = all(t); +%$ LOG = getReport(exception,'extended'); +%$ return +%$ end +%$ +%$ % Correct nodes for the covariance matrix +%$ for i=1:length(weights) +%$ nodes(:,i) = sqrt(2)*Omega*nodes(:,i); +%$ end +%$ +%$ % Check the results. +%$ +%$ % Compute (approximated) first order moments. +%$ m1 = nodes*weights; +%$ +%$ % Compute (approximated) second order moments. +%$ m2 = bsxfun(@times,nodes,transpose(weights))*transpose(nodes); +%$ +%$ % Compute (approximated) third order moments. +%$ m3 = nodes.^3*weights; +%$ +%$ t(2) = dyn_assert(m1,zeros(d,1),1e-12); +%$ t(3) = dyn_assert(m2(:),vec(Sigma),1e-12); +%$ t(4) = dyn_assert(m3,zeros(d,1),1e-12); +%$ T = all(t); +%@eof:4 + +%@test:5 +%$ % Set problem +%$ d = 5; +%$ +%$ t = zeros(6,1); +%$ +%$ % Call the tested routine +%$ try +%$ [nodes,weights] = cubature_with_gaussian_weight(d,5); +%$ t(1) = 1; +%$ catch exception +%$ t = t(1); +%$ T = all(t); +%$ LOG = getReport(exception,'extended'); +%$ return +%$ end +%$ +%$ % Check the results. +%$ nodes = sqrt(2)*nodes; +%$ +%$ % Compute (approximated) first order moments. +%$ m1 = nodes*weights; +%$ +%$ % Compute (approximated) second order moments. +%$ m2 = nodes.^2*weights; +%$ +%$ % Compute (approximated) third order moments. +%$ m3 = nodes.^3*weights; +%$ +%$ % Compute (approximated) fourth order moments. +%$ m4 = nodes.^4*weights; +%$ +%$ % Compute (approximated) fifth order moments. +%$ m5 = nodes.^5*weights; +%$ +%$ t(2) = dyn_assert(m1,zeros(d,1),1e-12); +%$ t(3) = dyn_assert(m2,ones(d,1),1e-12); +%$ t(4) = dyn_assert(m3,zeros(d,1),1e-12); +%$ t(5) = dyn_assert(m4,3*ones(d,1),1e-12); +%$ t(6) = dyn_assert(m5,zeros(d,1),1e-12); +%$ T = all(t); +%@eof:5