MATLAB: Evaluate advanced integral numerically

difficultnumerical integration

Hi
I am trying to evaluate the following integral in MATLAB: http://www.scribd.com/doc/100400549/mwe
However I haven't had any success, nor in a competing software (I'm not sure I'm allowed to mention its name). Do you guys know if it is even possible to evaluate such an integral in MATLAB? I have spent so many hours on this by now..
I would be very happy to get some feedback.
Best, Niles.

Best Answer

This integral can be calculated in MATLAB.
The integrand as written in your link is:
 J = @(x,y,z,f) sqrt(x.^2+y.^2)./sqrt(x.^2+y.^2+z.^2).*exp((-x.^2-y.^2-z.^2)/2) ./ ((f - 1e6*sqrt(x.^2+y.^2+z.^2)).^2 + (1e4)^2/4 )
But we're going to change this a bit.
Step 1. Convert it into cylindrical coordinates (r,theta,z).
 J = @(r,theta,z,f) r./sqrt(r.^2+z.^2).*exp((-r.^2-z.^2)/2) ./ ((f - 1e6*sqrt(r.^2+z.^2)).^2 + (1e4)^2/4 )
The limits of integration are now, r = [0, inf], z = [-inf,inf], theta = [0,2*pi].
Step 2. We don't want to integrate in 3d if we don't have to. Note that since theta is not in the equation, you just multiply by 2*pi, and you don't have to worry about it. Now you are just left with a 2d integral in r and z. Also note it is symmetric in z, so you can actually just integrate z = [0,inf] and multiply that by 2. Finally, since dx*dy*dz = r*dr*dtheta*dz, you need to thrown in an extra r.
Step 3 Put all that together and actually do the integral (You'll need to set the tolerance fairly small).
J = @(r,z,f) r./sqrt(r.^2+z.^2).*exp((-r.^2-z.^2)/2) ./ ((f - 1e6*sqrt(r.^2+z.^2)).^2 + (1e4)^2/4 )
f = 3e6;
Value = 2*(2*pi)*integral2(@(r,z) r.*J(r,z,f),0,inf,0,inf,'abstol',1e-16)
This is all it takes to evaluate the integral for a given value of f.
Step 4 Make a plot to show how the integral varies as a function of f
Value = [];
figure;
for f = linspace(-1e7,1e7,201);
  Value(end+1) = 2*(2*pi)*integral2(@(r,z) r.*J(r,z,f),0,inf,0,inf,'AbsTol',1e-16), semilogy(f,Value(end),'x');
  hold on;
  drawnow; 
end
Some of the points will take a while, but this loop could also be done quicker by parallel processing using a PARFOR if needed.
Related Question