specint(exp(-s*t)*t^(5/2)*bessel_j(-1/2,sqrt(t))^2,t)



We have the following known failure in rtest14.mac:

********************** Problem 57 ***************
Input:
specint(t^(5/2)*bessel_y(1/2,t^(1/2))^2*%e^-(p*t),t)


In $specint the expression with the bessel_y function is transformed to
the square of the bessel_j function. So we get the following integrand:

(1)   t^(5/2)*(bessel_j(-1/2,sqrt(t))^2

Furthermore, this expression is equivalent to:

(2)   2/%pi*cos(sqrt(t))^2

Maxima can do this transformation:

(%i17) t^(5/2)*bessel_j(-1/2,sqrt(t))^2,besselexpand:true;
(%o17) 2*cos(sqrt(t))^2*t^2/%pi

The problem is that the integrands (1) and (2) give different Laplace
transforms:

First the result for bessel_j(-1/2,sqrt(t))^2:

(%i14)
res1:factor(ratsimp(specint(exp(-s*t)*t^(5/2)*bessel_j(-1/2,sqrt(t))^2,t)));
(%o14) %e^-(1/s)*(8*s^3*%e^(1/s)-18*s^2*%e^(1/s)+4*s*%e^(1/s)
                                +15*sqrt(%pi)*%i*erf(%i/sqrt(s))*s^(5/2)
                                -20*sqrt(%pi)*%i*erf(%i/sqrt(s))*s^(3/2)
                                +4*sqrt(%pi)*%i*erf(%i/sqrt(s))*sqrt(s))
        /(2*%pi*s^6)

Next, the result for cos(sqrt(t))^2 (we use the flag besselexpand):

(%i15) res2 : factor(ratsimp(specint(exp(-s*t) * t^(5/2) *
bessel_j(-1/2, sqrt(t))^2, t))), besselexpand:true;
(%o15) %e^-(1/s)*(16*s^3*%e^(1/s)-18*s^2*%e^(1/s)+4*s*%e^(1/s)
                         +15*sqrt(%pi)*%i*erf(%i/sqrt(s))*s^(5/2)
                         -20*sqrt(%pi)*%i*erf(%i/sqrt(s))*s^(3/2)
                         +4*sqrt(%pi)*%i*erf(%i/sqrt(s))*sqrt(s))
        /(4*%pi*s^6)

The results differ by a factor 2 in most, but not in all terms.

I had a long search for the bug and I have found the problem in the
algorithm for the product of hypergeometric functions.

Maxima does the following transformation for our case of two
bessel_j(-1/2,sqrt(t)) functions:

bessel_j(-1/2,sqrt(t))^2 --> 2/%pi*2F3([0,1/2], [1/2,1/2,0], -t)

Next the hypergeometric function is reduced in two steps:

   2F3([0,1/2], [1/2,1/2,0], -t) --> 1F2([1/2], [1/2,1/2], -t) 
                                 --> 0F1([], [1/2], -t)

But, 0F1([],[1/2],-t) represents cos(2*sqrt(t)) and not cos(sqrt(t))^2
as expected. Therefore, we get the Laplace transform of cos(2*sqrt(t))
and not of cos(sqrt(t))^2, when we use the hypergeometric algorithm. We
can check this by doing the Laplace transform of cos(2*sqrt(t))
directly.

A correct hypergeometric representation of cos(sqrt(t))^2 is

   1/2*(0F1([],[1/2],-t) + 1).

The error is, that we do the following transformation for a parameter
a=0:

   2F3([a,1/2],[1/2,1/2,a],-t)  --> 1F2([1/2], [1/2,1/2], -t).

I think this transformation is not valid for a=0, because the
hypergeometric function 2F3 is not well defined for this case.

The parameter a is zero for v+u=-1, where v and u are the order of the
two Bessel J functions involved in the Laplace transformation.

For the square of bessel functions this is the case for v=-1/2.

Remark: There are a lot of more possibilities for wrong results. The
example in rtest14.mac is not the simplest one.

Dieter Kaiser