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

closed
nobody
5
2009-09-23
2009-09-22
No

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

## Discussion

• 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 - 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 - 2009-09-23
• status: open --> closed