Opened 2 years ago

Last modified 6 months ago

#28431 new defect

integrate(1/(a-c*b),a) incorrect

Reported by: gutow Owned by:
Priority: minor Milestone: sage-wishlist
Component: calculus Keywords: integration, 1/x, maxima
Cc: Merged in:
Authors: gutow Reviewers:
Report Upstream: N/A Work issues:
Branch: Commit:
Dependencies: Stopgaps: use algorithm='sympy' in the integrate call

Status badges

Description (last modified by gutow)

This appears to be a maxima related error.

integrate(1/(a-b),a)

=> log(a-b), which is correct

integrate(1/(a-c*b),a)

=> log(c*b-a), which is incorrect.

Sympy.integrate(1/(a-c*b),a)
=> log(a-c*b), which is correct.

As SymPy? does this correctly can we get a quick fix by rerouting 1/x type integrals through SymPy?? I am not familiar with the details of the backend translation used by the integrate function.

I assume this should also be reported upstream, but I do not know how to verify that it really is an upstream error.

Change History (25)

comment:1 Changed 2 years ago by gutow

  • Description modified (diff)
  • Stopgaps set to use `algorithm='sympy'` in the integrate call

I believe I have isolated this to maxima or how Sage is using it.

Incorrect:

>integrate(1/(a-b*n),a, algorithm='maxima')
log(b*n - a)

Correct:

>integrate(1/(a-b*n),a, algorithm='sympy')
log(-b*n + a)

comment:2 Changed 2 years ago by gutow

  • Description modified (diff)

comment:3 Changed 22 months ago by charpent

  • Keywords maxima added
  • Priority changed from major to critical

The problem may be with our translation from maxima, which, by itself, behaves correctly:

charpent@p-202-021:~$ sage -maxima
;;; Loading #P"/usr/local/sage-python3/local/lib/ecl/sb-bsd-sockets.fas"
;;; Loading #P"/usr/local/sage-python3/local/lib/ecl/sockets.fas"
;;; Loading #P"/usr/local/sage-python3/local/lib/ecl/defsystem.fas"
;;; Loading #P"/usr/local/sage-python3/local/lib/ecl/cmp.fas"
Maxima 5.42.2 http://maxima.sourceforge.net
using Lisp ECL 16.1.2
Distributed under the GNU Public License. See the file COPYING.
Dedicated to the memory of William Schelter.
The function bug_report() provides bug reporting information.
(%i1) integrate(1/(a-b*c),a);
(%o1)                            log(a - b c)

which is correct. And, BTW:

(%i5) diff(integrate(1/(a-b*c),a),a);

(%o5) 1/(a-b*c)

also correct.

EDIT : This is not yet another effect of domain:complex:

(%i7) domain:complex;              

(%o7) complex
(%i8) integrate(1/(a-b*c),a);

(%o8) log(a-b*c)

which is again correct...

Furthermore:

sage: var("a, b, c")
(a, b, c)
sage: integrate(1/(a-c*b),a)
log(b*c - a)                           /* Wrong */
sage: integrate(1/(a-c*b),a).diff(a)
-1/(b*c - a)                           /* CORRECT ??? */

Two wrongs make one right ???

This is, IMHO, a seemingly innocuous sign of a serious problem in a basic ability. I raise the priority to critical.

Last edited 22 months ago by charpent (previous) (diff)

comment:4 follow-up: Changed 22 months ago by mantepse

I think that integrate(1/(x-a), x) should actually be log(abs(x-a)), because for x < a the expression log(x-a) doesn't make much sense, does it?

comment:5 follow-up: Changed 22 months ago by nbruin

I think the simplification

sage: 1/(a-b*c)
-1/(b*c - a)

explains this. In maxima:

(%i1) integrate( -1/(b*c-a), a);
(%o1)                            log(b c - a)

We're not doing anything wrong. It's really down to branch choice. (having log(abs(x-a)) as antiderivative doesn't play nice with complex analytic computations) Furthermore,

sage: log(b*c-a).diff(a)
-1/(b*c - a)

is just correct.

This is just us using a monomial ordering and a normalization of rational functions that leads maxima to choosing a branch for an antiderivative that doesn't match what you might expect to be chosen naively.

Last edited 22 months ago by nbruin (previous) (diff)

comment:6 in reply to: ↑ 4 ; follow-up: Changed 22 months ago by charpent

Replying to mantepse:

I think that integrate(1/(x-a), x) should actually be log(abs(x-a)),

Nope. Whereas

sage: with assuming (x, "real", a, "real"): (log(abs(x-a))).diff(x)
-1/(a - x)

In general,

sage: (log(abs(x-a))).diff(x)
-1/2*(a - x + conjugate(a) - conjugate(x))/abs(-a + x)^2

because for x < a the expression log(x-a) doesn't make much sense, does it?

It does. In the complex plane, and at the cost of choosing one among an infinity of solutions (branch choice).

Seen another way, the definite integral integrate(x-a),x,b, c} does make sense, and can be trivially computed (use the change of variable t=x-a), except if a+b<0 and a+c>0 : the function can't be integrated in x==a, and the integration constant has absolutely no reason to be the same on the two sides of this discontinuity. If you suppose that the two integration constants are the same, you obtain a Cauchy principal value, which is not a value of the integral.

