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에서 Loops and Conditional Statements에 대해 자세히 알아보기
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!