simpadpt.m

Integraton routines in matlab

function [integ, xfv, nfv] = simpadpt(fun, a, b, tol, hmin);

integ=0; 
level=0; 
i=1; 
alpha(i)= a; 
beta(i) = b; 
step    = (beta(i)-alpha(i))/4; 
nfv     = 0;
for k=1:5, 
    x=a+(k-1)*step; 
    f(i,k)=feval(fun,x); 
    nfv=nfv+1; 
end
while (i > 0)
   % Compute Simson formula on interval [alpha(i), beta(i)]
   h = (beta(i)-alpha(i))/2; 
   S = (h/3)*(f(i,1)+4*f(i,3)+f(i,5));
   % Compute composite Simson formula on interval 
   % [alpha(i), (alpha(i)+beta(i))/2], [(alpha(i)+beta(i))/2, beta(i)]
   h  = h/2; 
   S2 = (h/3)*(f(i,1)+4*f(i,2)+f(i,3)); 
   S2 = S2+(h/3)*(f(i,3)+4*f(i,4)+f(i,5));
   
   tolrv = tol*(beta(i)-alpha(i))/(b-a); 
   errrv = abs(S-S2)/10;
   if (errrv > tolrv)
      % integral on [alpha(i), beta(i)] is not approximated accurately enough
      % Reduce interval to [alpha(i), (alpha(i)+beta(i))/2]
      i=i+1; 
      alpha(i) = alpha(i-1); 
      beta(i)  = (alpha(i-1)+beta(i-1))/2;
      f(i,1)   = f(i-1,1);
      f(i,3)   = f(i-1,2);
      f(i,5)   = f(i-1,3);
      len=abs(beta(i)-alpha(i));
      if (len >= hmin)
         step=len/4; 
         x=alpha(i)+step; 
         f(i,2)=feval(fun, x); 
         nfv=nfv+1; 
         x=beta(i)-step; 
         f(i,4)=feval(fun,x); 
         nfv=nfv+1;
      else, 
         xfv=xfv'; 
         disp(' Too small steplength   ')
         return 
      end
   else    
      % integral on [alpha(i), beta(i)] is approximated accurately enough
      % Work on next inteval [beta(i), beta(i-1)]
      integ=integ+S2; 
      level=level+1; 
      if (level==1)
         for k=1:5, 
             xfv(k)=alpha(i)+(k-1)*h; 
         end; 
         ist=5;
      else
         for k=1:4, 
             xfv(ist+k)=alpha(i)+k*h; 
         end; 
         ist=ist+4; 
      end;   
      if (beta(i)==b)
         return
      end; 
      i=i-1; 
      alpha(i)=beta(i+1); 
      f(i,1)=f(i+1,5); 
      f(i,3)=f(i,4); 
      step=abs(beta(i)-alpha(i))/4; 
      x=alpha(i)+step; 
      f(i,2)=feval(fun,x); 
      nfv=nfv+1; 
      x=beta(i)-step; 
      f(i,4)=feval(fun,x); 
      nfv=nfv+1;
   end
end