1function [nodes,weights] = gauss_hermite_weights_and_nodes(n) 2% Computes the weights and nodes for an Hermite Gaussian quadrature rule. 3 4%@info: 5%! @deftypefn {Function File} {@var{nodes}, @var{weights} =} gauss_hermite_weights_and_nodes (@var{n}) 6%! @anchor{gauss_hermite_weights_and_nodes} 7%! @sp 1 8%! Computes the weights and nodes for an Hermite Gaussian quadrature rule. designed to approximate integrals 9%! on the infinite interval (-\infty,\infty) of an unweighted smooth function. 10%! @sp 2 11%! @strong{Inputs} 12%! @sp 1 13%! @table @ @var 14%! @item n 15%! Positive integer scalar, number of nodes (order of approximation). 16%! @end table 17%! @sp 1 18%! @strong{Outputs} 19%! @sp 1 20%! @table @ @var 21%! @item nodes 22%! n*1 vector of doubles, the nodes (roots of an order n Hermite polynomial) 23%! @item weights 24%! n*1 vector of doubles, the associated weights. 25%! @end table 26%! @sp 2 27%! @strong{This function is called by:} 28%! @sp 2 29%! @strong{This function calls:} 30%! @sp 2 31%! @end deftypefn 32%@eod: 33 34% Copyright (C) 2011-2017 Dynare Team 35% 36% This file is part of Dynare. 37% 38% Dynare is free software: you can redistribute it and/or modify 39% it under the terms of the GNU General Public License as published by 40% the Free Software Foundation, either version 3 of the License, or 41% (at your option) any later version. 42% 43% Dynare is distributed in the hope that it will be useful, 44% but WITHOUT ANY WARRANTY; without even the implied warranty of 45% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 46% GNU General Public License for more details. 47% 48% You should have received a copy of the GNU General Public License 49% along with Dynare. If not, see <http://www.gnu.org/licenses/>. 50 51% Original author: stephane DOT adjemian AT univ DASH lemans DOT fr 52 53b = sqrt([1:n-1]/2); 54JacobiMatrix = diag(b,1)+diag(b,-1); 55[JacobiEigenVectors,JacobiEigenValues] = eig(JacobiMatrix); 56[nodes,idx] = sort(diag(JacobiEigenValues)); 57JacobiEigenVector = JacobiEigenVectors(1,:); 58JacobiEigenVector = transpose(JacobiEigenVector(idx)); 59weights = JacobiEigenVector.^2; 60nodes = sqrt(2)*nodes; 61 62%@test:1 63%$ n = 5; 64%$ [nodes,weights] = gauss_hermite_weights_and_nodes(n); 65%$ 66%$ sum_of_weights = sum(weights); 67%$ 68%$ % Expected nodes (taken from Judd (1998, table 7.4). 69%$ enodes = [-2.020182870; -0.9585724646; 0; 0.9585724646; 2.020182870]; 70%$ 71%$ % Check the results. 72%$ t(1) = dassert(1.0,sum_of_weights,1e-12); 73%$ t(2) = dassert(enodes,nodes/sqrt(2),1e-8); 74%$ T = all(t); 75%@eof:1 76 77%@test:2 78%$ n = 9; 79%$ [nodes,weights] = gauss_hermite_weights_and_nodes(n); 80%$ 81%$ sum_of_weights = sum(weights); 82%$ expectation = sum(weights.*nodes); 83%$ variance = sum(weights.*(nodes.^2)); 84%$ 85%$ % Check the results. 86%$ t(1) = dassert(1.0,sum_of_weights,1e-12); 87%$ t(2) = dassert(1.0,variance,1e-12); 88%$ t(3) = dassert(0.0,expectation,1e-12); 89%$ T = all(t); 90%@eof:2 91 92%@test:3 93%$ n = 9; 94%$ [nodes,weights] = gauss_hermite_weights_and_nodes(n); 95%$ 96%$ NODES = cartesian_product_of_sets(nodes,nodes); 97%$ WEIGHTS = cartesian_product_of_sets(weights,weights); 98%$ WEIGHTS = prod(WEIGHTS,2); 99%$ 100%$ sum_of_weights = sum(WEIGHTS); 101%$ expectation = transpose(WEIGHTS)*NODES; 102%$ variance = transpose(WEIGHTS)*NODES.^2; 103%$ 104%$ % Check the results. 105%$ t(1) = dassert(1.0,sum_of_weights,1e-12); 106%$ t(2) = dassert(ones(1,2),variance,1e-12); 107%$ t(3) = dassert(zeros(1,2),expectation,1e-12); 108%$ T = all(t); 109%@eof:3 110 111%@test:4 112%$ n = 9; sigma = .1; 113%$ [nodes,weights] = gauss_hermite_weights_and_nodes(n); 114%$ 115%$ sum_of_weights = sum(weights); 116%$ expectation = sum(weights.*nodes*.1); 117%$ variance = sum(weights.*((nodes*.1).^2)); 118%$ 119%$ % Check the results. 120%$ t(1) = dassert(1.0,sum_of_weights,1e-12); 121%$ t(2) = dassert(.01,variance,1e-12); 122%$ t(3) = dassert(0.0,expectation,1e-12); 123%$ T = all(t); 124%@eof:4 125