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

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

disp(sprintf('**** Example 8.1 Integral of Lagrange basis functions *'));

%%%Calculate weights of quadrate rule by integrals of interpolating 
%%%Lagrange polynomial. Test the result on f=exp(-x^2) on interval [0,1] with
%%%data points
xi=[0 0.5 1]';  %If xi are equidistant we get Newton-Cotes rule.
a=xi(1);    %Interval boundaries
b=xi(3);
%Define Lagrange polynomial of order 2, that interpolates 3 data points. 
%(Simpson's rule). It basis functions:
%l1 = ((x-xi(2))*(x-xi(3)))/((xi(1)-xi(2))*(xi(1)-xi(3)))
%l2 = ((x-xi(1))*(x-xi(3)))/((xi(2)-xi(1))*(xi(2)-xi(3)))
%l3 = ((x-xi(1))*(x-xi(2)))/((xi(3)-xi(2))*(xi(3)-xi(1)))
%Interpolating polynomial is not needed explicitly. Only weights are searched
%that are equal to integrals of basis functions. 
w1 = (b-a)/6    %integral(l1, a, b)
w2 = 4*(b-a)/6    %integral(l2, a, b)
w3 = (b-a)/6    %integral(l3, a, b)
%Integral of interpolating polynomial gives the approximate integral of
%the underlying function as a weighted sum of function evaluated at knot values.
%Evaluate function at nodes.
fa=exp(-xi(1)^2)
fb=exp(-xi(3)^2)
fm=exp(-xi(2)^2)
%And the approximation of the integral:
I=w1*fa + w2*fm + w3*fb
%Exact value of the integral is 0.746824


disp(sprintf('**** Example 8.1 Undetermined coefficients *'));

%%%Calculate weights of quadrature rule, for data from Exa. 8.1 using undetermined coefficients.
%Moment equations have to be solved for w:
%w1*1 + w2*1 + w3*1 = b-a    %integral(1, a, b)
%w1*xi(1) + w2*xi(2) + w3*xi(3) = (b^2-a^2)/2    %integral(x, a, b)
%w1*xi(1)^2 + w2*xi(2)^2 + w3*xi(3)^2 = (b^3-a^3)/3    %integral(x, a, b)
V = [1 1 1; xi(1) xi(2) xi(3); xi(1)^2 xi(2)^2 xi(3)^2;]
int = [b-a, (b^2-a^2)/2, (b^3-a^3)/3]';
w=V\int 
%same weights as in Exa. 8.1, so approximate integral has also the same value.

disp(sprintf('**** Example 8.2 Newton-Cotes Quadrature *'));

%%%Nodes of Newton-Cotes quadrature are equidistant.
%%%Compute the approximate value of the integral of exp(-x^2) from a to b.
%MatLab quadrature for the function exp(-x^2)can be calculated by quad() function.
a=xi(1)
b=xi(3)
quad ('exp(-x.^2)', a, b)

%Midpoint rule(a single node)
M =(b-a)*exp(-((a+b)/2)^2)
%Trapezoid rule(two nodes on interval border)
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 rule (tree nodes: border and centre of interval)
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 *'));

%%%Nodes of Clenshaw-Curtis quadrature are Chebyshev points (zeros or extrema)
%%%and defined on the interval [-1 1].
%If we take extrema including endpoints of Chebyshev polynomial T2
%ti=[cos(pi) cos(pi/2) cos(0)]'; we get the rule identical to Simpson's rule.
%Try wit zeros of Chebyshev polynomial T3. We have 3 zeros:
ti=[cos(5*pi/6) cos(3*pi/6) cos(pi/6)]'
an=-1;
bn=1;
%Using method of undetermined coefficients for polynomial of order 2, that interpolates
%3 data points we get weights in the same way as in Exa. 8.2
V = [1 1 1; ti(1) ti(2) ti(3); ti(1)^2 ti(2)^2 ti(3)^2;]
int = [bn-an, (bn^2-an^2)/2, (bn^3-an^3)/3]';
ew=V\int 
%We obtain Clenshaw-Curtis rule.
%Compute the approximate integral of exp(-t^2) from -1 to 1.
CC3 = ew(1)* exp(-ti(1)^2)+ ew(2)* exp(-ti(2)^2)+ ew(3)* exp(-ti(3)^2)
%Matlab value for this integral:
quad ('exp(-x.^2)', an, bn)

%%%Compute again the approximate value of the integral of exp(-t^2) on [0,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)
%Instead of xi use tti:
CC3 = ((b-a)/(bn-an))* (ew(1)* exp(-tti(1)^2)+ ew(2)* exp(-tti(2)^2)+ ew(3)* exp(-tti(3)^2))


disp(sprintf('**** Example 8.4 Gaussian Quadrature *'));

%%%Neither weights nor nodes of Gaussian quadrature rule are 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 example for n=2 (2 nodes gxi), using integrals of Lagrange basis functions:
%w1*1 + w2*1 = gint1
%w1*xi(1) + w2*xi(2) = gint2
%w1*xi(1)^2 + w2*xi(2)^2 = gint3
%w1*xi(1)^3 + w2*xi(2)^3 = gint4
gint = [bn-an, (bn^2-an^2)/2, (bn^3-an^3)/3, (bn^4-an^4)/4]';
%V = [1 1; gxi(1) gxi(2); gxi(1)^2 gxi(2)^2; gxi(1)^3 gxi(2)^3;]
%Nonlinear system should be solved by an iterative method for xw=[gx1; gx2; gw1; gw2].
%Needs one evaluation of multidimensional function
%F(gx,gw)=[gw1+gw2-gint(1); gw1*gx1 + gw2*gx2; gw1*gx1^2 + gw2*gx2^2-gint(3);gw1*gx1^3 + gw2*gx2^3]
%derivation of Jacobian matrix  
%J(gx,gw)=[0 0 1 1; gw1 gw2 gx1 gx2; 2*gw1*gx1 2*gw2*gx2 gx1^2 gx2^2; 2*gw1*gx1^2 3*gw2*gx2^2 gx1^3 gx2^3]
%its evaluation and solution of the linear system s=J\F at each iteration.
tol = 10^10*eps;	%Tolerance
xw0=[-0.5 0.5 1 2]';  	%Initial guesses
F0=[xw0(3)+xw0(4)-gint(1); xw0(3)*xw0(1) + xw0(4)*xw0(2); xw0(3)*xw0(1)^2 + xw0(4)*xw0(2)^2-gint(3); xw0(3)*xw0(1)^3 + xw0(4)*xw0(2)^3]
J0=[0 0 1 1; xw0(3) xw0(4) xw0(1) xw0(2); 2*xw0(3)*xw0(1) 2*xw0(4)*xw0(2) xw0(1)^2 xw0(2)^2; 3*xw0(3)*xw0(1)^2 3*xw0(4)*xw0(2)^2 xw0(1)^3 xw0(2)^3]
s=J0\-F0    %Solve the system
xw=xw0 + s 	%Update solution	
%Implement iterations to the tolerance tol.	
i=0;   
echo off
while (norm(xw-xw0)>tol & i<10)
    i=i+1;
    xw0=xw;
    F0=[xw0(3)+xw0(4)-gint(1); xw0(3)*xw0(1) + xw0(4)*xw0(2); xw0(3)*xw0(1)^2 + xw0(4)*xw0(2)^2-gint(3);xw0(3)*xw0(1)^3 + xw0(4)*xw0(2)^3]
    J0=[0 0 1 1; xw0(3) xw0(4) xw0(1) xw0(2); 2*xw0(3)*xw0(1) 2*xw0(4)*xw0(2) xw0(1)^2 xw0(2)^2; 3*xw0(3)*xw0(1)^2 3*xw0(4)*xw0(2)^2 xw0(1)^3 xw0(2)^3]
    s=J0\-F0		%Solve the system
    xw=xw0 + s 		%Update solution
    i
end
echo on
%Nodes and weights:
xg(1)=xw0(1)
xg(2)=xw0(2)
w1=xw0(3)
w2=xw0(4)
G2 = w1* exp(-xg(1)^2)+ w2* exp(-xg(2)^2)
%For other intervals, integrals ca be computed again using linear transformation as in
%Exa.8.2.1.