% Much of the hour was spent on the important topic of differentiating a
% program. Eventually, we returned to numerical quadrature, as discussed in
% the web notes for it.

% We are about to test composite quadratures, and need for that some test 
% functions. Rather than writing a function m-file, we use instead
% a very nice way available for defining functions on the fly, the
% socalled inline functions:

>> f = inline('x^5')
f =
     Inline function:
     f(x) = x^5
% We are ready to use a version of the program CompQNC from the book; this 
% version can be seen in the material on numerical quadrature on the web.

>> compqnc(f,0,1,3,10)
??? Error using ==> inline/feval
Error in inline expression ==> x^5
??? Error using ==> ^
Matrix must be square.

Error in ==> c:\courses\412\00\compqnc.m
On line 12  ==>    y = feval(fname,xx);

% Ok, the composite quadrature program expects to be given the name of a 
% function that works correctly even if the argument is a vector or a matrix.
% This means that we must make sure that all multiplications are done
% pointwise:

>> f = inline('x.^5')
f =
     Inline function:
     f(x) = x.^5
>> compqnc(f,0,1,3,10)
ans =
    0.1667
% we check how accurate this is by looking at more digits:

>> format long
>> ans
ans =
   0.16666875000000

% We increase the order of the method, from 3 to 4:
>> compqnc(f,0,1,4,10)
ans =
   0.16666759259259

%  ... to 5:
>> compqnc(f,0,1,5,10)
ans =
   0.16666666666667
% that is the exact value, even though we are using a rule that, offhand,
% is only exact for polynomials of degree < 5. What is happening???

% The answer: all NC rules are symmetric with respect to their interval,
% [0 .. 1], meaning that they give the same result for f(x) as for f(1-x).
% In particular, if  f(x) = -f(1-x), then they give the (correct) result 0.
% If now  f  is a quintic polynomial, then I can always write it as
%  f(x) = a*(x-1/2)^5 + q(x)
% with  q  some polynomial of degree < 5. The latter polynomial, the NC rule
% does exactly, by design, the former term, it does exactly, by symmetry.
% Therefore, the 5-point Newton-Cotes rule is actually of order 6.
%
% The same argument shows that any odd-point NC rule has order one more than
% the number of points used.

% Now we try the composite Gauss-Legendre rule. The m-point GL rule has maximal
% possible order, namely order 2m.
% (The following function m-file is also discussed and shown in the web notes
% on numerical quadrature).

% The QGL(3) has order 6, hence should be exact for our function, x^5 :

>> compqgl(f,0,1,3,10)
ans =
   0.16666666666667
% .., it is.

%  ...but even the QGL(2), only of order 4, does pretty well if we use ten
% subintervals:

>> compqgl(f,0,1,2,10)
ans =
   0.16666527777778

% The weights and nodes for these GL rules are usually provided by some
% table or m-file, like the following:
-
type glweights

  function [w,x] = GLWeights(m)
% [w,x] = GLWeights(m)
%
% w is a column m-vector consisting of the weights for the m-point
%  Gauss-Legendre rule.
% x is a column m-vector consisting of the abscissae.
% m is an integer that satisfies 2 <= m <= 6.

w = ones(m,1); x = ones(m,1);

switch m
case 1
case 2
   x(1)   = -0.577350269189626;
case 3
   w(1:2) = [5/9
             8/9              ];
   x(1:2) =[-0.774596669241483
             0                ];
case 4
   w(1:2) = [0.347854845137454
             0.652145154862546];
   x(1:2) =[-0.861136311594053
            -0.339981043584856];
case 5
   w(1:3) = [0.236926885056189
             0.478628670499366
             0.568888888888889];
   x(1:3) =[-0.906179845938644
            -0.538469310105683
             0                ];
case 6
   w(1:3) = [0.171324492379170
             0.360761573048139
             0.467913934572691];
   x(1:3) =[-0.932469514203152
            -0.661209386466265
            -0.238619186083197];
otherwise
   error('Sorry, I only know the m-point Gauss-Legendre rules for m=1:6.')
end

mo2 = 1:m/2;
x(m+1-mo2) = -x(mo2); w(m+1-mo2) =  w(mo2);

% the rest of the hour was spent discussing the use of the error term in
% a simple and a composite quadrature rule.


diary off