Computes the weights and nodes for an Hermite Gaussian quadrature rule.
0001 function [nodes,weights] = gauss_hermite_weights_and_nodes(n) 0002 % Computes the weights and nodes for an Hermite Gaussian quadrature rule. 0003 0004 %@info: 0005 %! @deftypefn {Function File} {@var{nodes}, @var{weights} =} gauss_hermite_weights_and_nodes (@var{n}) 0006 %! @anchor{gauss_hermite_weights_and_nodes} 0007 %! @sp 1 0008 %! Computes the weights and nodes for an Hermite Gaussian quadrature rule. designed to approximate integrals 0009 %! on the infinite interval (-\infty,\infty) of an unweighted smooth function. 0010 %! @sp 2 0011 %! @strong{Inputs} 0012 %! @sp 1 0013 %! @table @ @var 0014 %! @item n 0015 %! Positive integer scalar, number of nodes (order of approximation). 0016 %! @end table 0017 %! @sp 1 0018 %! @strong{Outputs} 0019 %! @sp 1 0020 %! @table @ @var 0021 %! @item nodes 0022 %! n*1 vector of doubles, the nodes (roots of an order n Hermite polynomial) 0023 %! @item weights 0024 %! n*1 vector of doubles, the associated weights. 0025 %! @end table 0026 %! @sp 2 0027 %! @strong{This function is called by:} 0028 %! @sp 2 0029 %! @strong{This function calls:} 0030 %! @sp 2 0031 %! @end deftypefn 0032 %@eod: 0033 0034 % Copyright (C) 2011 Dynare Team 0035 % stephane DOT adjemian AT univ DASH lemans DOT fr 0036 % 0037 % This file is part of Dynare. 0038 % 0039 % Dynare is free software: you can redistribute it and/or modify 0040 % it under the terms of the GNU General Public License as published by 0041 % the Free Software Foundation, either version 3 of the License, or 0042 % (at your option) any later version. 0043 % 0044 % Dynare is distributed in the hope that it will be useful, 0045 % but WITHOUT ANY WARRANTY; without even the implied warranty of 0046 % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 0047 % GNU General Public License for more details. 0048 % 0049 % You should have received a copy of the GNU General Public License 0050 % along with Dynare. If not, see <http://www.gnu.org/licenses/>. 0051 0052 b = sqrt([1:n-1]/2); 0053 JacobiMatrix = diag(b,1)+diag(b,-1); 0054 [JacobiEigenVectors,JacobiEigenValues] = eig(JacobiMatrix); 0055 [nodes,idx] = sort(diag(JacobiEigenValues)); 0056 JacobiEigenVector = JacobiEigenVectors(1,:); 0057 JacobiEigenVector = transpose(JacobiEigenVector(idx)); 0058 weights = JacobiEigenVector.^2; 0059 nodes = sqrt(2)*nodes; 0060 0061 %@test:1 0062 %$ n = 5; 0063 %$ [nodes,weights] = gauss_hermite_weights_and_nodes(n); 0064 %$ 0065 %$ sum_of_weights = sum(weights); 0066 %$ 0067 %$ % Expected nodes (taken from Judd (1998, table 7.4). 0068 %$ enodes = [-2.020182870; -0.9585724646; 0; 0.9585724646; 2.020182870]; 0069 %$ 0070 %$ % Check the results. 0071 %$ t(1) = dyn_assert(1.0,sum_of_weights,1e-12); 0072 %$ t(2) = dyn_assert(enodes,nodes/sqrt(2),1e-8); 0073 %$ T = all(t); 0074 %@eof:1 0075 0076 %@test:2 0077 %$ n = 9; 0078 %$ [nodes,weights] = gauss_hermite_weights_and_nodes(n); 0079 %$ 0080 %$ sum_of_weights = sum(weights); 0081 %$ expectation = sum(weights.*nodes); 0082 %$ variance = sum(weights.*(nodes.^2)); 0083 %$ 0084 %$ % Check the results. 0085 %$ t(1) = dyn_assert(1.0,sum_of_weights,1e-12); 0086 %$ t(2) = dyn_assert(1.0,variance,1e-12); 0087 %$ t(3) = dyn_assert(0.0,expectation,1e-12); 0088 %$ T = all(t); 0089 %@eof:2 0090 0091 %@test:3 0092 %$ n = 9; 0093 %$ [nodes,weights] = gauss_hermite_weights_and_nodes(n); 0094 %$ 0095 %$ NODES = cartesian_product_of_sets(nodes,nodes); 0096 %$ WEIGHTS = cartesian_product_of_sets(weights,weights); 0097 %$ WEIGHTS = prod(WEIGHTS,2); 0098 %$ 0099 %$ sum_of_weights = sum(WEIGHTS); 0100 %$ expectation = transpose(WEIGHTS)*NODES; 0101 %$ variance = transpose(WEIGHTS)*NODES.^2; 0102 %$ 0103 %$ % Check the results. 0104 %$ t(1) = dyn_assert(1.0,sum_of_weights,1e-12); 0105 %$ t(2) = dyn_assert(ones(1,2),variance,1e-12); 0106 %$ t(3) = dyn_assert(zeros(1,2),expectation,1e-12); 0107 %$ T = all(t); 0108 %@eof:3 0109 0110 %@test:4 0111 %$ n = 9; sigma = .1; 0112 %$ [nodes,weights] = gauss_hermite_weights_and_nodes(n); 0113 %$ 0114 %$ sum_of_weights = sum(weights); 0115 %$ expectation = sum(weights.*nodes*.1); 0116 %$ variance = sum(weights.*((nodes*.1).^2)); 0117 %$ 0118 %$ % Check the results. 0119 %$ t(1) = dyn_assert(1.0,sum_of_weights,1e-12); 0120 %$ t(2) = dyn_assert(.01,variance,1e-12); 0121 %$ t(3) = dyn_assert(0.0,expectation,1e-12); 0122 %$ T = all(t); 0123 %@eof:4