% started off trying out matlab's numerical quadrature commands quad and quad8
% Both use adaptive quadrature
>> help quad

 QUAD   Numerically evaluate integral, low order method.
    Q = QUAD('F',A,B) approximates the integral of F(X) from A to B to
    within a relative error of 1e-3 using an adaptive recursive
    Simpson's rule.  'F' is a string containing the name of the
    function.  Function F must return a vector of output values if given
    a vector of input values.  Q = Inf is returned if an excessive
    recursion level is reached, indicating a possibly singular integral.
 
    Q = QUAD('F',A,B,TOL) integrates to a relative error of TOL.  Use
    a two element tolerance, TOL = [rel_tol abs_tol], to specify a
    combination of relative and absolute error.
 
    Q = QUAD('F',A,B,TOL,TRACE) integrates to a relative error of TOL and
    for non-zero TRACE traces the function evaluations with a point plot
    of the integrand. 
 
    Q = QUAD('F',A,B,TOL,TRACE,P1,P2,...) allows parameters P1, P2, ...
    to be passed directly to function F:   G = F(X,P1,P2,...).
    To use default values for TOL or TRACE, you may pass in the empty
    matrix ([]).
 
    See also QUAD8, DBLQUAD.

%  note the possibility of getting a trace, which, in this case, means,
% a picture showing the actual points on the function used by the command

>> f = inline('sqrt(x)')
f =
     Inline function:
     f(x) = sqrt(x)

% integrate the function over the interval [0..1]. Note that the integrand
% has infinite slope at the left endpoint, making it difficult to integrate
% numerically.

>> quad(f,0,1,1e-6,1)
Warning: Recursion level limit reached in quad. Singularity likely.
> In C:\MATLABR11\toolbox\matlab\funfun\quad.m (quadstp) at line 102
  In C:\MATLABR11\toolbox\matlab\funfun\quad.m (quadstp) at line 141
  In C:\MATLABR11\toolbox\matlab\funfun\quad.m (quadstp) at line 141
  In C:\MATLABR11\toolbox\matlab\funfun\quad.m (quadstp) at line 141
  In C:\MATLABR11\toolbox\matlab\funfun\quad.m (quadstp) at line 141
  In C:\MATLABR11\toolbox\matlab\funfun\quad.m (quadstp) at line 141
  In C:\MATLABR11\toolbox\matlab\funfun\quad.m (quadstp) at line 141
  In C:\MATLABR11\toolbox\matlab\funfun\quad.m (quadstp) at line 141
  In C:\MATLABR11\toolbox\matlab\funfun\quad.m (quadstp) at line 141
  In C:\MATLABR11\toolbox\matlab\funfun\quad.m (quadstp) at line 141
  In C:\MATLABR11\toolbox\matlab\funfun\quad.m (quadstp) at line 141
  In C:\MATLABR11\toolbox\matlab\funfun\quad.m at line 72
Warning: Recursion level limit reached 40 times.
> In C:\MATLABR11\toolbox\matlab\funfun\quad.m at line 80
ans =
    0.6667

% let's let up on the tolerance a bit, using 1e-3 rather than 1e-6:

>> quad(f,0,1,1e-3,1)
Warning: Recursion level limit reached in quad. Singularity likely.
> In C:\MATLABR11\toolbox\matlab\funfun\quad.m (quadstp) at line 102
  In C:\MATLABR11\toolbox\matlab\funfun\quad.m (quadstp) at line 141
  In C:\MATLABR11\toolbox\matlab\funfun\quad.m (quadstp) at line 141
  In C:\MATLABR11\toolbox\matlab\funfun\quad.m (quadstp) at line 141
  In C:\MATLABR11\toolbox\matlab\funfun\quad.m (quadstp) at line 141
  In C:\MATLABR11\toolbox\matlab\funfun\quad.m (quadstp) at line 141
  In C:\MATLABR11\toolbox\matlab\funfun\quad.m (quadstp) at line 141
  In C:\MATLABR11\toolbox\matlab\funfun\quad.m (quadstp) at line 141
  In C:\MATLABR11\toolbox\matlab\funfun\quad.m (quadstp) at line 141
  In C:\MATLABR11\toolbox\matlab\funfun\quad.m (quadstp) at line 141
  In C:\MATLABR11\toolbox\matlab\funfun\quad.m (quadstp) at line 141
  In C:\MATLABR11\toolbox\matlab\funfun\quad.m at line 72
