integrate2

 higher order integration

0001 function [s]=integrate2(y,dx)
0002 %
0003 % higher order integration
0004 %
0005 %
0006 
0007 
0008 N=floor( (length(y)-5)/4);
0009 
0010 
0011 
0012 for (i=0:N);
0013   s(i+1)=2*dx/45*( 7*y(4*i+1) + 32*y(4*i+2) + 12*y(4*i+3) + 32*y(4*i+4) ...
0014          + 7*y(4*i+5));
0015 end
0016 
0017 last_grid_pt=4*N+5;
0018 L=length(y)-last_grid_pt;
0019 
0020 last_grid_pt,
0021 
0022 I=0;
0023 switch L
0024  case 1
0025   disp('one extra point');
0026  
0027   i=last_grid_pt-3;
0028   tx=0:4;
0029   ty=y(i:(i+4));
0030   p=polyfit(tx,ty,4);
0031   
0032   u=p(1)*4^5/5 + p(2)*4^4/4 + p(3)*4^3/3 + p(4)*4^2/2 + p(5)*4;
0033   l=p(1)*3^5/5 + p(2)*3^4/4 + p(3)*3^3/3 + p(4)*3^2/2 + p(5)*3;
0034   
0035 
0036   I=(u-l)*dx;
0037   
0038 
0039  case 2
0040   disp('two extra points');
0041 
0042   
0043 
0044   i=last_grid_pt-2;
0045   tx=0:4;
0046   ty=y(i:(i+4));
0047   p=polyfit(tx,ty,4);
0048   
0049   u=p(1)*4^5/5 + p(2)*4^4/4 + p(3)*4^3/3 + p(4)*4^2/2 + p(5)*4;
0050   l=p(1)*2^5/5 + p(2)*2^4/4 + p(3)*2^3/3 + p(4)*2^2/2 + p(5)*2;
0051  
0052   
0053   I=(u-l)*dx;
0054   
0055  case 3
0056   disp('three extra points');
0057 
0058   
0059 
0060   i=last_grid_pt-1;
0061   tx=0:4;
0062   ty=y(i:(i+4));
0063   p=polyfit(tx,ty,4);
0064   
0065   u=p(1)*4^5/5 + p(2)*4^4/4 + p(3)*4^3/3 + p(4)*4^2/2 + p(5)*4;
0066   l=p(1)*1^5/5 + p(2)*1^4/4 + p(3)*1^3/3 + p(4)*1^2/2 + p(5)*1;
0067  
0068   
0069   I=(u-l)*dx;
0070 
0071 
0072   
0073 end;
0074 
0075   
0076 
0077 s=sum(s)+I;