Evaluating integral numerically using quadgk

1 view (last 30 days)
Hi
I am trying to evaluate the following integral numerically in MatLAB: http://www.scribd.com/doc/100400549/mwe
However, the sqrt(x^2 + y^2) in the numerator has changed to x. Here is my attempt (based on this thread: http://mathworks.com/matlabcentral/answers/43929-evaluate-advanced-integral-numerically#answer_53917):
clc
clear all
myquad = @(fun,a,b,tol,trace,varargin)quadgk(@(x)fun(x,varargin{:}),a,b,'AbsTol',tol);
sigma = 1;
kappa = 1e4;
B = 1e6;
beta = 1.0;
f = 1e6;
fct = @(x,y,z,f) (...
(x)/sqrt(x.^2+y.^2+z.^2)).*(exp((-x.^2-y.^2-z.^2)/(2*sigma^2)))./ ((f - B*beta*sqrt(x.^2+y.^2+z.^2)).^2 + kappa^2/4 );
result = triplequad(@(x,y,z) fct(x,y,z,f), -Inf, Inf, -Inf, Inf, -Inf, Inf, 1e-16, myquad);
However this gives me a warning: "Rank deficient, rank = 0, tol = NaN.". I'm not sure what else to do here. Am I using triplequad correctly?
Best wishes, Niles.
  2 Comments
Niles Martinsen
Niles Martinsen on 14 Aug 2012
I noticed that the warning only enters when the limits are from -infinity to infinity. If I change them to 0..infinity, then the warning disappears. Unfortunately the integrand is odd, so I have to use -inf..inf.
Teja Muppirala
Teja Muppirala on 14 Aug 2012
If that changes to an x, then the integral is zero by symmetry.

Sign in to comment.

Answers (1)

Mike Hosea
Mike Hosea on 14 Aug 2012
Make sure that you are using elementwise operations. See the first / in fct? It should be ./ instead. Check all the other multiplications and divisions. If it's only a scalar multiplication or division, then you don't need to change it, but in fact I find it easier just to use all .* and ./ rather than * or /.

Community Treasure Hunt

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

Start Hunting!