comment:7 in reply to: ↑ 5 ; follow-up: Changed 22 months ago by charpent

Replying to nbruin:

I think the simplification

sage: 1/(a-b*c)
-1/(b*c - a)

explains this. In maxima:

(%i1) integrate( -1/(b*c-a), a);
(%o1)                            log(b c - a)

This is another (Maxima's) error.

We're not doing anything wrong.

So, why are we obtaining the opposite of the expected result (which is correctly given by Maxima) in the original case ?

Because we are passing a "simplified" operand, which sends us in Maxima's trap.

Not our fault ?

It's really down to branch choice. (having log(abs(x-a)) as antiderivative doesn't play nice with complex analytic computations)

Nope : unless I'm sorely mistaken (which is always possible), a branch choice problem would leave us with a constant (pure imaginary in this case, I believe) added to the correct value, not its opposite...

Last edited 22 months ago by charpent (previous) (diff)

comment:8 in reply to: ↑ 6 ; follow-up: Changed 22 months ago by charpent

Replying to charpent:

To illustrate, compare:

sage: integrate(1/x, x, -3, -1)
-log(3)
sage: integrate(1/x, x, -3, 1)
---------------------------------------------------------------------------
RuntimeError                              Traceback (most recent call last)

[ Snip... ]

ValueError: Integral is divergent.

comment:9 in reply to: ↑ 8 Changed 22 months ago by charpent

Replying to charpent:

Another illustration (to show that Maxima tries to give us a principal value, on whoch Sage chokes):

sage: %maxima integrate(1/x, x, -3, 1)
---------------------------------------------------------------------------

[ Snip ... ]

TypeError: error evaluating "integrate(1/x,x,-3,1)":
Error executing code in Maxima
CODE:
	integrate(1/x,x,-3,1);
Maxima ERROR:
	Principal Value
(%o67) -log(3)

comment:10 follow-up: Changed 22 months ago by mantepse

Sorry about my mistake, I thought that we were playing in the reals.

comment:11 in reply to: ↑ 7 ; follow-up: Changed 22 months ago by nbruin

Replying to charpent:

Nope : unless I'm sorely mistaken (which is always possible), a branch choice problem would leave us with a constant (pure imaginary in this case, I believe) added to the correct value, not its opposite...

log( b*c -a ) = log ( (-1) *(a-b*c))= log(-1) + log(a-b*c)

so, with certain branch choices they do differ by a constant and therefore have the same derivative; and the derivative has no complicated branch choices. Hence, log(a-b*c) and log(b*c-a) are both valid antiderivatives.

Last edited 22 months ago by nbruin (previous) (diff)

comment:12 in reply to: ↑ 10 Changed 22 months ago by charpent

Replying to mantepse:

Sorry about my mistake, I thought that we were playing in the reals.

Most of Sage's (and Maxima's) symbolic functions aim to play with complexes. Roughly speaking, they assume that SR element are complexes. Hence the difference in x.log().abs().diff() with and without assumptions about x...

For example:

sage: PolynomialRing(QQ,"x")(x^2+1).is_irreducible()
True
sage: PolynomialRing(QQbar,"x")(x^2+1).is_irreducible()
False
sage: PolynomialRing(SR,"x")(x^2+1).is_irreducible()
False

But note that:

sage: (x^2+1).is_irreducible()

[ Snip ... ]

AttributeError: 'sage.symbolic.expression.Expression' object has no attribute 'is_irreducible'

So these rules are not cast in stone...

comment:13 in reply to: ↑ 11 Changed 22 months ago by charpent

Replying to nbruin:

Replying to charpent:

Nope : unless I'm sorely mistaken (which is always possible), a branch choice problem would leave us with a constant (pure imaginary in this case, I believe) added to the correct value, not its opposite...

log( b*c -a ) = log ( (-1) *(a-b*c))= log(-1) + log(a-b*c)

so, with certain branch choices they do differ by a constant and therefore have the same derivative; and the derivative has no complicated branch choices. Hence, log(a-b*c) and log(b*c-a) are both valid antiderivatives.

Damn... You're right !

I stand corrected.

But still:

sage: with assuming(a, "real", b, "real", c, "real"): (1/(a - b*c)).integrate(a)
log(b*c - a)

