Trouble Evaluating an Integral

I'm having trouble evaluating an integral in MATLAB.

The integrand is defined in spherical coordinates, (r,theta,phi) and the limits are r->{0,inf}, theta ->{0,pi}, phi->{0,2*pi}.

integrand = @(r,theta,phi) exp(-(a(1)-(r.*sin(theta).*cos(phi))).^2).*...
    exp(-(a(2)-(r.*sin(theta).*sin(phi))).^2).*...
    r.*sin(theta).*exp(-(a(3)-(r.*cos(theta))).^2)

I want to evaluate this integral for different values of the 3-D vector, a.

In the example below, I keep the second and third elements of a to be fixed a to be 0 and only vary the first element, a(1).

I evaluate it using two different methods, the QUAD function in MATLAB and just rectangular method that I implemented myself.

My question is, why does quad perform so poorly compared to the few lines of code that implements the rectangular method. There are large discontinuities in the resulting plots at certain points. Is there anything I can do to improve it?

Is there a method that does better than both? You'll notice that the plot produced by the rectangular method is slightly asymmetrical and this shouldn't be the case.

%%%%%%%%%%%%QUAD FUNCTION%%%%%%%%%%%

clear all;
a1 = linspace(-50,50,101);

for i = 1:numel(a1)

      Countdown = numel(a1)-i     %Countdown iterations of for loop
      a =[a1(i); 0; 0];           %Setup the 3d vector a for this iteration

      %Define the integrand for this iteration
      integrand = @(r,theta,phi) exp(-(a(1)-(r.*sin(theta).*cos(phi))).^2).*...
      exp(-(a(2)-(r.*sin(theta).*sin(phi))).^2).*...
      r.*sin(theta).*exp(-(a(3)-(r.*cos(theta))).^2);

      %Perform quadrature
      QuadResult(i) = triplequad(integrand,0,100,0,pi,0,2*pi);

end

figure;
plot(a1,QuadResult); title('Quad Result');
pause

%%%%%%%%%%%RECTANGULAR METHOD%%%%%%%%%%%%%%%

clear all;
a1 = linspace(-50,50,101);

for i = 1:numel(a1)

      Countdown = numel(a1)-i     %Countdown iterations of for loop
      a =[a1(i); 0; 0];           %Setup the 3d vector a for this iteration

      %create 3-D rectangular grid
      rLimit = 100;
      rRes = rLimit/200;
      thetaRes = pi/200;
      phiRes = 2*pi/200;
      [r theta phi] = meshgrid(0:rRes:rLimit,0:thetaRes:pi,0:phiRes:2*pi);
      %Express Grid a 3 vectors
      rVec = r(:)';
      ThetaVec = theta(:)';
      PhiVec = phi(:)';

      %Define the integrand
      integrand = @(r,theta,phi) exp(-(a(1)-(r.*sin(theta).*cos(phi))).^2).*...
      exp(-(a(2)-(r.*sin(theta).*sin(phi))).^2).*...
      r.*sin(theta).*exp(-(a(3)-(r.*cos(theta))).^2);

      %Sample from integrand at grid coordinates
      Samples = integrand(rVec,ThetaVec,PhiVec);

      %Multiple each sample by the area of the grid cell
      RectResult(i)=sum(Samples)*rRes*thetaRes*phiRes

end

figure;
plot(a1,RectResult); title('RectResult');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Answer 1

[Teja Muppirala]

If you are using R2012a, you can try using INTEGRAL3

In any case, since you are integrating in spherical coordinates, you need to multiply your integrand by r.^2.*sin(theta)

a = [1 0 0];

integrand = @(r,theta,phi) ( exp(-(a(1)-(r.*sin(theta).*cos(phi))).^2).*...
    exp(-(a(2)-(r.*sin(theta).*sin(phi))).^2).*...
    r.*sin(theta).*exp(-(a(3)-(r.*cos(theta))).^2) ).*r.^2.*sin(theta);

integral3(integrand,0,inf,0,pi,0,2*pi)

[Simon]

Thanks for the reply.

Regarding multiplying the integrand by r.^2.*sin(theta), I am actually already doing that (that's where the r.*sin(theta)) comes from. There was a denominator of r so that cancelled out the square power in the numerator.

I have MATLAB 2011b. I am considering upgrading. Could you tell me if integral3 performs better? For example, what does your version return for:

a = [1 0 0];

integrand = @(r,theta,phi) ( exp(-(a(1)-(r.*sin(theta).*cos(phi))).^2).*...

exp(-(a(2)-(r.*sin(theta).*sin(phi))).^2).*...

r.*sin(theta).*exp(-(a(3)-(r.*cos(theta))).^2) );

integral3(integrand,0,inf,0,pi,0,2*pi)

Thanks again for your help.

[Teja Muppirala]

It returns 4.692434420087396 in about 0.4 seconds.

[Mike Hosea]

Just to add an answer to the "why" question, triplequad (and integral3) work to estimate and control the error. If you want to compare performance, always make sure you're paying for and getting the same accuracy in both cases. It's meaningless to note that one method is a lot faster than another when the faster method has a lot less accuracy. If you don't need more than a few digits of accuracy and want the integrator to be faster, set the tolerance(s) of the integrator accordingly, e.g. integral3(integrand,0,inf,0,pi,0,2*pi,'AbsTol',1e-5,'RelTol',1e-3).

[Simon]

Thank you both for your answers.

It looks like I'll be upgrading to R2012a!