Evaluate advanced integral numerically
6 visualizaciones (últimos 30 días)
Mostrar comentarios más antiguos
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.
Respuesta aceptada
Teja Muppirala
el 18 de Jul. de 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 comentario
Teja Muppirala
el 18 de Jul. de 2012
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.
Más respuestas (2)
Teja Muppirala
el 18 de Jul. de 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)
Ver también
Categorías
Más información sobre Data Type Identification en Help Center y File Exchange.
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!