bessel_i (2, %i*x) or bessel_i(2, x) with integrate --> noun form

On Oct. 24, 6:53 CDT, 2011, Barton Willis wrote:
---------------------------------------
>Oops...Maxima simplifies bessel_i(-2,x) --> bessel_i(2,x). So a = -2, isn't 
>a special case
>for integrate(bessel_i(a,x),x). The special case
>
> (%i7) integrate(bessel_i(0,x),x);
> (%o7) 
> ((bessel_i(0,x)*(%pi*struve_l(1,x)+2)-%pi*struve_l(0,x)*bessel_i(1,x))*x)/2
>
>is more complicated than using the hypergeometric representation, I think. 
>Thus let's try
>
>(defun bessel-i-integral (a z)
> [etc]
-------------------------------
This code works for my numerical needs:

1. integration over a large range:

(%i1) load(nint);
(%o1) "c:/work2/nint.mac"

(%i9) f(integrate(bessel_i(2,%i*x),x,1,500));
(%o9) -0.94996630633619
(%i10) f(integrate(bessel_i(3,%i*x),x,1,500));
(%o10) -0.96082031863327*%i
(%i11) f(integrate(bessel_i(4,%i*x),x,1,500));
(%o11) 1.009491962084913
(%i12) f(integrate(bessel_i(5,%i*x),x,1,500));
(%o12) 1.03430327860412*%i

replacing 500 >> 1000 leads to
Exceeded maximum allowed fpprec
error message.





2.  for numerical integration, ok to use bessel_i(2,x)
or bessel_i(2.0,x) which allows nint to avoid a
present quadpack trap with imagpart(bessel_i).

(%i16) f(integrate(bessel_i(0,%i*x),x,1,3));
(%o16) 0.4678368419201
(%i17) f(integrate(bessel_i(0.0,%i*x),x,1,3));
(%o17) 0.46783686924217

(%i18) f(integrate(bessel_i(1,%i*x),x,1,3));
(%o18) 1.0252496414599*%i
(%i19) f(integrate(bessel_i(1.0,%i*x),x,1,3));
(%o19) 1.0252496414599*%i

(%i20) f(integrate(bessel_i(2,%i*x),x,1,3));
(%o20) -0.6698200963581
(%i21) f(integrate(bessel_i(2.0,%i*x),x,1,3));
(%o21) -0.6698200859962

(%i22) f(integrate(bessel_i(3,%i*x),x,1,3));
(%o22) -0.28287409015192*%i
(%i23) f(integrate(bessel_i(3.0,%i*x),x,1,3));
(%o23) -0.28287408738349*%i

nint.mac loads nint.lisp which at present looks like:
-----------------------------------------------------------
;;;nint.lisp

;;;code from barton willis file
;;; topoly.lisp

(defun $complex_number_p (e)
  (complex-number-p e #'$numberp))


;;;code from barton oct. 24,2011, 6:53 CDT
 ;; 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)
 -----------------------------------------------------------------

Ted