Trouble Evaluating an Integral
2 次查看(过去 30 天)
显示 更早的评论
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 个评论
采纳的回答
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 个评论
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).
更多回答(0 个)
另请参阅
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!