% CLASS17NOV   Demonstrate various Newton-Cotes rules using COMPARERULES
% and then extract, by Romberg magic, great results from the Trapezoid rule.

% make sure comparerules.m is downloaded and in the current directory!

% run COMPARERULES, so we can extract Trapezoid results
figure(1)
comparerules
printf('figure 1 just compares Trapezoid, Simpsons, Midpoint\n')

% look at Trapezoid results and exact answer
format long g
exact = 1-exp(-1)
trapN16 = results(4,1)
printf('(error from N=16 Trapezoid) = %.5e\n',abs(trapN16-exact))
trapN128 = results(7,1)
printf('(error from N=128 Trapezoid) = %.5e\n',abs(trapN128-exact))

% plot h versus Trapezoid results
figure(2)
plot(hlist,results(:,1),'*','markersize',15)
xlabel h,  ylabel('Trapezoid result using h')
axis tight,  grid on
printf('figure 2 shows Trapezoid results as function of h\n')

% use polynomial interpolation to extrapolate to h=0  !!
x = hlist;
y = results(:,1)';
P = polyfit(x,y,6);
romberg1 = polyval(P,0)
printf('(error from 7 pt Romberg-like extrapolation) = %.5e\n',abs(romberg1-exact))

% subset of first 4 Trapezoid results, plus extrapolation
x = hlist(1:4);
y = results(1:4,1)';
P = polyfit(x,y,3);
romberg2 = polyval(P,0)
printf('(error from subset (4 pt) Romberg-like extrapolation) = %.5e\n',abs(romberg2-exact))

% actual Romberg result: extrapolate with h^2 as independent var
P = polyfit(x.^2,y,3);
romberg3 = polyval(P,0)
printf('(error from subset (4 pt) true Romberg extrapolation) = %.5e\n',abs(romberg3-exact))

