You do not say if the domain of integration is a simple one, perhaps a hypercube in R^6. It is not really relevant, but if the domain is more complex, that just makes it worse. And nothing that Meysam has posed will help you.
Numerical integration in a high number of dimensions is a nasty thing. And, sadly, 6 dimensions is what I would call moderately high. Why? Lets see what happens when you perform a single dimensional numerical integration. Commonly, tools like integral, or one of the quad tools, will require between 100 and 1000 function evaluations to compute an integral estimate in one dimension. That number of evals will depend on many things, but mainly I would think of the complexity of the integrand and the integration tolerance.
If your function is at all complicated, with locally high derivatives or even singularities, spikes, sharp bends, etc., then the integration will need to resolve where those problems lie, and use many points to pin down the best estimate of the integral. And naturally, the integration tolerance is important, since if you insist on a better estimate, then it will take more time and effort to get it for you. remember that any numerical integration, if it is to be accurate, must put many points in any region of importance. But adaptive numerical integrations do not know what your function looks like. All they do is evaluate the function, and then try to infer where the complexities lie from those function evals. Any such tool treats your objective asa simple black box. You pass numbers in, get numbers out. But inside the kernel, it is a black box that you cannot see into.
You might think that well, it makes sense then to just use a tool like Simpson's rule or trapezoidal rule. Just evaluate the function at a string of points. But then you must recognize that not much is happening at alost all of those points, so most of them are totally wasted function evals. So you need to accept the idea that 100-1000 points for a 1-d problem is about as good as you can do in a classical sense, IF the integration knows nothing about the integrand.
What happens in multiple dimensions? Now it gets exponentially nasty. An adaptive tool, in two dimensions will typically require at least on the order of 100^2 - 1000^2 function evaluations, so often at least 10000 function evaluations. In n dimensions (here n==6) we are now talking about potentially 100^6 = 10^12 function evals. And that may be on the low side. If your function is complicated, as I said, you may need to worry about something on the order of 1000^6 = 10^18 evals. This is the curse of dimensionality for nuemrical integration. 10^12 function evaluations is a TRILLION function evals. That will take some serious time to do.
So, what can you do? You do have several options. There is one numerical integration that does not suffer as badly from such a curse - Monte Carlo integration. In fact, the beauty of Monte Carlo is it can survive a 6 dimensional problem in a finite amount of time. And Monte Carlo is easy to perform too. The problem is that Monte Carlo integration will not give you a highly accurate estimate. In fact, you can estimate how the integration error for Monte Carlo performs with the number of evals. It tends to decrease as 1/sqrt(n), where n is the number of sample points, and that behavioral decrease will not depend on the number of dimensions!
A very important idea is often if you can do the integration analytically, in at least one of the necessary dimensions. If you can decrease the dimensions, of course then things get exponentially better.
Finally, you can consider tools like a tensor product Gaussian integration. But that typcially requires you know a very specific behavior of the integrand, and we have been told nothing in this respect. Sorry, but if you really want help, then I would need to know more about your problem.
So, without such specific information, I would strongly suggest you consider a Monte Carlo integration, which is actually trivially easy to compute. As a simple example, consider the function
Z = x1 + x2 + x3 + x^4 + x5 + x6
over the region from 0 to 1 in each of the 6 dimensions present. I've chosen it as an example, because it is a simple one where we know the answer to compare it to. Since we know
syms x
int(x,[0,1])
So 1/2 in one dimension. In 6 dimensions, the integral will be 3. A Monte carlo just has us sum the function values, divide by the area of the integration domain (here that ares is diff(limits)^6), and then divide by the number of sampels.
fun = @(X) sum(X,2);
intfun = @(N,limits) sum(fun(rand(N,6)*diff(limits) + limits(1)))/diff(limits)^6/N;
format long g
intfun(100,[0,1])
ans =
3.02889547585076
intfun(10000,[0,1])
ans =
3.00735201210059
intfun(10000000,[0,1])
ans =
2.99998533404415
So as we increased the sample size from 100 to 10 million, the error of integration goes down slowly. But you can see how trivially easy it is to write a Monte Carlo numerical integration.
I could have done far better if I knew the integrand itself of course.
