Opened 8 years ago

Closed 7 years ago

#11233 closed defect (fixed)

failing calculation of an integral

Reported by: casamayou Owned by: burcin
Priority: major Milestone: sage-5.1
Component: calculus Keywords: integrate
Cc: zimmerma Merged in: sage-5.1.beta0
Authors: Michael Orlitzky Reviewers: Paul Zimmermann
Report Upstream: Fixed upstream, in a later stable release. Work issues:
Branch: Commit:
Dependencies: Stopgaps:

Description

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

Attachments (1)

sage-trac_11233.patch (1.4 KB) - added by mjo 7 years ago.
Rebased patch on 5.0.beta14

Download all attachments as: .zip

Change History (16)

comment:1 Changed 8 years ago by kcrisman

  • Report Upstream changed from N/A to Not yet reported upstream; Will do shortly.

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.

comment:2 Changed 8 years ago by kcrisman

  • Report Upstream changed from Not yet reported upstream; Will do shortly. to Reported upstream. Little or no feedback.

This is now this bug report at the Maxima tracker.

comment:3 Changed 8 years ago by kcrisman

  • Report Upstream changed from Reported upstream. Little or no feedback. to Fixed upstream, but not in a stable release.

Apparently this is fixed upstream now.

comment:4 Changed 7 years ago by 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

comment:5 Changed 7 years ago by kcrisman

  • Report Upstream changed from Fixed upstream, but not in a stable release. to Fixed upstream, in a later stable release.

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

comment:6 Changed 7 years ago by mjo

  • Authors set to Michael Orlitzky
  • Status changed from new to needs_review

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 =)

comment:7 Changed 7 years ago by zimmerma

  • Reviewers set to Paul Zimmermann
  • Status changed from needs_review to positive_review

I give a positive review. Thank you Michael!

Paul

comment:8 Changed 7 years ago by jdemeyer

  • Milestone changed from sage-5.0 to sage-5.1

comment:9 Changed 7 years ago by jdemeyer

  • Status changed from positive_review to needs_work

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

comment:10 Changed 7 years ago by zimmerma

Michael, please can you rebase your patch?

Paul

Changed 7 years ago by mjo

Rebased patch on 5.0.beta14

comment:11 Changed 7 years ago by mjo

Sure, done.

comment:12 Changed 7 years ago by mjo

  • Status changed from needs_work to needs_review

comment:13 Changed 7 years ago by zimmerma

  • Status changed from needs_review to positive_review

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

comment:14 Changed 7 years ago by zimmerma

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

Paul

comment:15 Changed 7 years ago by jdemeyer

  • Merged in set to sage-5.1.beta0
  • Resolution set to fixed
  • Status changed from positive_review to closed
Note: See TracTickets for help on using tickets.