Home > K25104 > Integration > composite_rule.m

composite_rule

PURPOSE ^

Splits the interval [a,b] into n equal subintervals, then numerically

SYNOPSIS ^

function [ Q ] = composite_rule( f,a,b,n,m,rule )

DESCRIPTION ^

 Splits the interval [a,b] into n equal subintervals, then numerically
 integrates f on each interval according to the rule prescribed
 Input arguments:
   f, function handle
   a,b, integration bounds
   n, number of subintervals
   m, number of abscissae (if applicable)
   rule, rule for numerically integrating each subinterval
     possible values of rule are: 'midpoint', 'trapezium', 'simpsons', 
     'simpsons3/8', 'booles', 'milnes', 'gauss_legendre', 'gauss_lobatto' 
     and 'gauss_radau'.
 Output arguments:
   Q, value of integral

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

0001 function [ Q ] = composite_rule( f,a,b,n,m,rule )
0002 % Splits the interval [a,b] into n equal subintervals, then numerically
0003 % integrates f on each interval according to the rule prescribed
0004 % Input arguments:
0005 %   f, function handle
0006 %   a,b, integration bounds
0007 %   n, number of subintervals
0008 %   m, number of abscissae (if applicable)
0009 %   rule, rule for numerically integrating each subinterval
0010 %     possible values of rule are: 'midpoint', 'trapezium', 'simpsons',
0011 %     'simpsons3/8', 'booles', 'milnes', 'gauss_legendre', 'gauss_lobatto'
0012 %     and 'gauss_radau'.
0013 % Output arguments:
0014 %   Q, value of integral
0015 
0016 % first check user inputs
0017 if isa(f,'function_handle')==0;
0018     error('input must be a function handle');
0019 elseif a>=b;
0020     error('invalid interval');
0021 elseif mod(n,1)~=0 || n<=0; % if n does not equal zero modulo 1 then it
0022                             % is not an integer value
0023     error('input must be a positive integer');
0024 end
0025 
0026 % split interval into n equal subintervals
0027 x=linspace(a,b,n+1);
0028 
0029 y=zeros(n,1);% initialize vector to hold contribution from each subinterval
0030 if strcmp(rule,'midpoint');
0031     for i=1:n;
0032         y(i)=midpoint_rule(f,x(i),x(i+1));
0033     end
0034 elseif strcmp(rule,'trapezium');
0035     for i=1:n;
0036         y(i)=trapezium_rule(f,x(i),x(i+1));
0037     end
0038 elseif strcmp(rule,'simpsons');
0039     for i=1:n;
0040         y(i)=simpsons_rule(f,x(i),x(i+1));
0041     end    
0042 elseif strcmp(rule,'simpsons3/8');
0043     for i=1:n;
0044         y(i)=simpsons_alt_rule(f,x(i),x(i+1));
0045     end
0046 elseif strcmp(rule,'booles');
0047     for i=1:n;
0048         y(i)=booles_rule(f,x(i),x(i+1));
0049     end
0050 elseif strcmp(rule,'milnes');
0051     for i=1:n;
0052         y(i)=milnes_rule(f,x(i),x(i+1));
0053     end
0054 elseif strcmp(rule,'gauss_legendre');
0055     if mod(m,1)~=0||m<=0; 
0056         error('input must be a positive integer');
0057     end
0058     for i=1:n;
0059         y(i)=gauss_legendre_rule(f,x(i),x(i+1),m);
0060     end
0061 elseif strcmp(rule,'gauss_lobatto');
0062     if mod(m,1)~=0||m<=0; 
0063         error('input must be a positive integer');
0064     end
0065     for i=1:n;
0066         y(i)=gauss_lobatto_rule(f,x(i),x(i+1),m);
0067     end
0068 elseif strcmp(rule,'gauss_radau');
0069     if mod(m,1)~=0||m<=0; 
0070         error('input must be a positive integer');
0071     end
0072     for i=1:n;
0073         y(i)=gauss_radau_rule(f,x(i),x(i+1),m);
0074     end
0075 else 
0076     error('please enter a valid rule');
0077 end
0078 
0079 Q=sum(y);   % sum contribution from each subinterval to get y
0080 
0081 end
0082

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