Revision history [back]

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)
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,