bessel_i (2, %i*x) or bessel_i(2, x) with integrate --> noun form
- Subject: bessel_i (2, %i*x) or bessel_i(2, x) with integrate --> noun form
- From: Barton Willis
- Date: Sun, 23 Oct 2011 18:47:35 -0500
The conditional return is cute, but it will flummox diff, solve, and ...
Examples:
(%i62) integrate(bessel_i(a,x),x);
(%o62) if equal(a,-2) then (hypergeometric([3/2],[5/2,3],x^2/4)*x^3)/24 else (hypergeometric([(a+1)/2],[a+1,(a+3)/2],x^2/4)*x^(a+1))/(2^a*gamma(a+2))
(%i63) integrate(bessel_i(5,x),x);
(%o63) (hypergeometric([3],[4,6],x^2/4)*x^6)/23040
(%i64) integrate(bessel_i(-2,x),x);
(%o64) (hypergeometric([3/2],[5/2,3],x^2/4)*x^3)/24
What's the story with all this rat stuff....
(%i68) integrate(bessel_i(2,x),x,0,10.0);
rat: replaced 10.0 by 10/1 = 10.0
rat: replaced 10.0 by 10/1 = 10.0
rat: replaced 10.0 by 10/1 = 10.0
rat: replaced 2348.932064591061 by 241940/103 = 2348.932038834952
(%o68) 241940/103
The code passes a few tests such as
(%i94) integrate(bessel_i(-2,x),x);
(%o94) (hypergeometric([3/2],[5/2,3],x^2/4)*x^3)/24
(%i95) bessel_i(-2,x) - diff(%,x);
(%o95) -(hypergeometric([5/2],[7/2,4],x^2/4)*x^4)/240-(hypergeometric([3/2],[5/2,3],x^2/4)*x^2)/8+bessel_i(2,x)
(%i96) taylor(%,x,0,13);
(%o96)/T/ 0+...
;;----start---------------------------
(defun bessel-i-integral (a z)
(mfuncall 'mcond
(meqp a -2)
(mul (opcons 'mexpt z 3)
(div 1 24)
(opcons '$hypergeometric
(opcons 'mlist (div 3 2))
(opcons 'mlist (div 5 2) 3)
(div (mul z z) 4)))
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)
;;---end---------------------------------------------------
--Barton