%% Integration rules
% We want to integrate some function f\left(x\right)
% \int_a^b f\left(x\right)\;d\;x\;
% We can take any range x\in \left\lbrack a,\;b\right\rbrack and map it to the range \xi \in \left\lbrack 0,1\right\rbrack using iso-parametric mapping:
% x\left(\xi \right):=\left(1-\xi \right)a+\xi b and its derivative (or map) \frac{d\;x}{d\;\xi }=b-a then we have
% \int_a^b f\left(x\right)\;d\;x\;=\int_0^1 f\left(x\left(\xi \right)\right)\;\frac{d\;x}{d\;\xi }\;d\xi =\left(b-a\right)\int_0^1 f\left(x\left(\xi \right)\right)\;d\;\xi
% Let's say we want to integrate some polynomial of order n
% \int_a^b x^n \;d\;x\;=\left(b-a\right)\int_0^1 {x\left(\xi \right)}^n \;d\;\xi =\left(b-a\right)\int_0^1 {\left(\left(1-\xi \right)a+\xi b\right)}^n d\xi =-\frac{a^{n+1} -b^{n+1} }{n+1}
% Ok, so this can be used to integrate e.g.,
% \int_2^5 2x-4x^2 \;d\;x=-2\frac{a^{1+1} -b^{1+1} }{1+1}--4\frac{a^{2+1} -b^{2+1} }{2+1}=-135
% But in general it is awkward to use this rule to integrate e.g.,
% \int_{x_1 }^{x_2 } \left(1-x\right)\left(1-x\right)\;d\;x or \int_{x_1 }^{x_2 } {\underline{\varphi} }^T \underline{\varphi} \;d\;x without additional preparation which makes it cumbersome to implement in a computer program and run efficiently.
% The same goes for the integration rulse such as
% \int_{-1}^1 x^n \;d\;x=\left\lbrace \begin{array}{ll}
% \frac{2}{n+1} & \textrm{if}\;n\;\textrm{is}\;\textrm{even}\\
% 0 & \textrm{if}\;n\;\textrm{is}\;\textrm{odd}
% \end{array}\right.
% and
% \int_0^1 x^n \;d\;x=\frac{1}{n+1}
%% Numerical integration
% Let's consider a polynomial function f\left(x\right)=1+3x^2 -x^3 for x\in \left\lbrack 0,3\right\rbrack

clear

f = @(x) 1+3*x.^2-x.^3;
a = 0; b = 3;
figure
fplot(f,[a,b])
axis([a,b,0,5.1]); box off;
legend('show','Location','Northwest');

% We can integrate this analytically

I_e = int(f,sym("x"),a,b)
I_e = double(I_e)

% We can also use the numerical integration using the build in variable precision arithmatic integration

vpaintegral(f,sym("x"),a,b)

% This uses a numerical method for integration, is powerful and easy to use, drawback are a default high number of function evaluations.
% Now, let's go back to basics and take a look at the Riemann sum.
% \int_a^b f\left(x\right)\approx \frac{b-a}{N}\sum_{i=1}^N f\left(x_i \right)

N = 5;
xnod = linspace(a,b,N);
y = double(f(xnod));
figure
fplot(f,[a,b])
axis([a,b,0,5.1]); box off;
legend('show','Location','Northwest');hold on
hb = bar(xnod,y,1,'LineStyle','-',DisplayName="Riemann Sum");
dx = (b-a)/N;
I_R = sum(y*dx);
ht = title(['$I\approx$ ', sprintf('%0.2f', I_R),', $I = $', sprintf(' %0.2f', I_e),', $N=$ ',sprintf(' %d', N)]);
ht.Interpreter = "latex";
% exportgraphics(gca,"Riemann.gif","Append",false)

%
%
% Let's take a look at the error we're making numerically
% \epsilon ={\left|I\right.}_{\textrm{Riemann}} -I_{\textrm{exact}} \left|\right.

N = 200;
xnod = linspace(a,b,N);
dx = (b-a)/N;
I_R = sum(f(xnod)*dx);
err_R = abs(I_R-I_e)

