%% Euler Maclaurin and Trapezoidal Formulae
% Consider the sum S = \left(\frac 1{10}\right)^2 +  \left(\frac 1{11}\right)^2+ \left(\frac 1{12}\right)^2+\cdots+ \left(\frac 1{100}\right)^2. We would like an approximation without actually adding all of these terms. The exact value is

clear
format rat

f = @(x) x.^(-2);
S = sum(f(10:100))
format long
disp(['This is S = ', num2str(S)])

%
% Recall the Trapezoidal Formula from Calculus II: \int_a^b f(x)\, dx \approx \frac{\Delta x}{2}\left[f(x_0)+2f(x_1)+\cdots + 2 f(x_{n-1} )+ f(x_n)\right].
% Letting \Delta x=1and omitting half of the first and last terms in S, we see that S \approx \int_{10}^{100} \frac{dx}{x^2} +\frac 12 \left(\frac 1{10}\right)^2+\frac 12 \left(\frac 1{100}\right)^2 = \frac 1{10}-\frac 1{100} +\frac 12 \left(\frac 1{10}\right)^2+\frac 12 \left(\frac 1{100}\right)^2.
%

Sapprox =1/10-1/100+.5*1/10^2+.5*1/100^2

Error = abs(S-Sapprox)

%
% We are interested in finding the error in approximating an integral using the Trapezoidal Rule.
%

%f = @(x) cos(x).*sin(x);

a = 10;
b = 100;
n  = 90; 
h = (b-a)/n;
t = a:h:b;
t2 = a:h/10:b;
S2 = 0;
plot(t2,f(t2),'b','LineWidth',2);
hold on;
for i=1:n-1
    S2 = S2 + f(a+i*h);
    area([a+i*h,a+i*h,a+(i+1)*h,a+(i+1)*h],[0,f(a+i*h),f(a+(i+1)*h),0]);
end
hold off;
result = h/2*(f(b)+f(a)+2*S2)  % Integral approx.
result+.5*(f(a)+f(b))         % Sum approx.

%
% The error estimation comes from the Euler-Maclaurin Summation Formula:
% h\sum_{k=0}^n f(x_k)=\int_a^b f(x)\, dx +\frac{h} {2}\left[f(x_0)+ f(x_n)\right] + \sum_{\ell=1} \frac{B_{2\ell}}{(2\ell)!}h^{2\ell} \left(f^{(2\ell-1)}(b)-f^{(2\ell-1)}(a)\right).
% Here the B_n's are the Bernoulli numbers defined through the series \frac{x}{e^x-1}=\displaystyle\sum_{n=0}^\infty \frac{B_n}{n!}x^n.

format rat
for k=0:6
    B(k+1)=bernoulli(k);
end
B

%
% As a first bound on the error, we have \left| \int_a^b f(x)\, dx -T_f(h) \right| \approx  \frac{h^2}{12}\left|f'(b)-f'(a)\right|. We can verify this bound:

format long
bound=bernoulli(2)/factorial(2)*h^2*(b^(-3)-a^(-3))*factorial(2)
integral(f,10,100)-result

%
%% Trapezoidal Rule
% We approximate the integral of the function f(x)=e^{7x} over the interval [0,2].

f = @(x) exp(-7*x);
a = 0;  b = 2;

%
% In lieu of the exact value, we will use the built-in integral function to find an accurate result.

I = integral(f,a,b,'abstol',1e-14,'reltol',1e-14);
fprintf('Integral = %.15f\n',I)

%
% Here is the error at n=40 and the error bound from the Euler-Maclaurin Formula.

n=40;
h=(b-a)/n;
T = trapezoid(f,a,b,n);
err = I - T
bound=bernoulli(2)/factorial(2)*h^2*(-7)*(f(b)-f(a))

%
%% Extrapolation
% Double n and only use every other value. Then repeat two more times.

n = 2*n;  h = h/2;  t = h*(0:n)';
T(2) = T(1)/2 + h*sum( f(t(2:2:n)) );
err_2nd = I - T

%

n = 2*n;  h = h/2;  t = h*(0:n)';
T(3) = T(2)/2 + h*sum( f(t(2:2:n)) );
err_2nd = I - T

%

n = 2*n;  h = h/2;  t = h*(0:n)';
T(3) = T(2)/2 + h*sum( f(t(2:2:n)) );
err_2nd = I - T

%
% Use these results to construct Simpson approximation.

S = (4*T(2:3) - T(1:2)) / 3;
err_4th = I - S

%
% Repeat to obtain Romberg approximation.

R = (16*S(2) - S(1)) / 15;
err_6th = I - R

%
%% Only functions are below ....
%

function [T,t,y] = trapezoid(f,a,b,n)
%TRAPEZOID   Trapezoid formula for numerical integration.
% Input:
%   f     integrand (function)
%   a,b   interval of integration (scalars)
%   n     number of interval divisions
% Output:
%   T     approximation to the integral of f over (a,b)
%   t     vector of nodes used
%   y     vector of function values at nodes

h = (b-a)/n;
t = a + h*(0:n)';
y = f(t);
T = h * ( sum(y(2:n)) + 0.5*(y(1) + y(n+1)) ); 
end