Computes the weights and nodes for a Legendre Gaussian quadrature rule.
0001 function [nodes,weights] = gauss_legendre_weights_and_nodes(n,a,b) 0002 % Computes the weights and nodes for a Legendre Gaussian quadrature rule. 0003 0004 %@info: 0005 %! @deftypefn {Function File} {@var{nodes}, @var{weights} =} gauss_hermite_weights_and_nodes (@var{n}) 0006 %! @anchor{gauss_legendre_weights_and_nodes} 0007 %! @sp 1 0008 %! Computes the weights and nodes for a Legendre Gaussian quadrature rule. designed to approximate integrals 0009 %! on the finite interval (-1,1) 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 %! @item a 0017 %! Double scalar, lower bound. 0018 %! @item b 0019 %! Double scalar, upper bound. 0020 %! @end table 0021 %! @sp 1 0022 %! @strong{Outputs} 0023 %! @sp 1 0024 %! @table @ @var 0025 %! @item nodes 0026 %! n*1 vector of doubles, the nodes (roots of an order n Legendre polynomial) 0027 %! @item weights 0028 %! n*1 vector of doubles, the associated weights. 0029 %! @end table 0030 %! @sp 2 0031 %! @strong{Remarks:} 0032 %! Only the first input argument (the number of nodes) is mandatory. The second and third input arguments 0033 %! are used if a change of variables is necessary (ie if we need nodes over the interval [a,b] instead of 0034 %! of the default interval [-1,1]). 0035 %! @sp 2 0036 %! @strong{This function is called by:} 0037 %! @sp 2 0038 %! @strong{This function calls:} 0039 %! @sp 2 0040 %! @end deftypefn 0041 %@eod: 0042 0043 % Copyright (C) 2012 Dynare Team 0044 % 0045 % This file is part of Dynare. 0046 % 0047 % Dynare is free software: you can redistribute it and/or modify 0048 % it under the terms of the GNU General Public License as published by 0049 % the Free Software Foundation, either version 3 of the License, or 0050 % (at your option) any later version. 0051 % 0052 % Dynare is distributed in the hope that it will be useful, 0053 % but WITHOUT ANY WARRANTY; without even the implied warranty of 0054 % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 0055 % GNU General Public License for more details. 0056 % 0057 % You should have received a copy of the GNU General Public License 0058 % along with Dynare. If not, see <http://www.gnu.org/licenses/>. 0059 0060 % AUTHOR(S) stephane DOT adjemian AT univ DASH lemans DOT fr 0061 0062 bb = sqrt(1./(4-(1./transpose(1:n-1)).^2)); 0063 aa = zeros(n,1); 0064 JacobiMatrix = diag(bb,1)+diag(aa)+diag(bb,-1); 0065 [JacobiEigenVectors,JacobiEigenValues] = eig(JacobiMatrix); 0066 [nodes,idx] = sort(diag(JacobiEigenValues)); 0067 JacobiEigenVector = JacobiEigenVectors(1,:); 0068 JacobiEigenVector = transpose(JacobiEigenVector(idx)); 0069 weights = 2*JacobiEigenVector.^2; 0070 0071 if nargin==3 0072 weights = .5*(b-a)*weights; 0073 nodes = .5*(nodes+1)*(b-a)+a; 0074 end 0075 0076 %@test:1 0077 %$ [n2,w2] = gauss_legendre_weights_and_nodes(2); 0078 %$ [n3,w3] = gauss_legendre_weights_and_nodes(3); 0079 %$ [n4,w4] = gauss_legendre_weights_and_nodes(4); 0080 %$ [n5,w5] = gauss_legendre_weights_and_nodes(5); 0081 %$ [n7,w7] = gauss_legendre_weights_and_nodes(7); 0082 %$ 0083 %$ 0084 %$ % Expected nodes (taken from Judd (1998, table 7.2)). 0085 %$ e2 = .5773502691; e2 = [-e2; e2]; 0086 %$ e3 = .7745966692; e3 = [-e3; 0 ; e3]; 0087 %$ e4 = [.8611363115; .3399810435]; e4 = [-e4; flipud(e4)]; 0088 %$ e5 = [.9061798459; .5384693101]; e5 = [-e5; 0; flipud(e5)]; 0089 %$ e7 = [.9491079123; .7415311855; .4058451513]; e7 = [-e7; 0; flipud(e7)]; 0090 %$ 0091 %$ % Expected weights (taken from Judd (1998, table 7.2) and http://en.wikipedia.org/wiki/Gaussian_quadrature). 0092 %$ f2 = [1; 1]; 0093 %$ f3 = [5; 8; 5]/9; 0094 %$ f4 = [18-sqrt(30); 18+sqrt(30)]; f4 = [f4; flipud(f4)]/36; 0095 %$ f5 = [322-13*sqrt(70); 322+13*sqrt(70)]/900; f5 = [f5; 128/225; flipud(f5)]; 0096 %$ f7 = [.1294849661; .2797053914; .3818300505]; f7 = [f7; .4179591836; flipud(f7)]; 0097 %$ 0098 %$ % Check the results. 0099 %$ t(1) = dyn_assert(e2,n2,1e-9); 0100 %$ t(2) = dyn_assert(e3,n3,1e-9); 0101 %$ t(3) = dyn_assert(e4,n4,1e-9); 0102 %$ t(4) = dyn_assert(e5,n5,1e-9); 0103 %$ t(5) = dyn_assert(e7,n7,1e-9); 0104 %$ t(6) = dyn_assert(w2,f2,1e-9); 0105 %$ t(7) = dyn_assert(w3,f3,1e-9); 0106 %$ t(8) = dyn_assert(w4,f4,1e-9); 0107 %$ t(9) = dyn_assert(w5,f5,1e-9); 0108 %$ t(10) = dyn_assert(w7,f7,1e-9); 0109 %$ T = all(t); 0110 %@eof:1 0111 0112 %@test:2 0113 %$ nmax = 50; 0114 %$ 0115 %$ t = zeros(nmax,1); 0116 %$ 0117 %$ for i=1:nmax 0118 %$ [n,w] = gauss_legendre_weights_and_nodes(i); 0119 %$ t(i) = dyn_assert(sum(w),2,1e-12); 0120 %$ end 0121 %$ 0122 %$ T = all(t); 0123 %@eof:2 0124 0125 %@test:3 0126 %$ [n,w] = gauss_legendre_weights_and_nodes(9,pi,2*pi); 0127 %$ % Check that the 0128 %$ t(1) = all(n>pi); 0129 %$ t(2) = all(n<2*pi); 0130 %$ t(3) = dyn_assert(sum(w),pi,1e-12); 0131 %$ T = all(t); 0132 %@eof:3