Recommendations for evaluating a 6-D Integral

2 次查看（过去 30 天）
Simon
Simon 2012-6-6
I've have to integrate a 6-D function that doesn't seem to have an analytical solution.
Can anyone recommend a decent numerical integration/sampling technique or script that might be able to handle this. I've implemented a script that runs a 3-D quadrature on a function that calls another 3-D quadrature to obtain its value but it's using a for loop for the first quadrature and so it's very slow. It takes about 8 hours to evaluate the integral - the quadratures meshgrid is 70x70x70 and evaluates through each point individually.
I have access to a machine with quite a bit of memory (24GB) that should be able to handle fairly large matrices.
Any help at all would be greatly appreciated.

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

回答（3 个）

the cyclist
the cyclist 2012-6-6
I believe that the most efficient algorithms for the evaluation of high-dimensional integrals use Monte Carlo techniques. I recall that there is a good discussion of these techniques (and why they are superior in higher dimensions) in the book Numerical Recipes. I don't have my copy handy, so I can't be more specific.
Teja Muppirala
Teja Muppirala 2012-6-6
Here's a simple idea. If your integrand is more or less smooth. Evaluate it on a number of different coarse grids, and then extrapolate to estimate the answer.
If your grid spacing is h, then the error should converge with order h^2.
For example:
Integrate F = exp(x+y+z+u+v+w) from 0 to 1 on all 6 variables.
The analytic solution is (exp(1) - 1)^6. We will estimate this using 4 different grids, and then extrapolate to get the final result.
F = @(x,y,z,u,v,w) exp(x+y+z+u+v+w);
I = []; %Estimate ov each grid
h = []; %Spacing of each grid
for n = 5:8;
pts = linspace(0,1,n+1);
pts = conv(pts,[1 1]/2,'valid');
[X,Y,Z,U,V,W] = ndgrid(pts);
V = F(X,Y,Z,U,V,W);
h(end+1) = 1/n;
I(end+1) = sum(V(:))/(n^6);
end;
P = polyfit(h,I,numel(h));
EstValue = P(end) %The constant term is what we want
TrueValue = (exp(1)-1)^6
RelativeError = (EstValue - TrueValue)/TrueValue
  1 个评论
Simon
Simon 2012-6-6
That's an interesting approach.
Unfortunately, I seems like my function may not be smooth enough for it (or maybe I'm doing something wrong).
The variable names are r, theta, phi, s, gamma and delta. Their limits are zeros to infinity, pi. 2pi, infinity, pi and 2pi, respectively.
(I've approximated infinity as 70 in this case!).
Is dividing by n^6 still the right thing to do considering that I changed the intervals?
a = [1;0;0]
F = @(r,theta,phi,s,psi,delta) exp(-(a(1)-r.*sin(theta).*cos(phi)+s.*sin(psi).*cos(delta)).^2).*...
exp(-(a(2)-r.*sin(theta).*sin(phi)+s.*sin(psi).*sin(delta)).^2).*...
r.*sin(theta).*s.*sin(psi).*exp(-(a(3)-r.*cos(theta)+s.*cos(psi)).^2);
I = []; %Estimate ov each grid
h = []; %Spacing of each grid
for n = 5:20;
%Set Grid Limits and Generate Grid
rlimit = 70;
rpts = linspace(0,rlimit,n+1);
rpts = conv(rpts,[1 1]/2,'valid');
thetapts = linspace(0,pi,n+1);
thetapts = conv(thetapts,[1 1]/2,'valid');
phipts = linspace(0,2*pi,n+1);
phipts = conv(phipts,[1 1]/2,'valid');
slimit = 70;
spts = linspace(0,slimit,n+1);
spts = conv(spts,[1 1]/2,'valid');
gammapts = linspace(0,pi,n+1);
gammapts = conv(gammapts,[1 1]/2,'valid');
deltapts = linspace(0,2*pi,n+1);
deltapts = conv(deltapts,[1 1]/2,'valid');
[X,Y,Z,U,V,W] = ndgrid(rpts,thetapts,phipts,slimit,gammapts,deltapts);
%Evaluate Function
V = F(X,Y,Z,U,V,W);
h(end+1) = 1/n;
I(end+1) = sum(V(:))/(n^6);
end;
P = polyfit(h,I,numel(h));
EstValue = P(end) %The constant term is what we want

请先登录，再进行评论。

Teja Muppirala
Teja Muppirala 2012-6-6
Assuming that "gamma" and "psi" are the same thing, are you sure your integrand is correct? Look at when
theta = psi = pi/2
phi = delta = pi
r = s
Then your integrand turns into (after a bit of symbolic math)
s^2*exp(1)
Just try it: F(70,pi/2,pi,70,pi/2,pi)
you get exp(-1)*70^2
F(7000,pi/2,pi,7000,pi/2,pi)
you get exp(-1)*7000^2
So as r and s get big, it blows up to infinity.
  1 个评论
Simon
Simon 2012-6-6
Hmmm... I'll have a look at it closely tomorrow and get back to you.
And you are correct, gamma and psi are the same thing.

请先登录，再进行评论。

类别

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

标签

Community Treasure Hunt

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

Start Hunting!

Translated by