Different results for numerical and symbolic integration

200 views
Skip to first unread message

Pablo Vitoria

unread,
Oct 23, 2022, 6:14:09 PM10/23/22
to sage-devel
I am using Sage 9.7 running in Arch linux over WSL2

I get different results for an integral using numerical integration (which seems to agree with Mathematica) and symbolic integration:

numerical_integral(x^2*(log(x))^4/((x+1)*(1+x^2)),0,1)
(0.07608822217400527, 1.981757967172001e-07)

integral(x^2*(log(x))^4/((x+1)*(1+x^2)),x,0,1)
6*I*polylog(5, I) - 6*I*polylog(5, -I) + 765/64*zeta(5)
integral(x^2*(log(x))^4/((x+1)*(1+x^2)),x,0,1).n()
0.440633136273036

TB

unread,
Oct 23, 2022, 10:44:28 PM10/23/22
to sage-...@googlegroups.com
One more data point that uses maxima:
sage: sage.calculus.calculus.nintegral(x^2*(log(x))^4/((x+1)*(1+x^2)),
x, 0, 1)
(0.07608822235542824, 4.544650350490897e-10, 231, 0)

This can also be computed via
f = x^2*(log(x))^4/((x+1)*(1+x^2))
f.nintegral(x,0,1)

Regards,
TB
> --
> You received this message because you are subscribed to the Google
> Groups "sage-devel" group.
> To unsubscribe from this group and stop receiving emails from it, send
> an email to sage-devel+...@googlegroups.com
> <mailto:sage-devel+...@googlegroups.com>.
> To view this discussion on the web visit
> https://groups.google.com/d/msgid/sage-devel/4e75aaeb-5416-4ef8-826b-675c24fb297fn%40googlegroups.com <https://groups.google.com/d/msgid/sage-devel/4e75aaeb-5416-4ef8-826b-675c24fb297fn%40googlegroups.com?utm_medium=email&utm_source=footer>.

Frédéric Chapoton

unread,
Oct 24, 2022, 3:19:46 PM10/24/22
to sage-devel
more study of the bug (coming from maxima)

sage: C=x^2*(log(x))^4/((x+1)*(1+x^2))
sage: aa,bb=C.partial_fraction_decomposition()
sage: integral(aa,x,0,1)
-5/128*pi^5 + 45/64*zeta(5)
sage: integral(bb,x,0,1)
45/4*zeta(5)
sage: _+__
-5/128*pi^5 + 765/64*zeta(5)
sage: _.n()
0.440633136273039
sage: aa.nintegral(x,0,1)
(-11.58934902297507, 5.068708119893017e-08, 525, 0)
sage: bb.nintegral(x,0,1)
(11.66543724536065, 4.943314557692702e-08, 525, 0)
sage: integral(aa,x,0,1).n()
-11.2248041090899
sage: integral(bb,x,0,1).n()
11.6654372453629

Frédéric Chapoton

unread,
Oct 24, 2022, 3:24:30 PM10/24/22
to sage-devel
and one more step :

sage: integrate(x*log(x)^4/(x^2 + 1), x,0,1).n()
1.45817965567036
sage: (x*log(x)^4/(x^2 + 1)).nintegral(x,0,1)
(0.7290898278351722, 2.48288156701193e-09, 357, 0)
sage: integrate(-log(x)^4/(x^2 + 1), x,0,1).n()
-23.9077878738501
sage: (-log(x)^4/(x^2 + 1)).nintegral(x,0,1)
(-23.90778787384685, 1.267767046897461e-08, 483, 0)

Frédéric Chapoton

unread,
Oct 24, 2022, 3:33:41 PM10/24/22
to sage-devel
where we can see that there is a factor 2 between the wrong symbolic value and the correct numeric value

This should be filed as a bug in maxima.

Pablo Vitoria

unread,
Oct 24, 2022, 6:17:07 PM10/24/22
to sage-devel
Thank you very much for the analysis. If I understand correctly, the bug is in this part of the integral:  x*log(x)^4/(x^2 + 1)

I will try to report the bug to maxima. But, as far I can see https://sourceforge.net/p/maxima/bugs/ is very quiet, with few answers to the reports.

Emmanuel Charpentier

unread,
Oct 25, 2022, 9:46:13 AM10/25/22
to sage-devel

Essentially correct, bit it’s a tad less simpler :

sage: f(x) = x^2*(log(x))^4/((x+1)*(1+x^2)) 
sage: list(map(lambda u:u.integrate(x, 0, 1), f(x).partial_fraction_decomposition()))
[-5/128*pi^5 + 45/64*zeta(5), 45/4*zeta(5)]
sage: list(map(lambda u:mathematica.Integrate(u, (x, 0, 1)).sage(), f(x).partial_fraction_decomposition()))
[-5/128*pi^5 + 45/128*zeta(5), 45/4*zeta(5)]

HTH

Emmanuel Charpentier

unread,
Oct 25, 2022, 9:52:45 AM10/25/22
to sage-devel
Forgot : I have graphically checked the numerical quasi-equality of the symbolical integral given by Mathematica and the numerical integral given by Sage.

Frédéric Chapoton

unread,
Oct 30, 2022, 5:27:23 AM10/30/22
to sage-devel
The factor-2 problem is also present for the integral with no bounds. Please report in maxima's bug reporting site.

(%i20) integrate(x*log(x)^4/(x^2 + 1),x);
                     4         2
                2 log (x) log(x  + 1)           2                    2
(%o20) - (3 ((- ---------------------) + li (- x ) - 2 log(x) li (- x )
                          3                5                    4
                                                             3           2
                                                        4 log (x) li (- x )
                                       2           2                2
                                + 2 log (x) li (- x ) - -------------------))/2
                                              3                  3
(%i21) diff(%,x);  
                                         4
                                  2 x log (x)
(%o21)                            -----------
                                     2
                                    x  + 1

Pablo Vitoria

unread,
Oct 30, 2022, 6:49:31 PM10/30/22
to sage-devel
Ok, thanks a lot. I already reported the bug. 

Frédéric Chapoton

unread,
Oct 31, 2022, 2:53:22 AM10/31/22
to sage-devel
please post here the link to your report

Frédéric Chapoton

unread,
Oct 31, 2022, 5:44:12 AM10/31/22
to sage-devel
amusingly, one can see that some power of 2 is lurking around :

(%i36) expand(diff(integrate(x*log(x)**13/(1+x**2),x),x));                
                                          13
                                4096 x log  (x)
(%o36)                          ---------------
                                      2
                                  13 x  + 13

Le dimanche 30 octobre 2022 à 23:49:31 UTC+1, pvit...@gmail.com a écrit :

kcrisman

unread,
Nov 1, 2022, 7:01:33 AM11/1/22
to sage-devel
please post here the link to your report

Pablo Vitoria

unread,
Nov 1, 2022, 11:04:54 PM11/1/22
to sage-devel
I also added this to the Maxima bug report

Pablo Vitoria

unread,
Nov 1, 2022, 11:04:54 PM11/1/22
to sage-devel
Thanks a lot. I forgot to post it here

Frédéric Chapoton

unread,
Nov 10, 2022, 4:58:26 AM11/10/22
to sage-devel
note also that this works fine when the exponent is a variable ( I cannot post directly to sourceforge, sorry)


(%i2) expand (diff (integrate (x*log(x)^n/(1 + x^2), x), x));
                                        n
                                   x log (x)
(%o2)                              ---------
                                     2
                                    x  + 1
Reply all
Reply to author
Forward
0 new messages