Home > K25104 > Integration > gauss_lobatto_rule.m

gauss_lobatto_rule

PURPOSE ^

Employs Gauss-Lobatto rule with Legendre polynomials

SYNOPSIS ^

function [ Q ] = gauss_lobatto_rule( f,a,b,n )

DESCRIPTION ^

 Employs Gauss-Lobatto rule with Legendre polynomials 
 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

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

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

Generated on Mon 18-Jan-2016 10:25:49 by m2html © 2005