Wha calculations are behind these statements?

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

Answer 1

[Jan]

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.

Answer 2

[kubota555]

Thank you for your most helpful advice. It enabled me to understand the principles behind the notation, concepts and application of row vectors, reverse order elements, transpositions, dot products and realmin/class values (at least for the code sample above.) Your change to the error magnitude calculation is also most helpful. I will also try your suggestion to allow the error to occur and manage the potential NaN for insight into program behaviour.

You provided very clear responses to my questions on such basic issues. A great teacher.

Thank you once again.