MATLAB: Wha calculations are behind these statements

adaptivegautschiquadrature

As a beginner with Matlab, I'm trying to understand the numerical calculations behind three statements in the following code extract, and I'm hoping someone could explain them for me: – –
    %   Based on "adaptlob" by Walter Gautschi.
    %   Ref: W. Gander and W. Gautschi, "Adaptive Quadrature Revisited", 1998.
    %   http://www.inf.ethz.ch/personal/gander
    %   Copyright 1984-2004 The MathWorks, Inc. 
    %   $Revision: 1.8.4.7 $  $Date: 2004/12/06 16:35:09 $
    f = fcnchk(funfcn);
    % Initialize with 13 function evaluations.
    c = (a + b)/2;
    h = (b - a)/2;
    s = [.942882415695480 sqrt(2/3) .641853342345781 1/sqrt(5) .236383199662150];
    x = [a c-h*s c c+h*fliplr(s) b];
    y = feval(f,x,varargin{:}); y = y(:).';
    fcnt = 13;
    % Increase tolerance if refinement appears to be effective.
    Q1 = (h/6)*[1 5 5 1]*y(1:4:13).';
    Q2 = (h/1470)*[77 432 625 672 625 432 77]*y(1:2:13).';
    s = [.0158271919734802 .094273840218850 .155071987336585 ...
         .188821573960182  .199773405226859 .224926465333340];
    w = [s .242611071901408 fliplr(s)];
    Q0 = h*w*y.';
    r = abs(Q2-Q0)/abs(Q1-Q0+realmin(class(Q0)));
-
-
In a mathematical calculation sense, what are the numbers being used to calculate these terms:
1. the expression for w?
2. the expression for Q0?
3. the term realmin(class(Q0))?
I understand the remaining code. Any assistance appreciated. Thanks

Best Answer

You will get an error in this code.
f = fcnchk(funfcn);
c = (a + b)/2;
The values of a and b are not defined or provided as an input. This can only work, if you have M-files called a.m and b.m, but this would be a very bad choice of names.
1. the expression for w?
w = [s .242611071901408 fliplr(s)];
What exactly is your question concerning this line? w is defined as a row vector, which contains the elements of s, the value 0.242... (the leading zero can be omitted, but it is recommended to avoid this for clarity) and finally the elements of s in reversed order. See help fliplr.
    r = abs(Q2-Q0)/abs(Q1-Q0+realmin(class(Q0)));
2. the expression for Q0?
    Q0 = h*w*y.';
Q0 is h times w times the transposed of y. When w and y are vectors, this is the dot product.
3. the term realmin(class(Q0))?
Did you read the documentation already?
doc(realmin)
doc(class)
The smallest possible number of the class of Q0 is created, e.g. 2.2251e-308 if Q0 is a double. This avoid the division by 0 if Q1 and Q0 have exactly the same value. But this can still fail, when Q1 equals Q0+realmin. Better:
(abs(Q1 - Q0) + realmin(class(Q0))
This is still extremely large and the result is dominated by random rounding errors. Therefore I'd prefer to get a NaN from a division by zero, such that the problem is clear.
Related Question