failing calculation of an integral

[casamayou]

In the following calculation, Sage4.6.2 returns the opposite of the right result :

sage: var('a t')
(a, t)
sage: assume(a>0)
sage: assume(x>0)
sage: f = integrate(log(1+a/(x * t)^2), t, 1, oo)
sage: f
(sqrt(a)*x*log((x^2 + a)/x^2) - 2*a*arctan(sqrt(a)/x))/(sqrt(a)*x)
sage: f.subs(x=1, a=7).n()
-4.32025625668262

For information, Sage 4.6 gives the right result :

sage: var('a t')
(a, t)
sage: assume(a>0)
sage: assume(x>0)
sage: f = integrate(log(1+a/(x * t)^2), t, 1, oo)
sage: f
pi*sqrt(a)/x - (x*log((x2 + a)/x2) + 2*sqrt(a)*arctan(x/sqrt(a)))/x
sage: f.subs(x=1, a=7).n()
4.32025625668262

For information, Maple9 gives this :

> assume(a>0): assume(x>0):
> f := int(ln(1+a/(x * t)^2), t=1..infinity):
> lprint(f);
(2*ln(x)*x-2*a^(1/2)*arctan(x/a^(1/2))-ln(x^2+a)*x+a^(1/2)*Pi)/x
> evalf(subs(x=1, a=7, f));
bytes used=4000512, alloc=3341724, time=0.13
                                4.320256257

[kcrisman]

I can confirm this in Maxima 5.24.0:

Maxima 5.24.0 http://maxima.sourceforge.net
using Lisp SBCL 1.0.24
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) assume(a>0);
(%o1)                               [a > 0]
(%i2) assume(x>0);
(%o2)                               [x > 0]
(%i3) integrate(log(1+a/(x*t)^2),t,1,inf); 
                                  2
                                 x  + a             sqrt(a)
                   sqrt(a) x log(------) - 2 a atan(-------)
                                    2                  x
                                   x
(%o3)              -----------------------------------------
                                   sqrt(a) x

Looks like a different branch of log is being used. Also for comparison, here is what GSL says:

sage: numerical_integral(log(1+7/(x^2)),1,oo)
(4.3202562566855764, 2.2384237933839469e-07)

We should file something with the Maxima tracker.

[kcrisman]

This is now ​this bug report at the Maxima tracker.

[kcrisman]

Apparently this is fixed upstream now.

[jen]

Checking against 5.0.beta9, this looks good:

sage: var('a t')
(a, t)
sage: assume(a>0)
sage: assume(x>0)
sage: f = integrate(log(1+a/(x * t)^2), t, 1, oo)
sage: f
(sqrt(a)*x*log((x^2 + a)/x^2) - 2*a*arctan(sqrt(a)/x))/(sqrt(a)*x)
sage: f.subs(x=1, a=7).n()
4.32025625668262

[kcrisman]

Okay, so then all we need is a doctest referring to this ticket.

[mjo]

The scientific notation change in the abs_tol above is unrelated, but I'm 99% sure I did that and it's a little embarrassing =)

[zimmerma]

I give a positive review. Thank you Michael!

Paul

[jdemeyer]

The patch needs to be rebased to sage-5.0.beta14:

applying sage-trac_11233.patch
patching file sage/symbolic/integration/integral.py
Hunk #1 FAILED at 580
1 out of 1 hunks FAILED -- saving rejects to file sage/symbolic/integration/integral.py.rej
patch failed, unable to continue (try -v)
patch failed, rejects left in working dir
errors during apply, please fix and refresh sage-trac_11233.patch

[zimmerma]

Michael, please can you rebase your patch?

Paul

[mjo]

Rebased patch on 5.0.beta14

[mjo]

Sure, done.

[zimmerma]

I only checked that the rebased patch applies cleanly to sage 5.0.beta10, and that the integral returns the correct result, assuming the rest did not change.

Paul

[zimmerma]

oops, that should be "sage.5.0.beta14" above.

Paul