Ask Your Question
1

definite integral returns error but indefinite followed by equivalent subtraction does not

asked 4 years ago

louisgag gravatar image

Searching the other questions on definite integrals I did not find the same issue, maybe I missed it. I use SageMath 9.2

When I define these variables:

var('x,y', domain='real')
var('c,p_i', domain='positive')
f_x_pi = 1/5*(3*c - 0.06)*(0.2969*sqrt(x/(c - 0.02)) - 0.126*x/(c - 0.02) - 0.3516*x^2/(c - 0.02)^2 + 0.2843*x^3/(c - 0.02)^3 - 0.1036*x^4/(c - 0.02)^4)
c_i_u = c - 0.02

This:

I = integrate(integrate((x^2+y^2),y,-f_x_pi,f_x_pi),x,0,c_i_u)

return a Maxima requested additional constraints error but this:

I_indef = integrate(integrate((x^2+y^2),y,-f_x_pi,f_x_pi),x)
I = I_indef(x=c_i_u) - I_indef(x=0)

or relying on algorithm='sympy', which is slower, does not.

I am wondering if that's a bug and whether I should report it.

Preview: (hide)

1 Answer

Sort by » oldest newest most voted
1

answered 4 years ago

Emmanuel Charpentier gravatar image

updated 4 years ago

Having floating point constants makes extremely difficult to grasp the structure of the integration bounds. Furthermore, this introduces some absurdities :

sage: bool((c-0.02)==(3*c-0.06))
False

So we replace them :

# Problem data :
var('x,y', domain='real')
var('c,p_i', domain='positive')
f_x_pi = 1/5*(3*c - 0.06)*(0.2969*sqrt(x/(c - 0.02)) - 0.126*x/(c - 0.02) - 0.3516*x^2/(c - 0.02)^2 + 0.2843*x^3/(c - 0.02)^3 - 0.1036*x^4/(c - 0.02)^4)
c_i_u = c - 0.02


# Laziness : find and replace the floating point constants
def GetNums(expr):
    if expr.operator() is None:
        if expr.is_numeric(): return set([expr])
        return set()
    return reduce(union, list(map(GetNums, expr.operands())))

L = [u for u in GetNums(f_x_pi) if not u.is_integer()]
# Keep the rationals and 0.02 (see above)
L.remove(1/2)
L.remove(1/5)
D = dict(zip(L, [var("c_{}".format(u)) for u in range(len(L))]))
# Bounds of the inner integration :
D1 = copy(D)
var("c_7")
D1.update({c-0.02:c_7,3*c-0.06:3*c_7})
foo = f_x_pi.subs(D1)

Now, we can compute an exact symbolic form of your expression. The computation of the outer integral by Maxima is indeed hypothesis-dependent :

sage: %time with assuming(c_7>0): Intp = (x^2+y^2).integral(y, -foo, foo).integral(x, 0, c_7).subs({D[u]:QQ(u) for u in D}).subs({c_7:c-1/50}).expand().simplify()
CPU times: user 362 ms, sys: 12 ms, total: 374 ms
Wall time: 269 ms
sage: Intp
379461086555604037/20207687500000000000*c^4 - 379461086555604037/252596093750000000000*c^3 + 1138383259666812111/25259609375000000000000*c^2 - 379461086555604037/631490234375000000000000*c + 379461086555604037/126298046875000000000000000
sage: %time with assuming(c_7<0): Intn = (x^2+y^2).integral(y, -foo, foo).integral(x, 0, c_7).subs({D[u]:QQ(u) for u in D}).subs({c_7:c-1/50}).expand().simplify()
CPU times: user 307 ms, sys: 8.02 ms, total: 315 ms
Wall time: 262 ms

But those expressions are equal :

sage: bool(Intp==Intn())
True

Computation via Sympy is much slower but doesn't need hypotheses :

sage: %time Ints = (x^2+y^2).integral(y, -foo, foo).integral(x, 0, c_7, algorithm="sympy").subs({D[u]:QQ(u) for u in D}).subs({c_7:c-1/50}).expand().simplify()
CPU times: user 12.4 s, sys: 68 ms, total: 12.5 s
Wall time: 12.5 s

and leads to the same value :

sage: bool(Intp==Ints)

This is not a bug : Maxima and Sympy use different algorithms (that's the whole point of having them both in |Sage...). The algorithm used by Maxima depends on huypotheses on the sign of c-0.02, which appears in a square root in the integration bounds...). The expresssions obtained in the two cases happen to be equal.

FWIW : another case can be seen when solving the differential equation a*f(x).diff(x,2)+b*f(x).diff(x)+c*f(x)+d==0 for f(x) : Maxma will request an hypothesis on the sign of b^2-4*a*c, and its answer will be the difference of two exponentials if positive and the sum of two cos functions of x if negative (both multiplied by an exponential function of x)). These answers turn out to be identical up to an application of de Moivre formula, as well as being equal to the (exponentially-expressed, IIRC) answer given by Sympy. Again, that's not a bug, just different ways to reach or express the same result.

HTH,

Preview: (hide)
link

Your Answer

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

Add Answer

Question Tools

1 follower

Stats

Asked: 4 years ago

Seen: 659 times

Last updated: Mar 11 '21