Employs Gauss-Legendre rule to integrate f between a and b Input arguments: f, function handle a,b, integration bounds, a<b n, number of abscissae Output arguments: Q, value of integral
0001 function [ Q ] = gauss_legendre_rule( f,a,b,n ) 0002 % Employs Gauss-Legendre rule to integrate f between a and b 0003 % Input arguments: 0004 % f, function handle 0005 % a,b, integration bounds, a<b 0006 % n, number of abscissae 0007 % Output arguments: 0008 % Q, value of integral 0009 0010 % first check user inputs 0011 if isa(f,'function_handle')==0; 0012 error('input must be a function handle'); 0013 elseif a>=b; 0014 error('invalid interval'); 0015 elseif mod(n,1)~=0 || n<=0; % if n does not equal zero modulo 1 then it 0016 % is not an integer value 0017 error('input must be a positive integer'); 0018 end 0019 0020 syms x; % create a symbolic variable 0021 legendre=legendreP(n,x); % look up the n-th Legendre polynomial 0022 legendre=sym2poly(legendre);% convert symbolic expression to vector of 0023 % polynomial coefficents 0024 nodes=roots(legendre); % roots of the Legendre polynomial 0025 0026 % generate the weights by integrating the 0027 % Lagrange interpolating polynomials 0028 weights = zeros(n); 0029 syms L; % symbolic variable to hold Lagrange polynomial 0030 for k=1:n; 0031 L = Lagrangecardinal( nodes,k ); 0032 weights(k) = int(L,x,-1,1); 0033 end 0034 0035 Q=0; % initialize Q 0036 0037 % evaluate the quadrature 0038 % shifting the nodes to the interval [a,b] 0039 nodes = (b-a)/2* nodes + (a+b)/2 * ones(n,1); 0040 for i=1:n; 0041 Q=Q+weights(i)*feval(f,nodes(i)); 0042 end 0043 % adjusting the weights by a factor of (b-a)/2 0044 Q=(b-a)/2*Q; 0045 0046 end 0047