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
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