0
$\begingroup$

There are many similar questions but they all seem to be integrating functions that are more involved than just a product of polynomials and the issue has more to do with the integrand than the command Integrate[] or the domain.

Edit: This is just an example of my more general problem. I'm writing an application that will need to integrate the product of two functions over arbitrarily complex domains. The user will not be able to validate all or any of the integrals. The functions will primarily be of the following types: Polynomials, Interpolation, and sometimes Gaussians

Example:

I'm integrating the product of two polynomials over a domain made up of many Disk[]s and I get a small imaginary part!?

The functions:

Clear[p1, p2, reg];
p1[x_, y_] := 2 Sqrt[6] x y;
p2[x_, y_] := 
  40 x^3 - 120 x y^2 - 120 x^3 (x^2 + y^2) + 360 x y^2 (x^2 + y^2) + 
  84 x^3 (x^2 + y^2)^2 - 252 x y^2 (x^2 + y^2)^2;

The Region:

reg = BooleanRegion[(#1 && ! #2 && #3 && #4) || 
(#5 && ! #6 && #7 && #8) || 
(#9 && ! #10 && #11 && #12) || 
(#13 && ! #14 && #15 && #16) || 
(#17 && ! #18 && #19 && #20) ||
(#21 && ! #22 && #23 && #24) || 
(#25 && ! #26 && #27 && #28) || 
(#29 && ! #30 && #31 && #32) || 
(#33 && ! #34 && #35 && #36) &, 
{Disk[{4.0625`, 0.`}, 1], 
Disk[{4.0625`, 0.`}, 0.1`], Disk[{2.8125`, 0.`}, 2.25`], 
Disk[{0, 0}, 5.0625`], Disk[{2.1875`, 1.0825317547305482`}, 1], 
Disk[{2.1875`, 1.0825317547305482`}, 0.1`], 
Disk[{2.8125`, 0.`}, 2.25`], Disk[{0, 0}, 5.0625`], 
Disk[{2.1875`, -1.0825317547305482`}, 1], 
Disk[{2.1875`, -1.0825317547305482`}, 0.1`], 
Disk[{2.8125`, 0.`}, 2.25`], Disk[{0, 0}, 5.0625`], 
Disk[{-0.15625`, 2.4356964481437338`}, 1], 
Disk[{-0.15625`, 2.4356964481437338`}, 0.1`], 
Disk[{-1.40625`, 2.4356964481437338`}, 2.25`], 
Disk[{0, 0}, 5.0625`], Disk[{-2.03125`, 3.518228202874282`}, 1], 
Disk[{-2.03125`, 3.518228202874282`}, 0.1`], 
Disk[{-1.40625`, 2.4356964481437338`}, 2.25`], 
Disk[{0, 0}, 5.0625`], Disk[{-2.03125`, 1.3531646934131856`}, 1], 
Disk[{-2.03125`, 1.3531646934131856`}, 0.1`], 
Disk[{-1.40625`, 2.4356964481437338`}, 2.25`], 
Disk[{0, 0}, 5.0625`], Disk[{-0.15625`, -2.4356964481437338`}, 1],
Disk[{-0.15625`, -2.4356964481437338`}, 0.1`], 
Disk[{-1.40625`, -2.4356964481437338`}, 2.25`], 
Disk[{0, 0}, 5.0625`], Disk[{-2.03125`, -1.3531646934131856`}, 1],
Disk[{-2.03125`, -1.3531646934131856`}, 0.1`], 
Disk[{-1.40625`, -2.4356964481437338`}, 2.25`], 
Disk[{0, 0}, 5.0625`], Disk[{-2.03125`, -3.518228202874282`}, 1], 
Disk[{-2.03125`, -3.518228202874282`}, 0.1`], 
Disk[{-1.40625`, -2.4356964481437338`}, 2.25`], 
Disk[{0, 0}, 5.0625`]}];

I do the Intgral: 

Integrate[p1[x, y] p2[x, y] , {x, y} ∈ reg]

and get: 

(* output -3.99966 - 4.69926*10^-32 I *)

It is a very small Imaginary part and I'm tempted to just take the real part and pretend nothing happened. When I use NIntegrate[] I get the same result but it complains about converging too slowly suggesting that Integrate[] is not using NIntegrate[] to do the evaluation.

What should I do? Can I trust the real part?

Improve this question
$\endgroup$
5
  • $\begingroup$ If I insert Rationalize into your reg calculation thus reg=BooleanRegion[..., Rationalize[{YourListOfDisk},0]] and then do your integral otherwise unchanged then I get exactly precisely 0. Can you reproduce this? Can you explain this? I don't think I have made a mistake, but that is always a possibility. And I certainly can't explain this. Well, not yet. $\endgroup$
    – Bill
    Apr 28, 2019 at 1:11
  • $\begingroup$ Interesting. Yes, I get 0 also. As LouisB points out the example has a symmetry that will result in a zero integral. $\endgroup$
    – c186282
    Apr 28, 2019 at 14:31
  • $\begingroup$ If you want exact results you must use exact numbers. At a minimum you must avoid machine precision calculations and use arbitrary precision instead. This can also be done by Integrate[p1[x, y]*p2[x, y], {x, y} \[Element] Rationalize[reg, 0]] $\endgroup$
    – Bob Hanlon
    Apr 28, 2019 at 15:43
  • $\begingroup$ I like this Rationalize[reg, 0] Idea. I also can force or encourage the creation of the domain to use rational numbers. But still, I do not see how an imaginary number slipped in from integrating simple real functions over a real domain. $\endgroup$
    – c186282
    Apr 28, 2019 at 16:06
  • $\begingroup$ I appreciate what you cannot see. There are dozens, hundreds and perhaps thousands of pages of code and algorithms hiding behind Integrate. Some of those algorithms are triggered by the presence of a decimal point. Others are triggered by the lack of a decimal point. Every time I see "Why do I get a microscopic complex part?" I think "get rid of the decimal points and see what happens. Since you got 4-epsilon +I*epsilon I expected to get 4 and was surprised I got 0. But I served my role by getting someone to then look more deeply at your problem and see why the result should be 0. $\endgroup$
    – Bill
    Apr 28, 2019 at 17:11

2 Answers 2

Reset to default
2
$\begingroup$

When we look at the individual disks and the 9 boolean sub-regions, we see that the domain of the integral is symmetric about the x-axis. Use something like

Graphics[{Opacity[1/25], disks}, Axes -> True]

to visualize the disks and the boolean sub-regions. The variable disks is, of course, the list of disks from the question. That visualization is not too good. You will want to use the sub-regions listed below.

Next, we can look at the integrand, q = Simplify[ p1[x,y] p2[x,y] ]. We see that the integrand is an odd function of $y$. The integral of an odd function over a symmetric region is zero. So, why did we get something non-zero?

If we look at the integral over each of the 9 individual sub-regions we find one of them is almost zero and the others are 4 pairs of additive inverses, +/- versions of each other. That is what we would expect from our analysis.

Here are the sub-regions:

r[1] = BooleanRegion[(#1 && ! #2 && #3 && #4) &, disks];
r[2] = BooleanRegion[(#5 && ! #6 && #7 && #8) &, disks];
r[3] = BooleanRegion[(#9 && ! #10 && #11 && #12) &, disks];
r[4] = BooleanRegion[(#13 && ! #14 && #15 && #16) &, disks];
r[5] = BooleanRegion[(#17 && ! #18 && #19 && #20) &, disks];
r[6] = BooleanRegion[(#21 && ! #22 && #23 && #24) &, disks];
r[7] = BooleanRegion[(#25 && ! #26 && #27 && #28) &, disks];
r[8] = BooleanRegion[(#29 && ! #30 && #31 && #32) &, disks];
r[9] = BooleanRegion[(#33 && ! #34 && #35 && #36) &, disks];

We can visualize the sub-regions with 

Show[DiscretizeRegion@BooleanRegion[Or, Array[r, 9]], Axes -> True]

We can look at the integral over the 5th sub-region and 9th sub-region as examples:

Integrate[q, {x, y} ∈ r[5]]]
Integrate[q, {x, y} ∈ r[9]]]
Integrate[q, {x, y} ∈ BooleanRegion[Or, {r[5], r[9]}]]
    (*   -2.44193*10^8 + 2.21955*10^-15 I
          2.44193*10^8
          3.46263*10^-8   *)

The non-zero result for the integral over the entire region is a numerical artifact.

Improve this answer
$\endgroup$
1
  • $\begingroup$ Thank you, I completely agree with your analysis, and I must admit I did not try to work out what I expected the integral to be, but I knew it should be real. My real issue is that I'm putting together an application that will need to do many integrals of this type over arbitrarily complex domains. The user will not be able to analyze the validity of each integral. Really that is the whole point of using Mathematica. I will add to the question the motivation for the question. $\endgroup$
    – c186282
    Apr 28, 2019 at 14:41
0
$\begingroup$

In order to avoid zero values,

NIntegrate[p1[x,y]*p2[x, y]+1,{x, y}\[Element] reg,WorkingPrecision ->20,AccuracyGoa->3] - 

NIntegrate[1, {x, y} \[Element] reg, WorkingPrecision -> 20,AccuracyGoal -> 3]

0.000033943852052196 and warnings

It takes time.

Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.