% Note that the absolute integration error is still significantly large after 200 function evaluations!
%% Trapezoidal rule
% We can do a better numerical approximation
% First we approximate f\left(x\right)\;linearly: f\left(x\right)\approx \varphi_1 \left(x\right)f\left(a\right)+\varphi_2 \left(x\right)f\left(b\right)=\frac{b-x}{b-a}f\left(a\right)+\frac{x-a}{b-a}f\left(b\right) then we have
% I\approx \frac{f\left(a\right)}{b-a}\int_a^b b-x\;d\;x+\frac{f\left(b\right)}{b-a}\int_a^b x-a\;d\;x=\frac{f\left(a\right)}{b-a}{\left\lbrack b\;x-\frac{1}{2}x^2 \right\rbrack }_a^b +\frac{f\left(b\right)}{b-a}{\left\lbrack \frac{1}{2}x^2 -a\;x\right\rbrack }_a^b
% =\frac{f\left(a\right)}{b-a}\left(b^2 -\frac{1}{2}b^2 -a\;b-\frac{1}{2}a^2 \right)+\frac{f\left(b\right)}{b-a}\left(\frac{1}{2}b^2 -a\;b-\frac{1}{2}a^2 -a^2 \right)
% =\frac{f\left(a\right)}{b-a}\frac{1}{2}{\left(b-a\right)}^2 +\frac{f\left(b\right)}{b-a}\frac{1}{2}{\left(b-a\right)}^2
% =\frac{1}{2}\left(b-a\right)\left(f\left(a\right)+f\left(b\right)\right)
% This is called the trapezoidal rule. And for a chained approximation with N+1 evenly spaced points the approximation becomes:
% \int_a^b f\left(x\right)d\;x\approx \frac{b-a}{2N}\sum_{i=1}^N \left(f\left(x_i \right)+f\left(x_{i+1} \right)\right)
% For example for the function above using N=9  we have:

N = 9;

xnod = linspace(a,b,N+1)';

I_R = double(sum(f(xnod)*(b-a)/(N+1)))
err_R = abs(I_R-I_e)

xnod = linspace(a,b,N+1)';
x1 = xnod(1:end-1);
x2 = xnod(2:end);
I_T = sum(f(x1)+f(x2))*(b-a)/(2*N)
double(I_T)
I_T2 = trapz(xnod,f(xnod)) % Matlab built in trapezoidal integral
double(I_T2)
err_T = double(abs(I_T-I_e))

err_R/err_T

% 9 times better approximation!
%


