#3894 integrating besssel_k expressions

[Barton Willis]

There is almost surely a missing comma in front of z in bessel-k-integral-2. I don't know how to reliably make this error manifest, but sometimes with my own build, I get this bug after running the testsuite.

Also there is a call to $ev. Oh, my. Maybe it would be better to write code like this using the Maxima reader macro?

(defun bessel-k-integral-2 (n z)
  (cond 
   ((and ($integerp n) (<= 0 n))
    (cond
     (($oddp n)
      ;; integrate(bessel_k(2*N+1,z),z) , N > 0
      ;; = -(-1)^((n-1)/2)*bessel_k(0,z) 
      ;;   + 2*sum((-1)^(k+(n-1)/2-1)*bessel_k(2*k,z),k,1,(n-1)/2)
      (let* ((k (gensym))
         (answer `((mplus)
               ((mtimes) -1 ((%bessel_k) 0 ,z)
            ((mexpt) -1
             ((mtimes) ((rat) 1 2) ((mplus) -1 ,n))))
               ((mtimes) 2
            ((%sum)
             ((mtimes) ((%bessel_k) ((mtimes) 2 ,k) ,z)
              ((mexpt) -1
               ((mplus) -1 ,k
                ((mtimes) ((rat) 1 2) ((mplus) -1 ,n)))))
             ,k 1 ((mtimes) ((rat) 1 2) ((mplus) -1 ,n)))))))
    ;; Expand out the sum if n < 10.  Otherwise fix up the indices
    (if (< n 10) 
            (meval `(($ev) ,answer $sum))   ; Is there a better way?
      (simplify ($niceindices answer)))))
     (($evenp n)
      ;; integrate(bessel_k(2*N,z),z) , N > 0
      ;; = (1/2)*(-1)^(n/2)*%pi*z*(bessel_k(0,z)*struve_l(-1,z)
      ;;               +bessel_k(1,z)*struve_l(0,z))
      ;;    + 2 * sum((-1)^(k+n/2)*bessel_k(2*k+1,z),k,0,n/2-1)
      (let* 
      ((k (gensym))
       (answer `((mplus)
             ((mtimes) 2
              ((%sum)
               ((mtimes)
            ((%bessel_k) ((mplus) 1 ((mtimes) 2 ,k)) ,z)
            ((mexpt) -1
             ((mplus) ,k ((mtimes) ((rat) 1 2) ,n))))
               ,k 0 ((mplus) -1 ((mtimes) ((rat) 1 2) ,n))))
                 ((mtimes)
                  ((rat) 1 2)
                  $%pi
                  ((mexpt) -1 ((mtimes) ((rat) 1 2) ,n))
                  ,z ;<===========================was z
              ((mplus)
               ((mtimes)
                ((%bessel_k) 0 ,z)
            ((%struve_l) -1 ,z))
               ((mtimes)
                ((%bessel_k) 1 ,z)
            ((%struve_l) 0 ,z)))))))
    ;; expand out the sum if n < 10.  Otherwise fix up the indices
    (if (< n 10) 
            (meval `(($ev) ,answer $sum))  ; Is there a better way?
      (simplify ($niceindices answer)))))))
   (t nil)))

[Barton Willis]

I think there is another bug--with my own build and after running the testsuite & share testsuite, I get

    integrate(bessel_k(1,x),x);
ev: improper argument: sum

[Stavros Macrakis]

Yes, ev in Lisp code is always bad news. Unfortunately, adding noeval prevents the sum from having an effect. We need to build better functions for controlled re-evaluation. I have some Maxima code for this (as a prototype) which I should rewrite in Lisp.

[Barton Willis]

For the record, a bug caused by ev:

(%i1)   x : y;
(x) y
(%i2)   y : z;
(y) z
(%i3)   integrate(bessel_k(2,x),x);
(%o3)   -(%pi (struve_l(0,z)*bessel_k(1,z)+struve_l(-1,z)*bessel_k(0,z))*y)/2-2*bessel_k(1,z)

The result has both variables y & z.