isn't that funny? -Infinity is +Infinity (with double precision)

Discussion in 'Symbolic Algebra' started by Peter Pein, Mar 15, 2008.

  1. Peter Pein

    Peter Pein Guest

    Dear group,

    the following has been enterd into Maple 10. Is there a fix in later versions?
    sqrt(x) GAMMA(x + 1/2)
    f := x -> ----------------------
    GAMMA(x + 1)
    exp(-1) exp(1) exp(-1) exp(1)
    ------------------ - 1/8 --------------------
    exp(-1/2) exp(1/2) exp(-1/2) exp(1/2) x

    exp(-1) exp(1) 1
    + 1/128 --------------------- + O(----)
    2 3
    exp(-1/2) exp(1/2) x x
    1 1 1
    1 - --- + ------ + O(----)
    8 x 2 3
    128 x x

    the following integral should evaluate to a finite value (Mathematica gives
    0.00979687)
    infinity
    / 1/2
    | x GAMMA(x + 1/2) 1
    i0 := | ------------------- - 1 + --- dx
    | GAMMA(x + 1) 8 x
    /
    1
    Float(infinity)
    / 1/2
    | x GAMMA(x + 0.500000000000000000000000000000)
    | ------------------------------------------------
    | GAMMA(x + 1.)
    /
    1.

    0.125000000000000000000000000000
    - 1. + -------------------------------- dx
    x
    Classic Worksheet Interface, Maple 10.00, Windows, May 13 2005 B\
    uild ID 190196


    And why do I get 0.500... but not 1.000... in the evalf(i0,3*Digits)?

    I've got difficulties to understand the error messages which Maxima 5.14 gives me:
    (%i1) f(x):=sqrt(x)*gamma(x+1/2)/gamma(x+1)-1+1/8/x;
    (%o1) f(x):=(sqrt(x)*gamma(x+1/2))/gamma(x+1)-1+(1/8)/x
    (%i2) integrate(f(x),x,1,inf);
    `quotient' by `zero'
    -- an error. To debug this try debugmode(true);

    where can there be a division by zero?

    (%i3) quad_qagi(f(x),x,1,inf,1.e-8,1e6);
    Maxima encountered a Lisp error:
    Error in MACSYMA-TOP-LEVEL [or a callee]: 4000000.0 is not of type SEQUENCE.
    Automatically continuing.
    To reenable the Lisp debugger set *debugger-hook* to nil.

    Well, I know by myself, that 4000000.0 is not of type SEQUENCE. But how does
    it come to this error?

    And without specifications for accuracy and maximum steps I get the even more
    cryptic:

    (%i4) quad_qagi(f(x),x,1,inf);
    ***MESSAGE FROM ROUTINE DQAGI IN LIBRARY SLATEC.
    ***INFORMATIVE MESSAGE, PROG CONTINUES, TRACEBACK REQUESTED
    * ABNORMAL RETURN
    * ERROR NUMBER = 2
    *
    ***END OF MESSAGE

    Maxima encountered a Lisp error:

    Error in CONDITIONS::CLCS-UNIVERSAL-ERROR-HANDLER [or a callee]: Bind stack
    overflow.
    Fast links are on: do (use-fast-links nil) for debugging
    Broken at PCL::pRINT-STD-INSTANCE.
    1 (Continue) Maxima top-level
    2 (Abort) Return to debug level 1.
    3 Return to top level.
    dbl:MAXIMA>>

    What is going on here?
     
    Peter Pein, Mar 15, 2008
    #1
    1. Advertisements

  2. Dear Peter,

    Here come three Maple 11.02 defects for your convergent integral
    Int(x^(1/2)*GAMMA(x+.5000000000)/GAMMA(x+1.)-1.+.1250000000/x,
    x = 1. .. Float(infinity))
    Int(x^(1/2)*GAMMA(x+.50000000000000000000)/GAMMA(x+1.)-1.+
    .12500000000000000000/x,x = 1. .. Float(infinity))

    #
    # After a couple of minutes...
    #
    Float(undefined) # HOWLER!

    -------------------------------------------------------------

    Contrast with the instant perfect Mathematica 6.0.2 results

    NIntegrate[Sqrt[x] Gamma[x + 1/2]/Gamma[x + 1] - 1 + 1/(8 x),
    {x, 1, Infinity}]

    0.00979687

    NIntegrate[Sqrt[x] Gamma[x + 1/2]/Gamma[x + 1] - 1 + 1/(8 x),
    {x, 1, Infinity}, WorkingPrecision -> 30]

    0.00979686520191666594405859200437

    -------------------------------------------------------------

    PP> What is going on here?

    You did not upgrade to Mathematica 6.0.2 yet :)

    Cheers,

    Vladimir


     
    Vladimir Bondarenko, Mar 15, 2008
    #2
    1. Advertisements

  3. Peter Pein

    Peter Pein Guest

    And I do not want to while I'm reading about problems in the HelpBrowser and
    other problems.

    But at least the newer Maple does not imply that some finite value is -inf and
    calculated with double precision si +inf. It simply fails, which is OK to me,
    because no program of this complexity can be perfect.

    Cheerio!
    Peter
     
    Peter Pein, Mar 15, 2008
    #3
  4. Peter Pein

    Axel Vogt Guest

    This is when Maple finds, that the results has relative errors larger
    then requested (would be fine to have a mode, where it talks about it
    with suggestions ...)

    For your case best is to call with specific integration methods like

    Int(g(x),x=1..infinity, method = _NCrule);
    Int(g(x),x=1..infinity, method = _Dexp);

    and use evalf on it. The Newton Cotes method usually works in even
    very bad cases, the double exponential needs certain behavior in
    infinity (if I remember correctly).

    I always work with 14 Digits first, not 10 and get

    .97968652014535e-2

    This works in M10 and M11.
     
    Axel Vogt, Mar 15, 2008
    #4
  5. AV> .97968652014535e-2

    Could you please show the exact input in
    Maple 11/10 for which you got this number?

    Thanks in advance.
     
    Vladimir Bondarenko, Mar 15, 2008
    #5
  6. Peter Pein

    Axel Vogt Guest

    restart: interface(version); Digits:=14:

    Classic Worksheet Interface, Maple 11.02, Windows, Nov 10 2007 Build ID 330022


    f:= x -> sqrt(x)*GAMMA(x+1/2)/GAMMA(x+1):
    g:= x -> f(x)-1+1/8/x:

    Int(g(x),x=1..infinity, method = _NCrule): evalf(%);

    0.0097968652014535

    Int(g(x),x=1..infinity, method = _Dexp): evalf(%);

    0.0097968652014535
     
    Axel Vogt, Mar 15, 2008
    #6
  7. ......................................................

    Mathematica 6.0.2

    NIntegrate[Sqrt[x] Gamma[x + 1/2]/Gamma[x + 1] - 1 + 1/(8 x),
    {x, 1, Infinity}, WorkingPrecision -> 75]

    0.00979686520191666594405859200435067912845184025280088677877823557479\
    752620948

    ......................................................

    Maple 11.02

    restart; f:= x -> sqrt(x)*GAMMA(x+1/2)/GAMMA(x+1):g:= x -> f(x)-1+1/8/
    x:
    for j from 25 to 30 do
    print(j, evalf(Int(g(x),x=1..infinity, method = _Dexp,j))); od;


    25, Int(x^(1/2)*GAMMA(x+.5000000000000000000000000)/GAMMA(x+1.)-1.+.
    1250000000000000000000000/x,x = 1. .. Float(infinity))

    26, Int(x^(1/2)*GAMMA(x+.50000000000000000000000000)/GAMMA(x+1.)-1.+.
    12500000000000000000000000/x,x = 1. .. Float(infinity))

    27, Int(x^(1/2)*GAMMA(x+.500000000000000000000000000)/GAMMA(x+1.)-1.+.
    125000000000000000000000000/x,x = 1. .. Float(infinity))

    28, Int(x^(1/2)*GAMMA(x+.5000000000000000000000000000)/GAMMA(x+1.)-1.+.
    1250000000000000000000000000/x,x = 1. .. Float(infinity))

    29, Int(x^(1/2)*GAMMA(x+.50000000000000000000000000000)/GAMMA(x
    +1.)-1.+.12500000000000000000000000000/x,x = 1. .. Float(infinity))

    30, Int(x^(1/2)*GAMMA(x+.500000000000000000000000000000)/GAMMA(x
    +1.)-1.+.125000000000000000000000000000/x,x = 1. .. Float(infinity))

    restart; f:= x -> sqrt(x)*GAMMA(x+1/2)/GAMMA(x+1):g:= x -> f(x)-1+1/8/
    x:
    for j from 25 to 30 do
    print(j, evalf(Int(g(x),x=1..infinity, method = _NCrule,j))); od;

    25, Int(x^(1/2)*GAMMA(x+.5000000000000000000000000)/GAMMA(x+1.)-1.+.
    1250000000000000000000000/x,x = 1. .. Float(infinity))

    26, Int(x^(1/2)*GAMMA(x+.50000000000000000000000000)/GAMMA(x+1.)-1.+.
    12500000000000000000000000/x,x = 1. .. Float(infinity))

    27, Int(x^(1/2)*GAMMA(x+.500000000000000000000000000)/GAMMA(x+1.)-1.+.
    125000000000000000000000000/x,x = 1. .. Float(infinity))

    28, Int(x^(1/2)*GAMMA(x+.5000000000000000000000000000)/GAMMA(x+1.)-1.+.
    1250000000000000000000000000/x,x = 1. .. Float(infinity))

    29, Int(x^(1/2)*GAMMA(x+.50000000000000000000000000000)/GAMMA(x
    +1.)-1.+.12500000000000000000000000000/x,x = 1. .. Float(infinity))

    30, Int(x^(1/2)*GAMMA(x+.500000000000000000000000000000)/GAMMA(x
    +1.)-1.+.125000000000000000000000000000/x,x = 1. .. Float(infinity))
     
    Vladimir Bondarenko, Mar 15, 2008
    #7
  8. Peter Pein

    Axel Vogt Guest

    ....

    As quite often you are a bit noisy. Instead try to read what I have
    written before and you will find it out yourself.
     
    Axel Vogt, Mar 15, 2008
    #8
    1. Advertisements

Ask a Question

Want to reply to this thread or ask your own question?

You'll need to choose a username for the site, which only take a couple of moments (here). After that, you can post your question and our members will help you out.