is hard to swallow... Choosing I*pi as an integration constant in a problem involving real quantities is a stretch of imagination...

But that's a nice excuse for any undergrad committing an error of sign in an integration problem. Which, BTW, won't hold in the case of definite integration.

comment:14 Changed 22 months ago by gutow

So Sage needs to be clear about what conditions the solution holds for. It should also work correctly if forced to use real numbers and get the sign correct without a hidden constant of integration. If Maxima can handle this Sage should too.

comment:15 follow-up: Changed 22 months ago by gutow

Thinking on this a little more, another option that does not hide a mysterious choice of integration constants, such as log(-1), would be to provide multiple valid integrations. I'm thinking of behavior something like this:

>(1/(a-b)).integrate(a)
   [log(a-b),log(b-a)]

as when reporting other results with multiple valid solutions. This would leave the user to pick the correct one for their situation. If they are limited to real numbers then they would pick the version that provides positive input to the log function.

I am not sure about how this generalizes.

comment:16 in reply to: ↑ 15 ; follow-up: Changed 22 months ago by nbruin

Replying to gutow:

Thinking on this a little more, another option that does not hide a mysterious choice of integration constants, such as log(-1), would be to provide multiple valid integrations. I'm thinking of behavior something like this:

>(1/(a-b)).integrate(a)
   [log(a-b),log(b-a)]

as when reporting other results with multiple valid solutions. This would leave the user to pick the correct one for their situation. If they are limited to real numbers then they would pick the version that provides positive input to the log function.

I am not sure about how this generalizes.

It doesn't. This is exactly what integration constants express. You'd have to give an uncountable set of antiderivatives.

comment:17 in reply to: ↑ 16 Changed 22 months ago by gutow

Replying to nbruin:

Replying to gutow:

Thinking on this a little more, another option that does not hide a mysterious choice of integration constants, such as log(-1), would be to provide multiple valid integrations. I'm thinking of behavior something like this:

>(1/(a-b)).integrate(a)
   [log(a-b),log(b-a)]

as when reporting other results with multiple valid solutions. This would leave the user to pick the correct one for their situation. If they are limited to real numbers then they would pick the version that provides positive input to the log function.

I am not sure about how this generalizes.

It doesn't. This is exactly what integration constants express. You'd have to give an uncountable set of antiderivatives.

Duh! Of course, it doesn't generalize. I get my dunce hat for this one!

I am beginning to lean towards some kind of "domain" control as discussed on sage-devel. I can see algorithm='maxima_real' or 'sympy_real' as easy to use. The default should remain the complex domain. To facilitate usability by people like me (physical scientists and their students) the default domain needs to be stated up-front in the documentation for integrate.

As was discussed in the sage-devel thread, I am not clear this solves all the issues. However, it ought to help with the issue of log in the real domain requiring a positive real number input. I would expect it to trigger an error requiring an assumption about whether a > b or not.

comment:18 Changed 21 months ago by embray

  • Milestone changed from sage-8.9 to sage-9.1

Ticket retargeted after milestone closed

comment:19 Changed 17 months ago by mkoeppe

  • Milestone changed from sage-9.1 to sage-9.2

comment:20 Changed 11 months ago by mkoeppe

  • Milestone changed from sage-9.2 to sage-9.3

comment:21 Changed 6 months ago by kcrisman

  • Priority changed from critical to minor

I'm not even sure this is an error, but at any rate this is not critical, given that (up to integration constants and branch choice) this is correct.

As to documentation, there are a few places where we do mention this (IIRC). We haven't apologized for it, and I don't think we have to, but it's certainly true that it's confusing, and adding to the documentation in some visible places would be salutary.

comment:22 follow-up: Changed 6 months ago by kcrisman

Also, to Jonathan's suggestion, it seems like making the context in comment:6 more advertised would be useful. But since neither Sympy nor Maxima add the absolute value, I'm not sure if that aspect is really 'fixable'.

comment:23 in reply to: ↑ 22 Changed 6 months ago by gutow

Replying to kcrisman:

Also, to Jonathan's suggestion, it seems like making the context in comment:6 more advertised would be useful. But since neither Sympy nor Maxima add the absolute value, I'm not sure if that aspect is really 'fixable'.

From a practical usability perspective, I still think there should be a way to request the domain of the result, in addition to being much clearer about how the branch was chosen.

comment:24 Changed 6 months ago by kcrisman

That may be true, but if Maxima and Sympy don't, I'm not sure how we can :(

comment:25 Changed 6 months ago by mkoeppe

  • Milestone changed from sage-9.3 to sage-wishlist

Could open an issue with https://github.com/sympy/sympy expressing this feature request.

Note: See TracTickets for help on using tickets.