figure
axis([a,b,0,5.1]); box off;
legend('show','Location','Northwest');
hold on
nodes = [(1:N)', (2:N+1)'];
xm = mean(xnod(nodes),2);
if N < 15
    for iel = 1:N
        x1 = xnod(nodes(iel,1));
        x2 = xnod(nodes(iel,2));
        h = (b-a)/N;
        y1 = f(x1);
        y2 = f(x2);
        area([x1,x2],[y1, y2],'FaceColor','c','EdgeColor','k','HandleVisibility','off');
        plot([x1,x1],[0,y1],'k:','HandleVisibility','off')
        plot([x2,x2],[0,y2],'k:','HandleVisibility','off')
    end
    text(xm,xm*0+0.2,cellstr(num2str([1:N]')))
else
    area(xnod,f(xnod),'FaceColor','c','HandleVisibility','off');
end
fplot(f,[a,b],'r','LineWidth',2)
plot(xnod,f(xnod),'ok','HandleVisibility','off')

%
%
%% Exact numerical integration
% The trapezoidal rule integrates linear polynomials exactly:
% \int_a^b f\left(x\right)d\;x\approx \frac{b-a}{2}\left(f\left(a\right)+f\left(b\right)\right)

f = @(x)1+3*x;
a = 0; b = 3;
I1 = (f(a)+f(b))*(b-a)/2
I_exact = vpa(int(sym(f),a,b),5)

figure; hold on
area([a,b],[f(a), f(b)],'FaceColor','c');
fplot(f,[a,b])
plot([a,b],[f(a), f(b)],'b')

% But there is a cheaper way...
% We can evaluate the midpoint \frac{\left(a+b\right)}{2} instead:

xm = (a+b)/2;
I1 = (b-a)*f(xm)

% One function evaluation instead of 2!

figure; hold on

area([a,b],[f(xm), f(xm)],'FaceColor','c');
plot([a,b],[f(a), f(b)],'b-')
plot(xm,f(xm),'ko','MarkerSize',14,'MarkerFaceColor','k')

% This is exact because the error above the graph cancelles the error below the graph!
%
%% Gauss quadrature rule
% One rule to rule them all.
%
% This error cancelling property can be generelized so that n integration points integrate polynomials of degree 2n-1 exactly!
% Gauss quadrature states for x\in \left\lbrack 0,1\right\rbrack
% \int_0^1 p\left(x\right)\;d\;x\approx \sum_i^{\;} p\left(x_i \right)w_i
% \begin{array}{cccc}
% \textrm{Gauss}\;\textrm{order} & \textrm{Polynomial}\;\textrm{order} & \textrm{Weights} & \textrm{Points}\\
% 1 & 1 & 1 & \frac{1}{2}\\
% 2 & 3 & \frac{1}{2} & \frac{1}{2}-\frac{\sqrt{3}}{6}\\
% \; & \; & \frac{1}{2} & \frac{1}{2}+\frac{\sqrt{3}}{6}\\
% 3 & 5 & \frac{5}{18} & \frac{1}{2}-\frac{\sqrt{3}\sqrt{5}}{10}\\
% \; & \; & \frac{4}{9} & \frac{1}{2}\\
% \; & \; & \frac{5}{18} & \frac{1}{2}+\frac{\sqrt{3}\sqrt{5}}{10}
% \end{array}
% Example:

syms x
p = @(x) x;
I_e  = int(sym(p),x,0,1)

gp = 1/2; gw = 1;
I_G = p(gp)*gw

% Ok, not a big deal, after all, even the Riemann sum is exact for linear polynomials.
% How about a third order polynomial using only two function evaluations?
%
% \int_0^1 1+3x-{3x}^3 \;d\;x=\sum_i^2 f\left(x_i \right)w_i
% where x_1 =\frac{1}{2}-\frac{\sqrt{3}}{6} and x_2 =\frac{1}{2}+\frac{\sqrt{3}}{6}  and w_{1,2} =\frac{1}{2}

p = @(x)1+3*x.^2-3*x.^3;

gp = [1/2-sym(sqrt(3))/6, 1/2+sym(sqrt(3))/6]
gw = sym([1/2, 1/2])

I_Gauss = p(gp(1))*gw(1) + p(gp(2))*gw(2)
vpa(I_Gauss)
simplify(I_Gauss)

xnod = linspace(0,1,50);
I_trapz = vpa(trapz(xnod,p(xnod)))

I_exact = int(sym(p),x,0,1)

figure; hold on
area([0,1/2],[1,1]*p(gp(1)),'FaceColor','c');
area([1/2,1],[1,1]*p(gp(2)),'FaceColor','c');
fplot(p,[0,1], 'r', 'LineWidth',2)
plot(gp,p(gp),'k*')
plot([0.5, 0.5], [0, p(gp(2))], '-k')

% The figure visualizes the Gauss integration rule in action. The cyan and white areas between the curves are cancelling out!
% Higher order Gauss integration will not make the integral any more accurate, the computed value will be the same. We will just be wasting computational power.
%
%% Non trivial integrals
% Note that Gauss is not going to exactly integrade other functions in general.
%% Example - A tricky integral
% Let's say we need to integrate the function
% f\left(x\right)=\frac{1}{\sqrt{\sin \left(x\right)}}\;\;\textrm{for}\;x\in \left\lbrack 0,1\right\rbrack
% The resulting integral
% \int_0^1 \frac{1}{\sqrt{\sin \left(x\right)}}\;\;d\;x
% is non-elementary.

clear
syms x
f(x) = 1/sqrt(sin(x))
figure
fplot(f,[0,1],'k','LineWidth',2)
axis([0,1,0,10])
xnod = linspace(0.01,1,1000);
ynod = double(f(xnod));
hold on
area(xnod,ynod,'FaceColor','c','EdgeColor','none')
int(f)
int(f,x,0,1)

% We can of course use numerical integration to solve this

I_e = vpaintegral(f,x,0,1)

% We can apply Gauss integration using the tabulated values from above:

f = matlabFunction(f);

[gp, gw] = GaussTable(3);
I_G = f(gp)*gw
error_Gauss = I_e-I_G

% We can see that 3 Gauss points is not enough to sufficiently integrate this function. We need more points!
% Let's go down the rabit hole, there exists quadrature tables containing Gauss points for the typical range \left\lbrack -1,1\right\rbrack but there are also several algorithms for computing Gauss points for any Gauss order. Well use the well known Golub-Welsh algorithm taken from L. N. Trefethen, Spectral Methods in Matlab, SIAM, 2000 p.129

f = @(x)1./sqrt(sin(x))
a = 0; b = 1;
N = 17;
[gp, gw] = gauss(N, a, b);
I_G = vpa((b-a)*f(gp)*gw)

error_Gauss = abs(I_e-I_G)

xnod = linspace(a+0.01,b,N)';
I_T = trapz(xnod,f(xnod))
error_Trapz = abs(I_e-I_T)

% We can plot the error versus the number of evaluations


a = 0; b = 1;
err_G = []; err_Trapz = [];
NP = [2,4,8,16,32,64,128,256];
for N = NP
    [gp, gw] = gauss(N, a, b);
    I_G = vpa((b-a)*f(gp)*gw);
    xnod = linspace(a+0.003,b,N)';
    I_T = trapz(xnod,f(xnod));

    error_Gauss = abs(I_e-I_G);
    error_Trapz = abs(I_e-I_T);

    err_G(end+1) = error_Gauss;
    err_Trapz(end+1) = error_Trapz;
end

figure;
plot(NP,err_G,'-o','DisplayName','Gauss error'); hold on
plot(NP,err_Trapz,'-x','DisplayName','Trapz error')
gca().set(XScale='log',YScale='log')
legend show
xlabel('N evaluation points')
ylabel('Absolute error')
grid on
xticks(NP)
xticklabels({'2','4','8','16','32','64','128','256'})

% We can see that both methods converge to the numerical integral computed by vpaintegral. The trapezoidal integral here is volotile, its lower bound needs to be adjusted to avoid the obious error division by zero. The computation thus depends on the choice of the new lower bound, and for a careful choice of lower bound we can get the trapezoidal method to outperform Gauss for some intermediate number of evaluations. The problem is that we cannot know this value ahead of time and thus the method is not very suitable for this situation.
%
%% Example - A sinusoidal
% u\left(x\right)=2\;x\;\sin \left(2\;\pi \;x\right)+3\;\;\textrm{for}\;x\in \left\lbrack -0\ldotp 3,0\ldotp 9\right\rbrack
% If we want to apply Gaussian quadrature here, we need to make a variable change to integrate over ranges other that \left\lbrack 0,1\right\rbrack.
% We can take any range x\in \left\lbrack a,\;b\right\rbrack and map it to the range \xi \in \left\lbrack 0,1\right\rbrack using iso-parametric mapping:
% x\left(\xi \right):=\left(1-\xi \right)a+\xi b and its derivative (or map) \frac{d\;x}{d\;\xi }=b-a then we have
% \int_a^b f\left(x\right)\;d\;x\;=\int_0^1 f\left(x\left(\xi \right)\right)\;\frac{d\;x}{d\;\xi }\;d\xi =\left(b-a\right)\int_0^1 f\left(x\left(\xi \right)\right)\;d\;\xi
% and using Gauss quadrature approximation we have
% \left(b-a\right)\int_0^1 f\left(x\left(\xi \right)\right)\;d\;\xi \approx \left(b-a\right)\sum_i^N f\left(x_i \right)w_i

clear
syms x
u(x) = 2*x*sin(2*pi*x)+3
a = -0.3; b = 0.9;
figure
fplot(u,[a,b],'k','LineWidth',2)
% axis([a,b,0,10])
xnod = linspace(a,b,1000);
ynod = double(u(xnod));
hold on
area(xnod,ynod,'FaceColor','c','EdgeColor','none')

I_e = int(u,x,a,b);
I_e = vpa(I_e)

% And now let's approximate ut using Gauss quadrature

f = matlabFunction(u);
N = 8;
[gp, gw] = gauss(N, a, b);
I_G = vpa((b-a)*f(gp)*gw)

error_Gauss = I_e-I_G

xnod = linspace(a,b,N)';
I_T = vpa(trapz(xnod,f(xnod)))
error_Trapz = I_e-I_T

% Error analysis

error.N = [2,4,8,16,32,64,128,256].'
error.Gauss = zeros(length(error.N),1);
error.Trapz = zeros(length(error.N),1);

for i = 1:length(error.N)
    N = error.N(i);
    xnod = linspace(a,b,N)';
    I_Trapz = trapz(xnod,f(xnod)); % Matlab built in trapezoidal integral
    
    [gp, gw] = gauss(N, a, b);
    I_Gauss = (b-a)*f(gp)*gw;
    
    error_Gauss = abs(I_e-I_Gauss);
    error_Trapz = abs(I_e-I_Trapz);
    
    error.Gauss(i) = error_Gauss;
    error.Trapz(i) = error_Trapz;
    
end
figure; hold on
plot(error.N, error.Gauss,'-o')
plot(error.N, error.Trapz,'-o')
grid on
set(gca,'xscale','log','yscale','log')
legend('Gauss','Trapz')
xlabel('N evaluation points')
ylabel('Absolute error')
xticks(error.N)
xticklabels({'2','4','8','16','32','64','128','256'})

% Here we see that Gauss reaches the same precision as vpaintegral after 16 function evaluations.
%
%
%

function [gp, gw] = GaussTable(n)

    switch n
        case 1
            gp = 0.5;
            gw = 1;
        case 2
            gp = [0.5-sqrt(3)/6
                  0.5+sqrt(3)/6]';
            gw = [0.5
                  0.5];
        case 3
            gp = [0.5-sqrt(3)*sqrt(5)/10
                  0.5
                  0.5+sqrt(3)*sqrt(5)/10]';
            gw = [5/18
                  4/9
                  5/18];
    end

end

function [x, w, A] = gauss(n, a, b)
    
    %------------------------------------------------------------------------------
    % gauss.m
    %------------------------------------------------------------------------------
    %
    % Purpose:
    %
    % Generates abscissas and weigths on I = [ a, b ] (for Gaussian quadrature).
    %
    %
    % Syntax:
    %
    % [x, w, A] = gauss(n, a, b);
    %
    %
    % Input:
    %
    % n    integer    Number of quadrature points.
    % a    real       Left endpoint of interval.
    % b    real       Right endpoint of interval.
    %
    %
    % Output:
    %
    % x    real       Quadrature points.
    % w    real       Weigths.
    % A    real       Area of domain.
    %------------------------------------------------------------------------------
    
    
    % 3-term recurrence coefficients:
    n = 1:(n - 1);
    beta = 1 ./ sqrt(4 - 1 ./ (n .* n));
    
    % Jacobi matrix:
    J = diag(beta, 1) + diag(beta, -1); 
    
    
    % Eigenvalue decomposition:
    
    %
    % e-values are used for abscissas, whereas e-vectors determine weights.
    %
    
    [V, D] = eig(J);
    x = diag(D);
    
    
    % Size of domain:
    A = b - a;
    
    
    % Change of interval:
    
    %
    % Initally we have I0 = [ -1, 1 ]; now one gets I = [ a, b ] instead.
    %
    % The abscissas are Legendre points.
    %
    
    if ~(a == -1 && b == 1)
      x = 0.5 * A * x + 0.5 * (b + a);
    end
    
    
    % Weigths:
    w = V(1, :) .* V(1, :);
    w = w';
    x=x';
end