Recommendations for evaluating a 6-D Integral
8 次查看(过去 30 天)
显示 更早的评论
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.
0 个评论
回答(3 个)
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.
0 个评论
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
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.
另请参阅
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!