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);endfigure;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*phiResendfigure;plot(a1,RectResult); title('RectResult');%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Best Answer