Comparison of *integrate* with bessel functions using integer and float orders

comparison of integrate with bessel functions using integer and float orders
---------------------------
SUMMARY:

   Using new bessel.lisp and Barton Willis' new
  besselI.lisp (see below), both bessel_j and bessel_i
  yield a good answer for either integral or float
  order.

  Both bessel_k and bessel_y produce a noun
  form when a float order is used.
--------------------------------------------

DETAILS:
=========================
case bessel_j (n,z) as arg to integrate:
   produces good answers for either integral or float order
=========================

  for z = x = real,

(%i2) integrate(bessel_j(0,x),x);
(%o2) 
(bessel_j(0,x)*(2-%pi*struve_h(1,x))+%pi*struve_h(0,x)*bessel_j(1,x))*x
       /2
(%i3) integrate(bessel_j(0.0,x),x);
(%o3) 1.0*hypergeometric([0.5],[1.0,1.5],-x^2/4)*x
--------------------------------
   for z = %i*x = imaginary

(%i15) integrate(bessel_j(0,%i*x),x);
(%o15) (bessel_j(0,%i*x)*(2-%pi*struve_h(1,%i*x))
        +%pi*struve_h(0,%i*x)*bessel_j(1,%i*x))
        *x
        /2

(%i68) integrate(bessel_j(0.0,%i*x),x);
(%o68) 1.0*hypergeometric([0.5],[1.0,1.5],x^2/4)*x

======================================
case bessel_i (n,z) : get answer with either integral or float order
===============================
(MUST use Barton Willis' new besselI.lisp)

  for z = x:

(%i1) (display2d:false,ratprint:false,ratepsilon:2.0e-15)$

(%i2) integrate(bessel_i(4,x),x);
(%o2) 'integrate(bessel_i(4,x),x)

(%i3) load("bessel.lisp");
(%o3) "bessel.lisp"

(%i4) integrate(bessel_i(4,x),x);
(%o4) 'integrate(bessel_i(4,x),x)

(%i5) load("besselI.lisp");
(%o5) "besselI.lisp"

(%i6) integrate(bessel_i(4,x),x);
(%o6) hypergeometric([5/2],[7/2,5],x^2/4)*x^5/1920

(%i7) integrate(bessel_i(4.0,x),x);
(%o7) 5.2083333333333333E-4*hypergeometric([2.5],[3.5,5.0],x^2/4)*x^5.0

(%i25) integrate(bessel_i(0,x),x);
(%o25) hypergeometric([1/2],[1,3/2],x^2/4)*x
(%i26) integrate(bessel_i(0.0,x),x);
(%o26) 1.0*hypergeometric([0.5],[1.0,1.5],x^2/4)*x

  for z = %i*x:

(%i46) integrate(bessel_i(4,%i*x),x);
(%o46) hypergeometric([5/2],[7/2,5],-x^2/4)*x^5/1920
(%i47) integrate(bessel_i(4.0,%i*x),x);
(%o47) -5.2083333333333333E-4*%i*hypergeometric([2.5],[3.5,5.0],-x^2/4)
                             *(%i*x)^5.0

======================================
case bessel_k (n,z): only integral order gets answer
===============================

  for z = x  only integral order is done

(%i14) integrate(bessel_k(0,x),x);
(%o14) %pi*(struve_l(0,x)*bessel_k(1,x)+struve_l(-1,x)*bessel_k(0,x))*x/2

(%i15) integrate(bessel_k(0.0,x),x);
(%o15) 'integrate(bessel_k(0.0,x),x)
--------------------------------------
  for z = %i*x  only integral order is done

(%i9) integrate(bessel_k(0,%i*x),x);
(%o9) %pi*(struve_l(0,%i*x)*bessel_k(1,%i*x)
          +struve_l(-1,%i*x)*bessel_k(0,%i*x))*x
       /2

(%i10) integrate(bessel_k(0.0,%i*x),x);
(%o10) 'integrate(bessel_k(0.0,%i*x),x)

 =====================================
 case bessel_y (n,z):  only integral order gets answer
===============================

  for z = x

(%i20) integrate(bessel_y(0,x),x);
(%o20) %pi*(struve_h(0,x)*bessel_y(1,x)+struve_h(-1,x)*bessel_y(0,x))*x/2

(%i21) integrate(bessel_y(0.0,x),x);
(%o21) 'integrate(bessel_y(0.0,x),x)


--------------------------------------
for z = %i*x

(%i15) integrate(bessel_y(0,%i*x),x);
(%o15) %pi*(struve_h(0,%i*x)*bessel_y(1,%i*x)
           +struve_h(-1,%i*x)*bessel_y(0,%i*x))*x
        /2

(%i16) integrate(bessel_y(0.0,%i*x),x);
(%o16) 'integrate(bessel_y(0.0,%i*x),x)

--------------------------------------
 Ted Woollett

--------------------------------
p.s.

Barton Willis' code besselI.lisp
-----------------------------------
;;; file besselI.lisp

;;;code from barton willis, oct. 24, 2011
 ;; integrate(bessel_i(a,z),z) = 
hypergeometric([(a+1)/2],[a+1,(a+3)/2],x^2/4)*x^(a+1))/(2^a*gamma(a+2))
 ;; Maxima simplifies bessel_i(-2,z) --> bessel_(2,z), so a = -2 isn't a 
special case.

(defun bessel-i-integral (a z)
  "Integral of bessel_i(a,z) wrt z"
  (cond ((eq t (meqp a 1)) ;; integrate(bessel_i(1,z),z) = bessel_i(0,z)
  (opcons '%bessel_i 0 z))
 (t
  (mul
   (opcons 'mexpt z (add a 1))
   (opcons '$hypergeometric
    (opcons 'mlist (div (add a 1) 2))
    (opcons 'mlist (div (add a 3) 2) (add a 1))
    (div (mul z z) 4))
   (div
    1
    (mul
     (opcons 'mexpt 2 a)
     (opcons '%gamma (add a 2))))))))

(putprop '%bessel_i `((a z) nil ,#'bessel-i-integral) 'integral)

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