Faster code for double integration using Gauss-Legendre quadrature

Question

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

Answer 1

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