Home > K25104 > Integration > left_gauss_radau_rule.m

left_gauss_radau_rule

PURPOSE ^

Employs the left Gauss-Radau rule with n abscissae to integrate f between a and b

SYNOPSIS ^

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

DESCRIPTION ^

 Employs the left Gauss-Radau rule with n abscissae 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 ] = left_gauss_radau_rule( f,a,b,n )
0002 % Employs the left Gauss-Radau rule with n abscissae 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 leg1=legendreP(n,x);       % look up the necessary legendre polynomials
0022 leg0=legendreP(n-1,x);
0023 
0024 nodes=(leg0+leg1);  % construct function whose roots give the nodes
0025 nodes=sym2poly(nodes);     % convert symbolic expression to polynomial
0026 nodes=roots(nodes);        % find roots
0027 [~,index]=min(nodes);    % find index of element in nodes which equals -1
0028 % remove the end point entry, to get interior abscissae
0029 nodes=[nodes(1:(index-1));nodes((index+1):n)];  
0030 leg0=sym2poly(leg0);    % convert symbolic expression to polynomial
0031 
0032 weights=zeros(n-1,1);   % initialize vector to hold weights
0033 for k=1:n-1;
0034     temp=polyval(leg0,nodes(k)); % evaluate legendre polynomial at the knot
0035     weights(k)=(1-nodes(k))/((n*temp)^2);   % calculate weights
0036 end
0037 
0038 % evaluate the quadrature
0039 Q=2*feval(f,a)/(n*n);  % end point contribution
0040 
0041 for i=1:n-1;    % add the contributions from interior nodes
0042     Q=Q+weights(i)*feval(f,(b-a)*nodes(i)/2+(a+b)/2);   
0043 end
0044 % adjusting the weights by a factor of (b-a)/2
0045 Q=(b-a)/2*Q;
0046 
0047 end

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