%
%   ===   COMPARISON OF RESULTS FROM QUAD, GQUAD, GQUAD6   ===

> format long

> f=quad('sin',0,pi,1.e-5), % Integrate a simple function

f = 2.00000006453000

> % Use a ten point Gauss formula
> [b10,w10]=grule(10); f=gquad('sin',0,pi,1,b10,w10)

f = 2.00000000000000

> % Try a six point Gauss formula without having to generate
> % base points and weight factors.
> f=gquad6('sin',0,pi,1)

f = 1.99999999947727, % Accuracy is quite good even with a
                       % six point formula

> % Next consider a function having infinite slope at x=0
> f=3*quad('sqrt',0,1)
Recursion level limit reached in quad.  Singularity likely.
Recursion level limit reached in quad.  Singularity likely.
Recursion level limit reached in quad.  Singularity likely.
Recursion level limit reached in quad.  Singularity likely.
Recursion level limit reached in quad.  Singularity likely.
Recursion level limit reached in quad.  Singularity likely.
Recursion level limit reached in quad.  Singularity likely.
Recursion level limit reached in quad.  Singularity likely.

f =1.99998752859620

> f=3*gquad('sqrt',0,1,1,b10,w10)

f = 2.00026812880953, % The exact answer is 2

> [b20,w20]=grule(20); f=3*gquad('sqrt',0,1,1,b20,w20)

f = 2.00003589797091

> [b50,w50]=grule(50);  f=3*gquad('sqrt',0,1,100,b50,w50)

f = 2.00000000239878

> f=quad('log',0,1) ,   % Try to integrate a singular function

Warning: Log of zero

f =  -

> f= quad('log',eps,1) , % Stop just to the right of x=0
Recursion level limit reached in quad.  Singularity likely.
Recursion level limit reached in quad.  Singularity likely.
Recursion level limit reached in quad.  Singularity likely.
Recursion level limit reached in quad.  Singularity likely.
Recursion level limit reached in quad.  Singularity likely.
Recursion level limit reached in quad.  Singularity likely.
Recursion level limit reached in quad.  Singularity likely.
Recursion level limit reached in quad.  Singularity likely.

f=  -1.00421574636997

> % How well can a composite Gauss ten point rule do
> % using ten intervals?
> f= gquad('log',0,1,10,b10,w10)

f=  -0.99942637022162, % Exact answer is -1 so results are still
                         % rather poor. Perhaps a higher order
                         % will help.

> [b100,w100]=grule(100); f=gquad('log',0,1,1,b100,w100)

f= -0.99993748732245,  % This result is not much better

> % Let us try taking a large number of function values
> f=gquad('log',0,1,100,b100,w100)

f=  -0.99999937487322

> f=gquad6('log',0,1,1600), % Try the simpler six point formula

f=  -0.99999061950636

> %  The last several results show that we cannot just 'integrate
> %  through' a singularity by using many function values.

> format short