3 views (last 30 days)
Niles Martinsen on 18 Jul 2012
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.
Niles Martinsen on 18 Jul 2012
Note that f is a real variable running from -infinity to infinity.

Teja Muppirala on 18 Jul 2012
This integral can be calculated in MATLAB.
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.
Teja Muppirala on 18 Jul 2012
The INTEGRAL2 function was introduced in the most recent version of MATLAB, R2012a. If you don't have it, you can still do the same thing with DBLQUAD.

Niles Martinsen on 18 Jul 2012
Thanks for taking the time to write all that. That is very kind of you, and it gave me a big push in the right direction. Thanks!
I tried out your suggestion. I only have R2011a, so I have used dblquad, but there is a problem when I run
f = 3e2;
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 );
However I get an error message:
??? Error using ==>
fcnchk at 103
FUN must be a function,
a valid string
expression, or an
inline function object.
47
Do you know what is wrong?
Best, Niles.

Teja Muppirala on 18 Jul 2012
Ah, ok. Then you'll have to do it the hard way, like this:
f = 3e2;
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 );
Niles Martinsen on 18 Jul 2012
Thanks, it was extremely kind of you to push me forward. I learned a lot (and saved a lot of time!).
Best wishes!