# Real-valued integral over real variable returns complex result(s)

Posted 9 years ago
15324 Views
|
7 Replies
|
6 Total Likes
|
 Why does this real-valued integral over a real variable return a complex-valued result? fSine[x_] := Piecewise[{{1/(\[Pi]*Sqrt[1 - x^2]), -1 < x < 1}, {0, True}}]; Integrate[fSine[x]*fSine[1/10 - x], {x, -2, 2}]; N[%] 0.444207 - 1.815 I Sanity check: Integrate[fSine[x]*fSine[1/10 - x], {x, -2, 2}, Assumptions -> x \[Element] Reals]; N[%] 0.444207 - 1.815 I Compare to NIntegrate result: NIntegrate[fSine[x]*fSine[1/10 - x], {x, -2, 2}] 0.444207 - 5.75862*10^-11 I  Why is there a non-zero imaginary part for both Integrate and NIntegrate? Why are the two answers different? Attachments:
7 Replies
Sort By:
Posted 9 years ago
 This might be useful:  fSine[x_] := HeavisidePi[x/2]/(\[Pi]*Sqrt[1 - x^2]); Chop@N[Integrate[fSine[x]*fSine[1/10 - x], {x, -2, 2}], 20] and N[Quiet[NIntegrate[fSine[x]*fSine[1/10 - x], {x, -2, 2}, WorkingPrecision -> 45, MaxRecursion -> 1000]], 20] both give: 0.44420658153737745178 These work both with the Heaviside definition. We can achieve something similar with the original Piecewise definition. fSine[x_] := Piecewise[{{1/(\[Pi]*Sqrt[1 - x^2]), -1 < x < 1}, {0, True}}]; Integrate[fSine[x]*fSine[1/10 - x], {x, -2, 2}, GenerateConditions -> False, PrincipalValue -> True]; Chop@N[%, 20] which gives: 0.44420658153737745178 ,too.Cheers,MarcoPS: Oh, yes, and SetPrecision[N[Simplify@Integrate[fSine[x]*fSine[1/10 - x], {x, -2, 2}, GenerateConditions -> False, PrincipalValue -> True], 60], Infinity] suggests that the imaginary part is really zero...
Posted 9 years ago
 Seems like Mathematica should clean up after itself better, wrt the small imaginary part. Ill send a bug reportat least the symbolic version should get fixed. One trouble here is that the underlying 1/(\[Pi]*Sqrt[1 - x^2]) has branch points at ±1. If you ask for an approximate result, it's very hard to automatically deduce that the very large imaginary values just beyond the branch points don't matter. And then, the numerical approximation itself involves complex numbers. So, even if the cancellation seems perfect, Mathematica throws in a 0. I just to be sure that further processing does not assume that the imaginary part is guaranteed to be zero. Use Chop.I'd guess that the branch cut issue that apparently troubles the symbolic calculation must seem like a whack-a-mole game to the Mathematica developers: I don't believe there is any known algorithm here, so it must involve a lot of heuristics. Sometimes those will be wrong. FWIW this example comes out of finding the PDF for a sum of random sine variablesconvolution. I have a fallback method that doesnt do the convolutions. That kind of problem often works better in the domain of "generalized functions". By feeding HeavisidePi to Integrate, you're hinting that's what you want. However, using operations like Convolve and FourierTransform in place of explicit integrals, you'd push it harder in that direction. Integrate tries to give you a symbolic solution, but it doesn't really apply a particular definition of "integration" to a particular definition of "function". That's also heuristic.
