# Evaluate advanced integral numerically

6 views (last 30 days)

Show older comments

Hi

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.

### Accepted Answer

Teja Muppirala
on 18 Jul 2012

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.

##### 1 Comment

Teja Muppirala
on 18 Jul 2012

### More Answers (2)

Teja Muppirala
on 18 Jul 2012

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)

### See Also

### Categories

### Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!