Ask Your Question

Bug in integral?

asked 2013-09-02 09:20:49 +0100

petropolis gravatar image


def h(x): return 1/((x + 1)^(1/3) + 1)
hint = integral(h(x), (x, -6, -2))

Exception ValueError: ValueError('negative number to a fractional power not real',) in 'sage.gsl.integration.c_ff'


int(1/((x+1)^(1/3) + 1), x = -6..-2): 
1.562667835 - 1.105160054 I
edit retag flag offensive close merge delete

1 Answer

Sort by ยป oldest newest most voted

answered 2013-09-02 09:54:28 +0100

tmonteil gravatar image

updated 2013-09-02 10:04:17 +0100

If you consider that you work on the real field, Maple answers a complex number, which is just weird. Sage's Symbolic Ring does the same artificial choice since:

sage: h(-3.3)
0.408623955144076 - 0.281398327127290*I

So, i agree that there is a kind of bug here, see also this question.

If you consider that you work on the complex field, your function is not well defined since there are 3 different choices for the cubical root, leading to different integrals.

If you want the cube root to be the real one, your integral will go to -Infinity as x approaches -2. And Sage can see this. Let us specify the real cube root as follows for the negative numbers, so that you will get the real cube root (not a complex one):

sage: cuberoot(x) = -((abs(x))^(1/3))


sage: h(x) = 1/(cuberoot((x + 1)) + 1)
sage: integral_numerical(h, -6, -2)
(-inf, nan)
edit flag offensive delete link more


Not sure if Maple's answer is weird. Mathematica gives the same answer:

petropolis gravatar imagepetropolis ( 2013-09-02 11:10:46 +0100 )edit

It depends on your needs. Where is the integral coming from ? What is its meaning ? When you write `(-8)^(1/3)`, do you except the real cube root of `-8` (that is `-2`) or do you expect a complex cube root (like `1.00 + 1.73*I`) ?

tmonteil gravatar imagetmonteil ( 2013-09-02 14:27:54 +0100 )edit

tmonteil, thanks for your comments! (1) Working on the complex field I assume that the principal root is chosen by default; Maple and Mathematica seem to do this also. Do I have to declare this in Sage explicitly? (2) When choosing the real-valued root the Cauchy principal value is -oo, as you say. How do I assure that Sage uses the real-valued root in the above integral? (3) Now, what is the meaning of the Maple and Mathematica solutions? I think they compute an indefinite integral (antiderivative) and than take the difference at the end points. This leads to the value 1.562667835 - 1.105160054*I which both programs give. Let's try to do this with Sage: def h(x): return 1/((x + 1)^(1/3) + 1) hint = integral(h(x), x) (hint(x=-6) - hint(x=-2)).n() sage: 2.76637207714311 - 1.394212302715

petropolis gravatar imagepetropolis ( 2013-09-03 07:49:47 +0100 )edit

For your question (2), you can use the function `cuberoot()` i defined in my answer. For your question (1), you should know that the other cube roots come from the real one by multiplying by `j` or `j^2` as explained in [this answer]( Concerning your question (3), i have no idea on how does Maple and Mathematica work here (i do not have them). Concerning the choice of a "principal root" in general, you have to know that there are some issues on chosing one root globally, if a curve winds around zero, chosing the "principal root" let you jump from one branch of the cube root to another one on the way... This may explain the different results, and why i do not like how the `Symbolic Ring` handles this.

tmonteil gravatar imagetmonteil ( 2013-09-04 14:47:28 +0100 )edit

Your Answer

Please start posting anonymously - your entry will be published after you log in or create a new account.

Add Answer

Question Tools


Asked: 2013-09-02 09:20:49 +0100

Seen: 431 times

Last updated: Sep 02 '13