1function [nodes,weights] = gauss_legendre_weights_and_nodes(n,a,b) 2% Computes the weights and nodes for a Legendre Gaussian quadrature rule. 3 4%@info: 5%! @deftypefn {Function File} {@var{nodes}, @var{weights} =} gauss_hermite_weights_and_nodes (@var{n}) 6%! @anchor{gauss_legendre_weights_and_nodes} 7%! @sp 1 8%! Computes the weights and nodes for a Legendre Gaussian quadrature rule. designed to approximate integrals 9%! on the finite interval (-1,1) 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%! @item a 17%! Double scalar, lower bound. 18%! @item b 19%! Double scalar, upper bound. 20%! @end table 21%! @sp 1 22%! @strong{Outputs} 23%! @sp 1 24%! @table @ @var 25%! @item nodes 26%! n*1 vector of doubles, the nodes (roots of an order n Legendre polynomial) 27%! @item weights 28%! n*1 vector of doubles, the associated weights. 29%! @end table 30%! @sp 2 31%! @strong{Remarks:} 32%! Only the first input argument (the number of nodes) is mandatory. The second and third input arguments 33%! are used if a change of variables is necessary (ie if we need nodes over the interval [a,b] instead of 34%! of the default interval [-1,1]). 35%! @sp 2 36%! @strong{This function is called by:} 37%! @sp 2 38%! @strong{This function calls:} 39%! @sp 2 40%! @end deftypefn 41%@eod: 42 43% Copyright (C) 2012-2017 Dynare Team 44% 45% This file is part of Dynare. 46% 47% Dynare is free software: you can redistribute it and/or modify 48% it under the terms of the GNU General Public License as published by 49% the Free Software Foundation, either version 3 of the License, or 50% (at your option) any later version. 51% 52% Dynare is distributed in the hope that it will be useful, 53% but WITHOUT ANY WARRANTY; without even the implied warranty of 54% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 55% GNU General Public License for more details. 56% 57% You should have received a copy of the GNU General Public License 58% along with Dynare. If not, see <http://www.gnu.org/licenses/>. 59 60% AUTHOR(S) stephane DOT adjemian AT univ DASH lemans DOT fr 61 62bb = sqrt(1./(4-(1./transpose(1:n-1)).^2)); 63aa = zeros(n,1); 64JacobiMatrix = diag(bb,1)+diag(aa)+diag(bb,-1); 65[JacobiEigenVectors,JacobiEigenValues] = eig(JacobiMatrix); 66[nodes,idx] = sort(diag(JacobiEigenValues)); 67JacobiEigenVector = JacobiEigenVectors(1,:); 68JacobiEigenVector = transpose(JacobiEigenVector(idx)); 69weights = 2*JacobiEigenVector.^2; 70 71if nargin==3 72 weights = .5*(b-a)*weights; 73 nodes = .5*(nodes+1)*(b-a)+a; 74end 75 76%@test:1 77%$ [n2,w2] = gauss_legendre_weights_and_nodes(2); 78%$ [n3,w3] = gauss_legendre_weights_and_nodes(3); 79%$ [n4,w4] = gauss_legendre_weights_and_nodes(4); 80%$ [n5,w5] = gauss_legendre_weights_and_nodes(5); 81%$ [n7,w7] = gauss_legendre_weights_and_nodes(7); 82%$ 83%$ 84%$ % Expected nodes (taken from Judd (1998, table 7.2)). 85%$ e2 = .5773502691; e2 = [-e2; e2]; 86%$ e3 = .7745966692; e3 = [-e3; 0 ; e3]; 87%$ e4 = [.8611363115; .3399810435]; e4 = [-e4; flipud(e4)]; 88%$ e5 = [.9061798459; .5384693101]; e5 = [-e5; 0; flipud(e5)]; 89%$ e7 = [.9491079123; .7415311855; .4058451513]; e7 = [-e7; 0; flipud(e7)]; 90%$ 91%$ % Expected weights (taken from Judd (1998, table 7.2) and http://en.wikipedia.org/wiki/Gaussian_quadrature). 92%$ f2 = [1; 1]; 93%$ f3 = [5; 8; 5]/9; 94%$ f4 = [18-sqrt(30); 18+sqrt(30)]; f4 = [f4; flipud(f4)]/36; 95%$ f5 = [322-13*sqrt(70); 322+13*sqrt(70)]/900; f5 = [f5; 128/225; flipud(f5)]; 96%$ f7 = [.1294849661; .2797053914; .3818300505]; f7 = [f7; .4179591836; flipud(f7)]; 97%$ 98%$ % Check the results. 99%$ t(1) = dassert(e2,n2,1e-9); 100%$ t(2) = dassert(e3,n3,1e-9); 101%$ t(3) = dassert(e4,n4,1e-9); 102%$ t(4) = dassert(e5,n5,1e-9); 103%$ t(5) = dassert(e7,n7,1e-9); 104%$ t(6) = dassert(w2,f2,1e-9); 105%$ t(7) = dassert(w3,f3,1e-9); 106%$ t(8) = dassert(w4,f4,1e-9); 107%$ t(9) = dassert(w5,f5,1e-9); 108%$ t(10) = dassert(w7,f7,1e-9); 109%$ T = all(t); 110%@eof:1 111 112%@test:2 113%$ nmax = 50; 114%$ 115%$ t = zeros(nmax,1); 116%$ 117%$ for i=1:nmax 118%$ [n,w] = gauss_legendre_weights_and_nodes(i); 119%$ t(i) = dassert(sum(w),2,1e-12); 120%$ end 121%$ 122%$ T = all(t); 123%@eof:2 124 125%@test:3 126%$ [n,w] = gauss_legendre_weights_and_nodes(9,pi,2*pi); 127%$ % Check that the 128%$ t(1) = all(n>pi); 129%$ t(2) = all(n<2*pi); 130%$ t(3) = dassert(sum(w),pi,1e-12); 131%$ T = all(t); 132%@eof:3