Learn how easy it is to sync an existing GitHub or Google Code repo to a SourceForge project! See Demo

Close

#1405 bessel_i(1/2,0) -> divide by zero error

closed
nobody
Lisp Core (472)
5
2008-07-09
2008-04-30
Raymond Toy
No

bessel_i(1/2,0) produces an error. I think it should be 0 since bessel_i(1/2,x) is sqrt(2/%pi)*sinh(x)/sqrt(x).

Discussion

  • Raymond Toy
    Raymond Toy
    2008-06-20

    Logged In: YES
    user_id=28849
    Originator: YES

    Doesn't seem get a divide by zero error anymore in latest CVS. It returns the noun form.

    But bessel_i(v,x) ~ (1/z*x)^v/gamma(v+1) for small x, it should be 0.

    I think the bug is in simp-bessel-i in the lines:

           \(\(or \(and \(numberp order\) \(> order 0\)\) \(integerp order\)\)
            0\)
           \(\(and \(numberp order\) \(< order 0\)\)
            '$infinity\)
    

    Perhaps numberp should be changed to mnump and (> order 0) should be (great order 0)

     
  • Raymond Toy
    Raymond Toy
    2008-06-20

    Logged In: YES
    user_id=28849
    Originator: YES

    Doesn't seem get a divide by zero error anymore in latest CVS. It returns the noun form.

    But bessel_i(v,x) ~ (1/z*x)^v/gamma(v+1) for small x, it should be 0.

    I think the bug is in simp-bessel-i in the lines:

           \(\(or \(and \(numberp order\) \(> order 0\)\) \(integerp order\)\)
            0\)
           \(\(and \(numberp order\) \(< order 0\)\)
            '$infinity\)
    

    Perhaps numberp should be changed to mnump and (> order 0) should be (great order 0)

     
  • Raymond Toy
    Raymond Toy
    2008-06-20

    Logged In: YES
    user_id=28849
    Originator: YES

    Doesn't seem get a divide by zero error anymore in latest CVS. It returns the noun form.

    But bessel_i(v,x) ~ (1/z*x)^v/gamma(v+1) for small x, it should be 0.

    I think the bug is in simp-bessel-i in the lines:

           \(\(or \(and \(numberp order\) \(> order 0\)\) \(integerp order\)\)
            0\)
           \(\(and \(numberp order\) \(< order 0\)\)
            '$infinity\)
    

    Perhaps numberp should be changed to mnump and (> order 0) should be (great order 0)

     
  • Dieter Kaiser
    Dieter Kaiser
    2008-07-04

    Logged In: YES
    user_id=2039760
    Originator: NO

    If the gloabal flag BESSELEXPAND is FALSE we get for all Bessel functions with a Rational order and a zero as argument a noun form with the CVS code 1.58. (In previous versions a lot of the zero checks were not implemented.)

    But, if we set the global flag BESSELEXPAND to TRUE the Bessel function will be expanded and then evaluated for the argument zero. We get for all Bessel functions the described error.

    I think it is a question of convention and the desired behaviour of Maxima what we should do.

    The code to catch the zero argument is only implemented for CL numbers. The last time the code was extended to look more carefully for zero arguments didn't extend this to Maxima Rational, Bigfloat numbers or even expression which might simplify.

    In my opinion the reason to look only for CL numbers is, that only for such numbers we do the numerical evaluation. We check for a zero argument to prevent Lisp Errors in the numerical slatec routines. At this point we need to do the simplification for the Bessel function with zero argument.

    In the case of Rational numbers it depends on the flag BESSELEXPAND what we should do. In my opinion it is the best to check for the zero argument before we actually do the expansion.

    If we try to simplify all Bessel functions with zero argument in advance we get a lot of problems and questions, too. The reason is that the result of the simplification with an argument of zero depends on the sign of the order and is different for complex and real order of the Bessel function.

    Here an example what we can do in simp-bessel-i:

    ((and $besselexpand (setq rat-order (max-numeric-ratio-p order 2)))
    (cond
    ((= arg 0)
    ;; We don't do the expansion for a zero argument.
    (if (> rat-order 0) 0 '$infinity))
    (t
    ;; When order is a fraction with a denominator of 2, we
    ;; can express the result in terms of elementary
    ;; functions.
    ;;
    ;; I[1/2](z) = sqrt(2/%pi/z)*sinh(z)
    ;; I[-1/2](z) = sqrt(2/%pi/z)*cosh(z)
    (bessel-i-half-order rat-order arg))))

    We no longer get the divison by zero error. Here the results:

    (%i7) besselexpand:false;
    (%o7) false

    (%i8) bessel_i(1/2,0);
    (%o8) bessel_i(1/2,0)

    (%i9) bessel_i(-1/2,0);
    (%o9) bessel_i(-1/2,0)

    (%i10) besselexpand:true;
    (%o10) true

    (%i11) bessel_i(1/2,0);
    (%o11) 0

    (%i12) bessel_i(-1/2,0);
    (%o12) infinity

    This check for a zero argument before we expand the Bessel function could be implemented for all other Bessel functions too.

    Dieter Kaiser

     
  • Dieter Kaiser
    Dieter Kaiser
    2008-07-04

    Logged In: YES
    user_id=2039760
    Originator: NO

    Sorry, but I have omitted to check arg to be an CL number. The correct code would be:

    ((and $besselexpand (setq rat-order (max-numeric-ratio-p order 2)))
    (cond
    ((and (numberp arg) (= arg 0))
    (if (> rat-order 0) 0 '$infinity))
    ( t
    ;; When order is a fraction with a denominator of 2, we
    ;; can express the result in terms of elementary
    ;; functions.
    ;;
    ;; I[1/2](z) = sqrt(2/%pi/z)*sinh(z)
    ;; I[-1/2](z) = sqrt(2/%pi/z)*cosh(z)
    (bessel-i-half-order rat-order arg))))

    Dieter Kaiser

     
  • Dieter Kaiser
    Dieter Kaiser
    2008-07-09

    Logged In: YES
    user_id=2039760
    Originator: NO

    For all Bessel functions a test has been added, which check for a zero argument before the expansion in terms of elementary functions.

    Example for Bessel I:

    besselexpand:false;
    bessel_i(1/2,0) -> bessel_i(1/2,0)

    besselexpand:true;
    bessel_i(1/2,0) -> 0

    Close the bug.

     
  • Dieter Kaiser
    Dieter Kaiser
    2008-07-09

    • status: open --> closed