Warning: Recursion level limit reached 8 times.
> In C:\MATLABR11\toolbox\matlab\funfun\quad.m at line 80
ans =
    0.6667
% still difficult. Look at the answer in more detail:
>> format long
ans
ans =
   0.66666250953207
% the graph shows high density of points near the left endpoint.


% Now try out quad8 which is supposedly more efficient.
>> quad8(f,0,1,1e-3,1)
Warning: Recursion level limit reached in quad8. Singularity likely.
> In C:\MATLABR11\toolbox\matlab\funfun\quad8.m (quad8stp) at line 129
  In C:\MATLABR11\toolbox\matlab\funfun\quad8.m (quad8stp) at line 123
  In C:\MATLABR11\toolbox\matlab\funfun\quad8.m (quad8stp) at line 123
  In C:\MATLABR11\toolbox\matlab\funfun\quad8.m (quad8stp) at line 123
  In C:\MATLABR11\toolbox\matlab\funfun\quad8.m (quad8stp) at line 123
  In C:\MATLABR11\toolbox\matlab\funfun\quad8.m (quad8stp) at line 123
  In C:\MATLABR11\toolbox\matlab\funfun\quad8.m (quad8stp) at line 123
  In C:\MATLABR11\toolbox\matlab\funfun\quad8.m (quad8stp) at line 123
  In C:\MATLABR11\toolbox\matlab\funfun\quad8.m (quad8stp) at line 123
  In C:\MATLABR11\toolbox\matlab\funfun\quad8.m (quad8stp) at line 123
  In C:\MATLABR11\toolbox\matlab\funfun\quad8.m (quad8stp) at line 123
  In C:\MATLABR11\toolbox\matlab\funfun\quad8.m at line 70
Warning: Recursion level limit reached 2 times.
> In C:\MATLABR11\toolbox\matlab\funfun\quad8.m at line 77
ans =
   0.66666663842017
%  still difficulty, but the answer is slightly better; in fact, way better
%  than the tolerance we specified. This may point to a not very sophisticated
%  assignment of the tolerance when subintervals are considered.

% notice the points are not uniformly spaced, in fact, there
% are strange gaps. 
% These are due to use of Gauss rules which are open rules. In fact, from
% the graph, we guess that 3-point Gauss-Legendre rules are being used.

% Let's try an even less smooth function, namely a step function:

>> f = inline('x>.3')
f =
     Inline function:
     f(x) = x>.3
>> quad8(f,0,1,1e-3,1)
Warning: Recursion level limit reached in quad8. Singularity likely.
> In C:\MATLABR11\toolbox\matlab\funfun\quad8.m (quad8stp) at line 129
  In C:\MATLABR11\toolbox\matlab\funfun\quad8.m (quad8stp) at line 123
  In C:\MATLABR11\toolbox\matlab\funfun\quad8.m (quad8stp) at line 124
  In C:\MATLABR11\toolbox\matlab\funfun\quad8.m (quad8stp) at line 123
  In C:\MATLABR11\toolbox\matlab\funfun\quad8.m (quad8stp) at line 123
  In C:\MATLABR11\toolbox\matlab\funfun\quad8.m (quad8stp) at line 124
  In C:\MATLABR11\toolbox\matlab\funfun\quad8.m (quad8stp) at line 124
  In C:\MATLABR11\toolbox\matlab\funfun\quad8.m (quad8stp) at line 123
  In C:\MATLABR11\toolbox\matlab\funfun\quad8.m (quad8stp) at line 123
  In C:\MATLABR11\toolbox\matlab\funfun\quad8.m (quad8stp) at line 124
  In C:\MATLABR11\toolbox\matlab\funfun\quad8.m (quad8stp) at line 123
  In C:\MATLABR11\toolbox\matlab\funfun\quad8.m at line 70
