Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

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: Histogram of <code>NIntegrate</code> timing

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. Example functions that takes less than 15m 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. Example functions that takes more than 15m 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.

share|improve this question
add comment

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Browse other questions tagged or ask your own question.