Trouble Evaluating an Integral

4 次查看（过去 30 天）
Simon
Simon 2012-6-1
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');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

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

采纳的回答

Teja Muppirala
Teja Muppirala 2012-6-1
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更早的评论
Mike Hosea
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
Simon 2012-6-1
Thank you both for your answers.
It looks like I'll be upgrading to R2012a!

请先登录，再进行评论。

更多回答（0 个）

类别

在 Help Center 和 File Exchange 中查找有关 Stability Analysis 的更多信息

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by