Warning: Recursion level limit reached 2 times.
> In C:\MATLABR11\toolbox\matlab\funfun\quad8.m at line 77
ans =
   0.69991207492835
% the integral should be .7, so, this is pretty good!

% The book also considers integration of spline functions, but does it wrong,
% I think. Better to do the following: Suppose we have in hand the command
% fnint(pp) that returns to us the piecewise polynomial function that is the
% integral or primitive of the piecewise polynomial function  f  in pp.
% Then we can compute  integral f(x) dx  over [a..b] in one statement:
%   diff( ppval( fnint(pp) , [a b] ))
% Let's try that out:

>> x = linspace(0,2*pi,5);
>> sp = spline(x,sin(x));
>> diff(ppval(fnint(sp),[0,2*pi]))
ans =
     0
% That's remarkable! The exact answer! (Actually, it's just symmetry.)
% Let's try it instead on the interval 
>> diff(ppval(fnint(sp),[0,pi]))
ans =
  2.09439510239320
%  .. while the exact value should be  2 .
%
% Check this against the spline interpolant to 5 points on the interval
% [0 .. pi]. That makes the breakpoint space half of what it was, hence,
% since the error in spline interpolation behaves like spacing^4, this should
% cut the error by 2^4 = 16.

>> x = linspace(0,pi,5);
>> sp = spline(x,sin(x));
>> diff(ppval(fnint(sp),[0,pi]))
ans =
  2.00455975498442
%  compare this with 1/16 of the earlier error:
>> .094395/16
ans =
   0.00589968750000
% ... pretty close!

%%  the rest of the class time was spent on matrices in matlab, especially
% matrix multiplication. See the notes under the story so far.

%  get some matrices to play with:
>> a = reshape(1:10,2,5)
a =
     1     3     5     7     9
     2     4     6     8    10
>> b = reshape(1:6,3,2)
b =
     1     4
     2     5
     3     6
%  well, to multiply a with b, I need a to have only three columns.
% Here is one way to get rid of the last two columns:
>> a(:,[4 5]) = []
a =
     1     3     5
     2     4     6
% now I compute the outer product of the first column of a with the first
% row of b , to make clear that I don't get a scalar but a full matrix:

>> a(:,1)*b(1,:)
ans =
     1     4
     2     8
% Of course, if I multiply the row with the column, i.e., the other way around
% I get the scalar or inner product:
>> b(1,:)*a(:,1)
ans =
     9
%%  Then we generated the Hilbert matrix, the one with A(i,j) = 1/(i+j-1).
%% We don't want to use any loops !!!!

function H = hilb(n)
%HILB   Hilbert matrix.
%   HILB(N) is the N by N matrix with elements 1/(i+j-1),
%   which is a famous example of a badly conditioned matrix.
%   See INVHILB for the exact inverse.
%
%   This is also a good example of efficient MATLAB programming
%   style where conventional FOR or DO loops are replaced by
%   vectorized statements.  This approach is faster, but uses
%   more storage.

%   C. Moler, 6-22-91.
%   Copyright (c) 1984-98 by The MathWorks, Inc.
%   $Revision: 5.6 $  $Date: 1997/11/21 23:28:56 $

%   I, J and E are matrices whose (i,j)-th element
%   is i, j and 1 respectively.

J = 1:n;
J = J(ones(n,1),:);
I = J';
E = ones(n,n);
H = E./(I+J-1);

%%%%%   this can be done slightly more swiftly now. E.g., for n=3:
>> n = 3;
>> J = repmat(1:n,n,1);
>> H = 1./(J'+J-1)
H =
   1.00000000000000   0.50000000000000   0.33333333333333
   0.50000000000000   0.33333333333333   0.25000000000000
   0.33333333333333   0.25000000000000   0.20000000000000

% I contacted Cleve Moler about this. His present solution is:
>> [I,J] = ndgrid(1:n);
>> H = 1./(I+J-1)
% ... very neat!
diary off