Evaluate advanced integral numerically

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]

Note that f is a real variable running from -infinity to infinity.

Answer 1

[Teja Muppirala]

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.

[Teja Muppirala]

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.

Answer 2

[Niles Martinsen]

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 );

Value = 2*(2*pi)*dblquad(@(r,z) r.*J(r,z,f),0,inf,0,inf,'abstol',1e-16)

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.

Error in ==> dblquad at
47
    quadf =
    fcnchk(quadf);

Do you know what is wrong?

Best, Niles.

Answer 3

[Teja Muppirala]

Ah, ok. Then you'll have to do it the hard way, like this:

myquad = @(fun,a,b,tol,trace,varargin)quadgk(@(x)fun(x,varargin{:}),a,b,'AbsTol',tol);

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 );
Value = 2*(2*pi)*dblquad(@(r,z) r.*J(r,z,f),0,Inf,0,Inf,1e-16,myquad)

[Niles Martinsen]

Thanks, it was extremely kind of you to push me forward. I learned a lot (and saved a lot of time!).

Best wishes!