echo on
%%%%% ch08_exa.m %%%%%
% Run script and study the effects of each line.

echo off
format compact;
format short;
disp(sprintf('**** Examples from Chapter 08 ****'));


disp(sprintf('**** Example 8.1 Integral of Lagrange basis functions *'));
echo on
%%% Calculate weights of quadrature rule by integrals of interpolating 
%%% Lagrange polynomial.
% Data points
xi=[0 0.5 1]'; % If xi are predifined as eqidistant this results in Newton-Cotes method
x=sym('x');
% Let us define Lagrange polynomial of order 2, that interpolates 3 data points. 
% Basis functions
l1 = simplify(((x-xi(2))*(x-xi(3)))/((xi(1)-xi(2))*(xi(1)-xi(3))))
l2 = simplify(((x-xi(1))*(x-xi(3)))/((xi(2)-xi(1))*(xi(2)-xi(3))))
l3 = simplify(((x-xi(1))*(x-xi(2)))/((xi(3)-xi(2))*(xi(3)-xi(1))))
% Interpolating polynomial is not needed explicitely. Only weights are searched.
% Integral of interpolating polynomial, gives the aproximate integral of
% the underlying function as a weigted sum of function evaluated at knot values.
a=sym('a');
b=sym('b');
w1 = int(l1, a, b)
w2 = int(l2, a, b)
w3 = int(l3, a, b)
% evaluated for a=xi(1) and b=xi(3) we obtain Simpson's rule.
a=xi(1);
b=xi(3);
ew1=eval(w1)
ew2=eval(w2)
ew3=eval(w3)

disp(sprintf('**** Example 8.1 Undetermined coefficients *'));
echo on
%%% Calculate weights of quadrature rule by undetermined coefficients.
a=sym('a');
b=sym('b');
im1 = int(1, a, b) % exact for the monomial of degree 1
im2 = int(x, a, b) % exact for the monomial of degree 2
im3 = int(x^2, a, b) % exact for the monomial of degree 3
a=xi(1);
b=xi(3);
eim1=eval(im1); % evaluate for a and b
eim2=eval(im2);
eim3=eval(im3);
% moment equations have to be solved for w
%w1*1 + w2*1 + w3*3 = im1
%w1*xi(1) + w2*xi(2) + w3*xi(3) = im2
%w1*xi(1)^2 + w2*xi(2)^2 + w3*xi(3)^2 = im3
V = [1 1 1; xi(1) xi(2) xi(3); xi(1)^2 xi(2)^2 xi(3)^2;]
w=V\[eim1 eim2 eim3]' % same weights as in the previous example.


disp(sprintf('**** Example 8.2 Newton-Cotes Quadrature *'));
echo on
%%% Nodes of Newton-Cotes quadrature are eqidistant.
%%% Compute the approximate value of the integral of exp(-x^2) from a to b.
%MatLab quadrature for the function exp(-x^2) from file ch08_f.m
a=xi(1)
b=xi(3)
quad ('ch08_f', a, b)

% Midpoint
M =(b-a)*exp(-((a+b)/2)^2)
% Trapezoid
T = (b-a)/2*(exp(-a^2)+ exp(-b^2))
% Error estimate:
E=(T-M)/3
%I=M+E and I=T-2*E
% Simpson's
S = (b-a)/6*(exp(-a^2)+ 4*exp(-((a+b)/2)^2) + exp(-b^2))
% and by definition with the calculated weigths 
N3 = w(1)* exp(-xi(1)^2)+ w(2)* exp(-xi(2)^2)+ w(3)* exp(-xi(3)^2)
% and finally, the alternative formula:
S1=2/3*M + 1/3*T

