Evaluate advanced integral numerically

6 次查看（过去 30 天）
Niles Martinsen
Niles Martinsen 2012-7-18
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.

请先登录，再回答此问题。

采纳的回答

Teja Muppirala
Teja Muppirala 2012-7-18
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 个评论
Teja Muppirala
Teja Muppirala 2012-7-18
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.

请先登录，再进行评论。

更多回答（2 个）

Niles Martinsen
Niles Martinsen 2012-7-18
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.
Teja Muppirala
Teja Muppirala 2012-7-18
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)
  1 个评论
Niles Martinsen
Niles Martinsen 2012-7-18
Thanks, it was extremely kind of you to push me forward. I learned a lot (and saved a lot of time!).
Best wishes!

请先登录，再进行评论。

类别

在 Help Center 和 File Exchange 中查找有关 Data Type Identification 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by