This is no longer a problem with Mathematica 10.* where I can use regions to do the integrations
Overview
I have a very large matrix I need to calculate where each element requires the integration of a product of two functions over some complicated domain. I would like to find a set of options for NIntegrate that will work well for all of my integrations. Most of the integrations take only seconds but others take up to 4 hours. Below is a histogram of timing for the integrations. It is clear that there are some very fast one and some very slow ones:
The following set of plots shows six integrands that took less than 15m along with a visualization of the domain. The plot label shows how long it took to integrate.
The following set of plots shows six integrands that took more than 15m along with a visualization of the domain. The plot label shows how long it took to integrate.
As you can see the two groups do not look much different however all of the ones that take a long time results in integrated values that are very close to zero ($\approx 10^{-16}$).
Details
Here are the functions for the forth plot in each set shown above:
A good one (fast):
rg={{-3.35513, 5.30513}, {-1.51554, 7.14471}};
funa=Function[{x, y}, E^(1/2 (-(0.39 + x/5) (100. (0.39 + x/5) -
6.12323*10^-15 (0. + y/5)) - (-6.12323*10^-15 (0.39 + x/5) +
50. (0. + y/5)) (0. + y/5)))]
funb=Function[{x, y}, E^(1/2 (-(-0.195 + x/5) (55.3571 (-0.195 + x/5) +
15.4647 (-0.562917 + y/5)) - (15.4647 (-0.195 + x/5) +
94.6429 (-0.562917 + y/5)) (-0.562917 + y/5)))]
region=Function[{x,y}, ((1/2 (4.54663 - 0.909327 x - 0.525 y) >= 0 &&
1/2 (0.954793 - 1.05 y) >= 0 &&
1/2 (-2.63705 + 0.909327 x - 0.525 y) >= 0 &&
1/2 (-2.63705 + 0.909327 x + 0.525 y) >= 0 &&
1/2 (0.954793 + 1.05 y) >= 0 &&
1/2 (4.54663 - 0.909327 x + 0.525 y) >=
0) || (1/2 (3.45544 - 0.909327 x - 0.525 y) >= 0 &&
1/2 (2.04599 - 1.05 y) >= 0 &&
1/2 (-0.454663 + 0.909327 x - 0.525 y) >= 0 &&
1/2 (-1.54586 + 0.909327 x + 0.525 y) >= 0 &&
1/2 (-0.136399 + 1.05 y) >= 0 &&
1/2 (2.36425 - 0.909327 x + 0.525 y) >=
0) || (1/2 (2.36425 - 0.909327 x - 0.525 y) >= 0 &&
1/2 (-0.136399 - 1.05 y) >= 0 &&
1/2 (-1.54586 + 0.909327 x - 0.525 y) >= 0 &&
1/2 (-0.454663 + 0.909327 x + 0.525 y) >= 0 &&
1/2 (2.04599 + 1.05 y) >= 0 &&
1/2 (3.45544 - 0.909327 x + 0.525 y) >=
0) || (1/2 (2.04599 - 0.909327 x - 0.525 y) >= 0 &&
1/2 (3.45544 - 1.05 y) >= 0 &&
1/2 (2.36425 + 0.909327 x - 0.525 y) >= 0 &&
1/2 (-0.136399 + 0.909327 x + 0.525 y) >= 0 &&
1/2 (-1.54586 + 1.05 y) >= 0 &&
1/2 (-0.454663 - 0.909327 x + 0.525 y) >=
0) || (1/2 (0.954793 - 0.909327 x - 0.525 y) >= 0 &&
1/2 (4.54663 - 1.05 y) >= 0 &&
1/2 (4.54663 + 0.909327 x - 0.525 y) >= 0 &&
1/2 (0.954793 + 0.909327 x + 0.525 y) >= 0 &&
1/2 (-2.63705 + 1.05 y) >= 0 &&
1/2 (-2.63705 - 0.909327 x + 0.525 y) >=
0) || (1/2 (-0.136399 - 0.909327 x - 0.525 y) >= 0 &&
1/2 (2.36425 - 1.05 y) >= 0 &&
1/2 (3.45544 + 0.909327 x - 0.525 y) >= 0 &&
1/2 (2.04599 + 0.909327 x + 0.525 y) >= 0 &&
1/2 (-0.454663 + 1.05 y) >= 0 &&
1/2 (-1.54586 - 0.909327 x + 0.525 y) >=
0) || (1/2 (-0.454663 - 0.909327 x - 0.525 y) >= 0 &&
1/2 (-1.54586 - 1.05 y) >= 0 &&
1/2 (-0.136399 + 0.909327 x - 0.525 y) >= 0 &&
1/2 (2.36425 + 0.909327 x + 0.525 y) >= 0 &&
1/2 (3.45544 + 1.05 y) >= 0 &&
1/2 (2.04599 - 0.909327 x + 0.525 y) >=
0) || (1/2 (-1.54586 - 0.909327 x - 0.525 y) >= 0 &&
1/2 (-0.454663 - 1.05 y) >= 0 &&
1/2 (2.04599 + 0.909327 x - 0.525 y) >= 0 &&
1/2 (3.45544 + 0.909327 x + 0.525 y) >= 0 &&
1/2 (2.36425 + 1.05 y) >= 0 &&
1/2 (-0.136399 - 0.909327 x + 0.525 y) >=
0) || (1/2 (-2.63705 - 0.909327 x - 0.525 y) >= 0 &&
1/2 (-2.63705 - 1.05 y) >= 0 &&
1/2 (0.954793 + 0.909327 x - 0.525 y) >= 0 &&
1/2 (4.54663 + 0.909327 x + 0.525 y) >= 0 &&
1/2 (4.54663 + 1.05 y) >= 0 &&
1/2 (0.954793 - 0.909327 x + 0.525 y) >=
0)) && ! ((-2.75 + x)^2 + (0. + y)^2 <=
0.316406 || (1.375 + x)^2 + (-2.38157 + y)^2 <=
0.316406 || (1.375 + x)^2 + (2.38157 + y)^2 <= 0.316406)]
A Bad one (slow):
rg={{2.9, 5.}, {-1.05, 1.05}};
funa=Function[{x, y}, E^(1/2 (-(0. + 100. (-0.39 + x/5)) (-0.39 + x/5)
- (0. + 50. (0. + y/5)) (0. + y/5)))]
funb=Function[{x, y}, 1.90476 (0. + y)]
region=Function[{x, y},
1/2 (4.54663 - 0.909327 x - 0.525 y) >= 0 &&
1/2 (0.954793 - 1.05 y) >= 0 &&
1/2 (-2.63705 + 0.909327 x - 0.525 y) >= 0 &&
1/2 (-2.63705 + 0.909327 x + 0.525 y) >= 0 &&
1/2 (0.954793 + 1.05 y) >= 0 &&
1/2 (4.54663 - 0.909327 x + 0.525 y) >=
0 && (-2.75 + x)^2 + (0. + y)^2 > 0.316406]
The ugly integral
NIntegrate[
Boole[region[x, y]] funa[x, y] funb[x, y],
{x, rg[[1, 1]], rg[[1, 2]]}, {y, rg[[2, 1]], rg[[2, 2]]},
Method->{"SymbolicPiecewiseSubdivision","ExpandSpecialPiecewise"->False,
Method->"LocalAdaptive"},
MaxRecursion->50,MinRecursion->20,
AccuracyGoal->7,PrecisionGoal->6
]
I have picked MinRecursion to be 20 to make sure that the smaller features of the domain are found. (This is just one example of a domain others are more involved with smaller items.) I have used "LocalAdaptive" because I seam to get fewer complaints about slow convergence than when I use "GlobleAdaptive" and finally I use "ExpandSpecialPiecewise"->False because there are a few cases where I get an error that I have exceeded the number of piecewise divisions (This dose not appear to be connected with the number of conditions in the domain but instead the divisions done by NIntegrate)
I have tried many combinations of options for NIntegrate but I have not found one that will work for all gazillion matrix elements. I realize it is not possible with the code give to test for all of my cases but I will give good ideas a run in my code and report back.