disp(sprintf('**** Example 8.2.1 Clenshaw-Curtis Quadrature *'));
echo on
%%% Nodes of Clenshaw-Curtis quadrature are Chebyshev points (zeros or extrema)
%%% and defined on the interval [-1 1].
%%% Calculate weights of Clenshaw-Curtis quadrature by integrals of interpolating 
%%% Lagrange polynomial.
% Data points
ti=[cos(pi) cos(pi/2) cos(0)]'; % Chebyshev points - extrema including endpoints
t=sym('t');
% Let us define Lagrange polynomial of order 2, that interpolates 3 data points. 
% Basis functions
l1 = simplify(((t-ti(2))*(t-ti(3)))/((ti(1)-ti(2))*(ti(1)-ti(3))))
l2 = simplify(((t-ti(1))*(t-ti(3)))/((ti(2)-ti(1))*(ti(2)-ti(3))))
l3 = simplify(((t-ti(1))*(t-ti(2)))/((ti(3)-ti(2))*(ti(3)-ti(1))))
% Interpolating polynomial is not needed explicitely. Only weights are searched.
% Integral of interpolating polynomial, gives the aproximate integral of
% the underlying function as a weigted sum of function evaluated at knot values.
an=sym('an');
bn=sym('bn');
w1 = int(l1, an, bn)
w2 = int(l2, an, bn)
w3 = int(l3, an, bn)
% evaluated for an=ti(1) and bn=ti(3) we obtain Clenshaw-Curtis, that is
% identical to the Simpson's rule (because we have only 3 nodes), for interval [-1 1].
an=ti(1);
bn=ti(3);
ew1=eval(w1)
ew2=eval(w2)
ew3=eval(w3) 
% Compute the approximate integral of exp(-t^2) from -1 to 1.
CC3 = ew1* exp(-ti(1)^2)+ ew2* exp(-ti(2)^2)+ ew3* exp(-ti(3)^2)

%%% Compute again the approximate value of the integral of exp(-t^2) from 0 to 1.
% The transformation of the interval [an bn]=[-1 1] to the interval [a b]=[0 1] can be 
% performed by a transformation 
tti=((b-a)*ti + a*bn - b*an)/(bn-an)
CC3 = ((b-a)/(bn-an))* (ew1* exp(-tti(1)^2)+ ew2* exp(-tti(2)^2)+ ew3* exp(-tti(3)^2))


disp(sprintf('**** Example 8.4 Gaussian Quadrature *'));
echo on
%%% Neither wights nor nodes of Gaussian quadrature rule are not defined in advance
%%% Integration interval is given i.e. [-1 1]=[an bn]
% We have 2n moment equations that have to be solved for w and xi
an=sym('an')
bn=sym('bn')
im1 = int(1, an, bn) % exact for the monomial of degree 1
im2 = int(x, an, bn) % exact for the monomial of degree 2
im3 = int(x^2, an, bn) % exact for the monomial of degree 3
im4 = int(x^3, an, bn) % exact for the monomial of degree 4
an=-1;
bn=1;
eim1=eval(im1) % evaluate for a and b
eim2=eval(im2)
eim3=eval(im3)
eim4=eval(im4)
%w1*1 + w2*1 = eim1
%w1*xi(1) + w2*xi(2) = eim2
%w1*xi(1)^2 + w2*xi(2)^2 = eim3
%w1*xi(1)^3 + w2*xi(2)^3 = eim4
% Nonlinear system should be solved by an iterative method.
x0=nonlin;
w1=x0(1)
w2=x0(2)
xg(1)=x0(3)
xg(2)=x0(4)
G2 = w1* exp(-xg(1)^2)+ w2* exp(-xg(2)^2)

%%% Compute again the approximate value of the integral of exp(-t^2) from 0 to 1.
% The transformation of the interval [an bn]=[-1 1] to the interval [a b]=[0 1] can be 
% performed by a transformation 
tg=((b-a)*xg + a*bn - b*an)/(bn-an)
G2 = ((b-a)/(bn-an)) * (w1* exp(-tg(1)^2)+ w2* exp(-tg(2)^2))