Ask Your Question
1

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

asked 2021-03-09 12:16:33 +0100

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.

edit retag flag offensive close merge delete

1 Answer

Sort by ยป oldest newest most voted
1

answered 2021-03-11 15:07:44 +0100

Emmanuel Charpentier gravatar image

updated 2021-03-11 16:02:45 +0100

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,

edit flag offensive delete link more

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: 2021-03-09 12:16:33 +0100

Seen: 495 times

Last updated: Mar 11 '21