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