Evaluate advanced integral numerically
3 次查看(过去 30 天)
显示 更早的评论
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.
采纳的回答
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
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 个)
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)
另请参阅
类别
在 Help Center 和 File Exchange 中查找有关 Function Creation 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!