Faster code for double integration using Gauss-Legendre quadrature

I came with a the following code to evaluate a double integral using Gauss Legendre quadrature in MatLab

m=100;

%generate weights and abscissas
[wx,xx]=leg(-1,1,m);
[wtheta,xtheta]=leg(0,2.*pi,m);

%define function
psi=@(x,theta) hypergeom(-3./4,1./2,x.^2.*exp(1i.*theta));

%integrate with respect to x
intx=zeros(1,m);
for num=1:m
intx(num)=sum(wx.*psi(xx,xtheta(num)));
end

sum(wtheta.*intx)


I defined the function leg(x1,x2,m) in a different script to generate the weights and abscissas and I just call it in my code. My MatLab code runs significanlty slowly compared to NIntegrate of Mathematica.

I'd like to make my code faster because my idea of using MatLab is that it is faster than Mathematica. Is there any way I can make my code run faster?

Attached is a profile summary when i ran the code

• Your profiling reveals that mupadmex consumes most of the time. I suggest to go deeper in your profiling. This routine is a MEX files, it means that it is compiled from a low-level language like Fortran or C which is faster so the optimization is already done. That also means that you can't improve this routine because it is already compiled. But I don't understand where it is called in your code, that's where you should look into. Why do you need this routine and what is its purpose ? – Coriolis Aug 8 '16 at 19:13
• I just clicked the run and time button because I wanted to know how long the code ran :) – PhilCsar Aug 8 '16 at 22:39
• Are you using symbolic algebra inside your function leg? If so, that might be the reason for your algorithm to be much slower. – nicoguaro Aug 10 '16 at 13:28
• @nicoguaro no it doesnt – PhilCsar Aug 10 '16 at 14:50
• Your hypergeom is slow. See math.stackexchange.com/q/478052/24717 – Memming Sep 9 '16 at 9:20

The reason might be that your method is not as sophisticated. I would bet that NIntegrate uses an adaptive method and splits the integrals in parts which are then refined according to error estimates. If you do not improve your algorithm, you probably won't beat it. As someone pointed out, mupadex is the costly part. This comes from the number of function evaluations (specifically hypergeom). You can reduce them by some divide-and-conquer algorithm. It is not a good idea to use such a high number of Gauss-nodes anyway. A better one is as mention to resort to a lower order integrator and refine the grid. If you want code, look at the matlab function integral2Calc.m for the actual algorithm (integral2t) or at its wrapper integral2.m

In case anyone wants to run the code above, here is a working leg function.

function [x, w] = leg(a,b,N)
% Generates the abscissa and weights for a Gauss-Legendre quadrature.
beta = (1:N-1)./sqrt(4*((1:N-1)).^2-1);
[w,x] = eig(diag(beta,-1)+diag(beta,1));
absc = diag(x);
wght = 2*w(1,:)'.^2;

% Linear map from[-1,1] to [a,b]
x=(a*(1-absc)+b*(1+absc))/2;
w=wght;

x=x';
w=w';

end

• This method of generating the abscissae and weights of GL quadrature is totally new to me, I'm only aware of generating them from the roots of the Legendre function iteratively (e.g. Newton's method). What is the relationship here between the eigenvalues/vectors of this matrix of off-diagonals? – cbcoutinho Dec 7 '16 at 10:19
• @cbcoutinho: this appears to be a specific implementation (for the case of the Gauss-Legendre quadrature) of the more general algorithm proposed by Golub and Welsch (1969) for quadrature rules for orthogonal polynomials. See their paper: ams.org/journals/mcom/1969-23-106/S0025-5718-69-99647-1 – okrzysik Dec 7 '16 at 10:41
• @okrzysik: Thanks for the quick reply and the link to the article. I'm curious if there has been a comparison between these two (or more?) methods of determining the abscissae and weights. Should there be a new question for that? – cbcoutinho Dec 7 '16 at 12:17
• @cbcoutinho: each method has its pros and cons (some of which they share); there are plenty of more recent articles, introducing sophisticated approaches for computing quadrature nodes and weights, that discuss some of the pros and cons. I'm unsure if any such posts, on this site, discuss comparisons between the methods. Perhaps you might like to look at some of the links given in this question on MathOverflow: mathoverflow.net/questions/203863/… – okrzysik Dec 7 '16 at 12:38