% Script 3 for Week 8

% We test here the cumulative integral functions like cumtrapz but
% using Simpsons rule and allowing  unequal data spacing. 

% They are both WAY better than MATLABs cumtrapz. 

% You can copy and run the code immediately beneath as a block several times
%  to see how cumtrapz, cumsimp13 and cumsimp38 behave relative to each other, 
% at least for this particular test function.


format long

clear all
t=10*randn([1,34]);
t=sort(t);
max(abs(diff(t)))

f=@(t) t.^2+sin(t)
y=f(t);

% y=yf+0.1*randn([1,34]);

for k=1:34
exact(k)=integral(@(x) f(x), t(1),t(k));
end

A=cumtrapz(t,y);

D=cumsimp13(t,y);

E=cumsimp38(t,y);

et=abs(exact-A);
es13=abs(exact-D);
es38=abs(exact-E);

met=max(et)
mes38=max(es38)
mes13=max(es13)

setr=sum(et)
ses13=sum(es13)
ses38=sum(es38)


figure(43)
plot(t,et,'*',t,es13,'*',t,es38,'*')
legend('et','es13','es38')


figure(44)
plot(t,es13,'*',t,es38,'*')
legend('es13','es38')

tic
for k=1:100
cumtrapz(t,y);
end
cumtrapzspeed=toc


tic
for k=1:100
cumsimp13(t,y);
end
cumsimp13speed=toc

tic
for k=1:100
cumsimp38(t,y);
end
cumsimp38speed=toc




% Here is another way of testing these tools using a function 
% I call TheComparalator.m

%*****************************************************
function [E R]=TheComparalator(fun,a,b)

% This script compares the order of the error of various numerical 
% integration functions. Input function to be integrated and interval.
% 

exact=integral(fun,a,b);
if exact==0,error('The integral over the selected interval is 0. Choose another interval.'),end

fprintf('\n The rows provide performance stats (log base 6 of the error)  for:\n trapuneq.m, OneThirdSimp.m, simp13.m, ThreeEighthSimp.m and simp38.m\n respectively and the columns test these methods by breaking\n the interval into 12, 72, 432, 2592 and finally 15552 sub-intervals.\n') 

for k=1:5
s=linspace(a,b,2*6^k+1);
ys=fun(s);
A1(k)=trapuneq(s,ys);
A2(k)=OneThirdSimp(s,ys);
A3(k)=simp13(s,ys);
A4(k)=ThreeEighthSimp(s,ys);
A5(k)=simp38(s,ys);

E(1,k)=log(abs((A1(k)-exact)/exact))/log(6);
E(2,k)=log(abs((A2(k)-exact)/exact))/log(6);
E(3,k)=log(abs((A3(k)-exact)/exact))/log(6);
E(4,k)=log(abs((A4(k)-exact)/exact))/log(6);
E(5,k)=log(abs((A5(k)-exact)/exact))/log(6);

end
%*****************************************************

If you test it by

TheComparalator(@(t) sin(t)+t.^4+1./t, 0.2,3)

you get:

 The rows provide performance stats (log base 6 of the error)  for:
 trapuneq.m, OneThirdSimp.m, simp13.m, ThreeEighthSimp.m and simp38.m
 respectively and the columns test these methods by breaking
 the interval into 12, 72, 432, 2592 and finally 15552 sub-intervals.

ans =

   -2.5199   -4.5100   -6.5097   -8.5096  -10.5096
   -4.3731   -7.8072  -11.7641  -15.7625  -16.8613
   -4.0597   -7.1876  -10.9478  -14.9027  -18.1543
   -4.1228   -7.3910  -11.3128  -15.3105  -16.8544
   -4.1356   -7.3677  -11.2875  -15.3021  -18.4166
   
   The actual integral we are approximating is
   
   integral(@(t) sin(t)+t.^4+1./t, 0.2,3)
   ans = 53.2780

and  eps(ans) =  7.1054e-15

The table contains the exponent on 6 of the error of the method. The last column corresponds to a delta x of

(3-0.2)/(2*6^5+1)=1.8003e-04

so the delta x is REALLY small.

So the last entry corresponds to the error of simp38.m at this delta x.

6^(-18.4166)= 4.6677e-15

so the error there corresponds to the epsilon of the answer. You can't do better than that.

The trapezoid rule error is much bigger:

6^(-10.5096) = 6.6365e-09