Posted 9 years ago
 One trouble here is that the underlying 1/([Pi]*Sqrt[1 - x^2]) has branch points at ±1. Indeed. I had thought that there would be a small chance to avoid that problem, at least on the numerical side, by defining my fSine function using < rather than <=, so that its value at ±1 is defined to be zero. But no such luck. However, using operations like Convolve and FourierTransform in place of explicit integrals, you'd push it harder in that direction. Thanks for the suggestion. Unfortunately, in my particuar case, Mathematica punts with Convolve even for the sum of two sines (convolving fSine with itself once); the characteristic function for the sum of N sines is BesselJ[0, \[Omega]]^N and Mathematica punts for N > 2 when doing the inverse Fourier transform, attempting to get the PDF for sums of more than two sines. Apparently my only hope is to get some plots of the first few PDFs, say, for up to six summed sines. My approach for that is to compute the characteristic functions at discrete points, then inverse DFT. Its the best I can think of.Finally: Dont quadrature packages usually have ways of dealing with infinities? A quick look at my old copy of Numerical Recipes, as well as my even older copy of the HP-15C Advanced Functions Handbook (now theres a reference for you) indicate that this is true. It seems that possibly Mathematica has gotten itself into a bit of trouble by relying on its symbolics here. The help page for NIntegrate indicates that it looks for singularities at the boundaries. In any case, I doubt that a straight-up integrator would never return a result with an imaginary component.
Anonymous User
Anonymous User
Posted 9 years ago
 The symbolic indefinite integral is quite a beast, and includes imaginary terms. This is common: mathematics often forces Mathematica to express a real result in terms of complex numbers. Unfortunately, in this difficult case, Mathematica also seems to misplace a branch cut when computing the definite integral. That's fairly common. You could report it as a bug, and maybe someone will figure out a way to fix this case.The break points in Piecewise are at the singularities. That's tough for both symbolic and numerical methods. Here, NIntegrate appears by default to do some symbolic processing at the singularities to help out. It apparently has to go complex here, and the poorly cancelled imaginary part results. If you set Method -> {Automatic, "SymbolicProcessing" -> False} you'll get a real result of poor accuracy after complaints of slow convergence.You can get a better symbolic result by using HeavisidePI instead of Piecewise here: fSine[x_] := HeavisidePi[x/2]/(\[Pi]*Sqrt[1 - x^2]) Integrate[fSine[x]*fSine[1/10 - x], {x, -2, 2}] // N NIntegrate[fSine[x]*fSine[1/10 - x], {x, -2, 2}] (*0.444207 + 0. I*) (*NIntegrate::slwcon: Numerical integration converging too slowly; suspect one of the following: singularity, value of the integration is 0, highly oscillatory integrand, or WorkingPrecision too small. >>*) (*NIntegrate::ncvb: NIntegrate failed to converge to prescribed accuracy after 9 recursive bisections in x near {x} = {-0.8985}. NIntegrate obtained 0.4371297669288421 and 0.0061209610092931716 for the integral and error estimates. >>*) (*0.43713*) HeavisidePI tends to be more friendly to calculus than Piecewise. You still get the +0. I in the numerical result from the symbolic integral because it involves complex numbers whose imaginary parts cancel. If that's a problem, use Chop or Re when you know the result is real. The NIntegrate result isn't so good with HeavisidePI. It's similar to the result using Piecewise without symbolic processing: real but innaccurate. Some Method setting might work better: there are so many because there's no predicting what method might work best.
Posted 9 years ago
 Thanks for the tip on preferring HeavisidePi over Piecewise.Seems like Mathematica should clean up after itself better, wrt the small imaginary part. Ill send a bug reportat least the symbolic version should get fixed.FWIW this example comes out of finding the PDF for a sum of random sine variablesconvolution. I have a fallback method that doesnt do the convolutions.
Posted 9 years ago
 But in my example, NIntegrate did not return a non-zero imaginary part. It is small, and one may suspect machine precision problems, but It is not zero. Having a non-zero imaginary part messes up any number of subsequent uses of the result.I would further maintain that your result, while numerically correct, is at least marginally wrong or inconvenient because it returns any imaginary part. If I perform nearly any other integral such as this, such as In[6]:= NIntegrate[Sin[x], {x, 0, 1}] Out[6]= 0.459698 I dont get + 0. I appended to the end.I think both of the examples I give, Integrate and NIntegrate, are wrongbuggy. Should I file a bug report?
Posted 9 years ago
 You can get the same answer as Nintegrate by finding the indefinite integral and then using limits int = Integrate[fSine0[x]*fSine0[1/10 - x], x]; (Limit[int, x -> 2] - Limit[int, x -> -2]) // N (*0.444207 + 0. I*) Which is even more accurate than Nintegrate. I do not know now why it gave 0.444207 - 1.815 I . may be branch cut issue.