Trouble Evaluating an Integral

Question

Simon 2012-6-1

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/40033-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');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Teja Muppirala 2012-6-1

1
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/40033-trouble-evaluating-an-integral#answer_49665

在 MATLAB Online 中打开

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)

4 个评论
显示 2更早的评论隐藏 2更早的评论

Mike Hosea 2012-6-1

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 2012-6-1