#1772 specint(exp(-s*t)*bessel_i(1/2,t)^2,t) not correct


I think we have some problems with the Laplace transform of Bessel functions with a half integral index.

This is an example for the square of the bessel_i(1/2,t) function:

(%i1) assume(s>0)$

We do the Laplace transform for the square of the expanded bessel_i(1/2,t) function and use the algorithm for the sinh function to get the Laplace transform:

(%i2) expr1:bessel_i(1/2,t)^2,besselexpand:true;
(%o2) 2*sinh(t)^2/(%pi*t)

(%i3) specint(exp(-s*t)*expr1,t);
(%o3) -log(1-4/s^2)/(2*%pi)

Aagain, but we do not expand the function. Now, the hypergeometric algorithm is used.

(%i4) expr2:bessel_i(1/2,t)^2;
(%o4) bessel_i(1/2,t)^2
(%i5) specint(exp(-s*t)*expr2,t);
(%o5) (%i-1)*%i*log(1-4/s^2)/(2^(3/2)*%pi)

The answers differ by a factor sqrt(%i). First the ratio of the two answers:

(%i17) specint(exp(-s*t)*expr1,t)/specint(exp(-s*t)*expr2,t);
(%o17) sqrt(2)*%i/(%i-1)

The ratio differs by a factor sqrt(%i):

(%i18) rectform(sqrt(%i)*%);
(%o18) 1

I think we have more of such problems with the Laplace transform of Bessel functions.

There is one expected failure for the Laplace transform of bessel_y(1/2,sqrt(t))^2 in the testfile rtest14.mac which might be related to this problem. I had already a long search to find the bug, but had no success. Perhaps it is a mathematical problem related to the hypergeometric transformation and integration for a half integral index.

Dieter Kaiser


  • Dieter Kaiser

    Dieter Kaiser - 2009-09-22

    I have found the bug for the square of the Bessel I function. We transform to two Bessel J functions. In the transformation is missing a factor %i^-v. (In the following code I have already replaced the function 1fact and have inserted the powers of %i).

    ;; Laplace transform of square of Bessel I function
    (cond ((setq l (onei^2 u))
    (setq index1 (cdras 'v l)
    arg1 (mul '$%i (cdras 'w l))
    rest (mul (power '$%i (neg index1))
    (power '$%i (neg index1)) ; the missing factor
    (cdras 'u l)))
    (return (lt1j^2 rest arg1 index1))))

    Now we get:

    (%i3) assume(s>0)$

    We expand the Bessel function:

    (%i4) specint(exp(-s*t)*bessel_i(1/2,t)^2,t),besselexpand:true;
    (%o4) -log(1-4/s^2)/(2*%pi)

    Now we use the hypergeometric code and get the same result:

    (%i5) specint(exp(-s*t)*bessel_i(1/2,t)^2,t);
    (%o5) -log(1-4/s^2)/(2*%pi)

    Dieter Kaiser

  • Dieter Kaiser

    Dieter Kaiser - 2009-09-23

    The problem for bessel_i(1/2,t)^2 is fixed in hypgeo.lisp revsion 1.62.
    The problem for bessel_y(1/2,sqrt(t)) is still open. But it is another bug.
    Closing this bug report.

    Dieter Kaiser

  • Dieter Kaiser

    Dieter Kaiser - 2009-09-23
    • status: open --> closed

Log in to post a comment.

Get latest updates about Open Source Projects, Conferences and News.

Sign up for the SourceForge newsletter:

JavaScript is required for this form.

No, thanks