%  scripts and functions for MAT 331, Thursday, Sept 14, 2017

%-------------------------------------------------------------
% function to compute different approximations for integrals:
% left-hand-rule, right-hand-rule, trapezoid rule and midpoint rule

function  out = integrate(f,a,b,n)
   delta=(b-a)/n;
   x=[a:delta:b];
   x1=x(1:end-1);
   x2=x(2:end);
   lhr=delta*sum(f(x1));
   rhr=delta*sum(f(x2));
   trap=(lhr+rhr)/2;
   mid=delta*sum(f((x1+x2)/2));
   % n should be even for Simpson's rule
   simp=(delta/3)*(4*sum(f(x(2:2:end)))+...
      2*sum(f(x(1:2:end)))-f(x(1))-f(x(end)));
   out=[lhr,rhr,trap,mid,simp]
return


%-------------------------------------------------------------

% script to compare accuracy of different methods of integration

f = @(x) x.^2;    % choose function to test
answer = 1/3;     % give exact answer (from calculus)
a=0;              % choose integration interval
b=1;

count=1;
clear results;
test=[10:10:100]   % choose number of divisions to test
for  n=test        
    results(count,:)=integrate(f,a,b,n);
    count=count+1;
end % for n
results

figure;
hold on;
grid on;
title('Plot results for different methods of integration');
plot(test,results(:,1),'b');   % plot left hand rule
plot(test,results(:,2),'r');   % plot right hand rule
plot(test,results(:,3),'k');   % plot trapezoid rule
plot(test,results(:,4),'m');   % plot midpoint rule
plot(test,results(:,4),'g');   % plot Simpson's rule
legend('LHR','RHR','TRAP','MID','SIMP'); % creates labels on plot

%return;

accuracy =abs(results-answer)
figure;
hold on;
grid on;
title('Plot accuracy for different methods of integration');
plot(test,accuracy(:,1),'b');   % plot left hand rule
plot(test,accuracy(:,2),'r');   % plot right hand rule
plot(test,accuracy(:,3),'k');   % plot trapezoid rule
plot(test,accuracy(:,4),'c');   % plot midpoint rule
plot(test,accuracy(:,4),'g');   % plot Simpson's rule
legend('LHR','RHR','TRAP','MID','SIMP'); % creates labels on plot

%return;

log_accuracy=log(accuracy)
figure;
hold on;
grid on;
title('Plot log-accuracy for different methods of integration');
plot(test,log_accuracy(:,1),'b');   % plot left hand rule
plot(test,log_accuracy(:,2),'r');   % plot right hand rule
plot(test,log_accuracy(:,3),'k');   % plot trapezoid rule
plot(test,log_accuracy(:,4),'c');   % plot midpoint rule
plot(test,log_accuracy(:,5),'c');   % plot Simpson's rule
legend('LHR','RHR','TRAP','MID','SIMP'); % creates labels on plot