#3712 Symbolic integration may fail when called with numerical constants in function and/or limits.

[Emmanuel Charpentier]

This has been detected in Sage's version of Maxima (currently 5.44.0, running under ECL 20.4.24) in a Debian testing system tunning on x86-64 ; the imputation to maxima is done by calling Sage's Maxima interpreter in an Emacs session and reproducing the problem :

Consider the following (Maxima) code, attempting to find the (numerical values of) a couple of (symbolic) integrals defined with numerical constants :

E1:200;
nu:0.33;
rho1:7850;
Pi:%pi,numer;
ah:5;
e0:0.25;
h:1/ah;
em:1-sqrt(1-e0);
aa:Pi*z/(2*h)+Pi/(4); 
Ez:E1*(1-e0*(cos(aa)));
rhoz:rho1*(1-em*(cos(aa)));
Q11:Ez;
A11:integrate(Q11,z,-h/2,h/2),numer;

All of this, including the first integral, works. But the second one :

B11:integrate(Q11*z,z,-h/2,h/2),numer;

never returns, and needs an interrupt (which may fail) to get back control.

However, the integration with literal constant works : replace them by literal constants :

(%i23) L1:[a=200, b=1, c=0.25, d=7.853981633974483, e=0.7853981633974483,f=1/10];

The integral becomes :

(%i26) Q:a*(b-c*cos(d*z+e));
(%o26)                      a (b - c cos(d z + e))
(%i27) integrate(Q*z,z,-f,f);
                                                       2  2      2
          2 c d f sin(d f - e) + 2 c cos(d f - e) - b d  f  + b e
(%o27) a (--------------------------------------------------------
                                       2
                                    2 d
                                                                   2  2      2
                      2 c d f sin(d f + e) + 2 c cos(d f + e) - b d  f  + b e

                    - --------------------------------------------------------)
                                                   2
                                                2 d

And numerical substitution of these literal constants gives the sought result (obtained/checked by other methods) :

(%i28) subst(L1, integrate(Q*z,z,-f,f)); (%o28) 0.1739496967711208

The original Sage [ticket]((https://trac.sagemath.org/ticket/31351#ticket) shows that Sympy's integrator doesn't present this quirk.

[Stavros Macrakis]

1) Sage is probably misusing Maxima here -- for numerical integration, it should not be using integrate(...), numer at all.
2) Maxima integration was not designed to handle floating-point coefficients, though it does work in some cases.
3) The "numer" is not required to cause this bug, so this isn't a minimal example.
4) All that being said, there is a bug in Maxima, where different coefficients change the running time:

 integrate(z*(1-cos(7.85*z)),z,0,1) => 0 seconds
 integrate(z*(1-cos(7.8539*z)),z,0,1) => 3.2 seconds
 integrate(z*(1-cos(7.85398*z)),z,0,1) => 77 seconds

I suspected that this had something to do with rationalization of the coefficients (with increasing denominators), and sure enough, coefficients like float(1/9) cause it to run fast. But that's not all that's going on:

B11:'integrate(Q11*z,z,-h/2,h/2)$
B11r: scanmap('ratsimp,B11)$
ev(B11r,nouns) => result in 0.16 seconds

[Robert Dodier]

Tracing the minimal example integrate(z*cos(7.853981633974483*z + 0.7853981633974483),z) shows that integrate calls RISCHINT which eventually calls ROOTFAC which calls PGCD. With the constant = 7.853981633974483, PGCD is called over and over starting with a large number in its argument (not sure what it represents) and counting down. With a smaller number of digits e.g. 7.853, the counting starts at a smaller number and quickly gets to 1 where it stops, and ROOTFAC returns (1 1 1 1 1 ...) to its caller.

My first guess is that that the problem is related to factor(exp(%i*z*p/q) + 1) --> lots of terms, and the number of terms is related to p and q. I think there are other bug reports about trouble caused by trying to factor such expressions. Maybe it's possible to detect a potential problem in ROOTFAC and refuse to try it. However I'm not sure what ROOTFAC is trying to do.

[Robert Dodier]

FWIW I'm working on a bug fix in which ROOTFAC just gives up if it loops too many times; the lack of a result from ROOTFAC doesn't seem to prevent Maxima from getting a result by some other path (I can't tell that the result from ROOTFAC is being used even in cases when ROOTFAC succeeds). I'll try to push that soon.

[Emmanuel Charpentier]

Le mardi 16 février 2021 à 07:52 +0000, Robert Dodier a écrit :

• status: open --> closed
• Comment:
  Fixed by commit 2138ce9. I put a loop counter into ROOTFAC to just
  give up (returning the input unchanged) if the loops exceed
  factor_max_degree.

Sage's ticket status updated accordingly

Thanks a lot ! 

--
Emmanuel Charpentier

[bugs:#3712] Symbolic integration may fail when called with numerical
constants in function and/or limits.
Status: closed
Group: None
Labels: integrate floating point rational numer
Created: Sun Feb 07, 2021 04:49 PM UTC by Emmanuel Charpentier
Last Updated: Tue Feb 09, 2021 05:26 PM UTC
Owner: nobody
This has been detected in Sage's version of Maxima (currently 5.44.0,
running under ECL 20.4.24) in a Debian testing system tunning on x86-
64 ; the imputation to maxima is done by calling Sage's Maxima
interpreter in an Emacs session and reproducing the problem :
Consider the following (Maxima) code, attempting to find the
(numerical values of) a couple of (symbolic) integrals defined with
numerical constants :
E1:200;nu:0.33;rho1:7850;Pi:%pi,numer;ah:5;e0:0.25;h:1/ah;em:1-
sqrt(1-e0);aa:Piz/(2h)+Pi/(4); Ez:E1(1-e0(cos(aa)));rhoz:rho1(1-
em(cos(aa)));Q11:Ez;A11:integrate(Q11,z,-h/2,h/2),numer;
All of this, including the first integral, works. But the second one
:
B11:integrate(Q11z,z,-h/2,h/2),numer;
never returns, and needs an interrupt (which may fail) to get back
control.
However, the integration with literal constant works : replace them
by literal constants :
(%i23) L1:[a=200, b=1, c=0.25, d=7.853981633974483,
e=0.7853981633974483,f=1/10];
The integral becomes :
(%i26) Q:a(b-ccos(dz+e));(%o26) a (b - c cos(d z + e))(%i27)
integrate(Qz,z,-f,f); 2 2 2 2 c d f sin(d f - e) + 2 c cos(d f - e)
- b d f + b e(%o27) a (----------------------------------------------
---------- 2 2 d 2 2 2 2 c d f sin(d f + e) + 2 c cos(d f + e) - b d
f + b e - --------------------------------------------------------) 2
2 d
And numerical substitution of these literal constants gives the
sought result (obtained/checked by other methods) :
(%i28) subst(L1, integrate(Qz,z,-f,f)); (%o28) 0.1739496967711208
The original Sage
[ticket]((https://trac.sagemath.org/ticket/31351#ticket) shows that
Sympy's integrator doesn't present this quirk.
Sent from sourceforge.net because you indicated interest in
https://sourceforge.net/p/maxima/bugs/3712/
To unsubscribe from further messages, please visit
https://sourceforge.net/auth/subscriptions/

[Robert Dodier]

Fixed by commit 2138ce9. I put a loop counter into ROOTFAC to just give up (returning the input unchanged) if the loops exceed factor_max_degree.