sdiffgrad - two small extensions



We have a bug report "derivative of gamma_incomplete - ID: 2793294". The
problem is, that the implemented derivative is only valid for a
parameter a>0 for gamma_incomplete(a,x). 

I think we have two possibilities: 
a) return a more general form in terms of meijer-g functions or 
b) check the parameter a to be positive.

I am interested in solution b) and have added some code to sdiffgrad
with the following functionality:

1. Allow a lambda expression to define a derivative.
2. Allow the value NIL, when we do not know the derivative.

I have taken the main idea from the code for the integral-property. 
This is an example for the grad-property of the Incomplete Gamma
function with the extensions:

(putprop '%gamma_incomplete
  `((a z)
    ,(lambda (a unused)
       (declare (ignore unused))
       (cond ((member ($sign a) '($pos $pz))
              ;; The derivative wrt a in terms of                 
              ;; hypergeometric_regularized 2F2
              ;; function and the Generalized Incomplete Gamma function 
              ;; (functions.wolfram.com), only for a>0.
              '((mplus)
                 ((mtimes)
                   ((mexpt) ((%gamma) a) 2)
                   ((mexpt) z a)
                   (($hypergeometric_regularized)
                      ((mlist) a a)
                      ((mlist) ((mplus) 1 a) ((mplus) 1 a))
                      ((mtimes) -1 z)))
                 ((mtimes) -1
                   ((%gamma_incomplete_generalized) a 0 z)
                   ((%log) z))
                 ((mtimes)
                   ((%gamma) a)
                   ((mqapply) (($psi array) 0) a))))
             (t
              ;; No derivative. Maxima generates a noun form.  
              nil)))
    ;; The derivative wrt z
    ((mtimes) -1
      ((mexpt) $%e ((mtimes) -1 z))
      ((mexpt) z ((mplus) -1 a))))
  'grad)

Now we get the following:

The sign of the parameter a is not known. Return a noun form:

(%i3) diff(gamma_incomplete(a,x),a);
(%o3) 'diff(gamma_incomplete(a,x),a,1)

The parameter a has a positive sign:

(%i4) assume(a>0)$

(%i5) diff(gamma_incomplete(a,x),a);
(%o5) -(gamma(a)-gamma_incomplete(a,x))*log(x)
       +gamma(a)^2*hypergeometric_regularized([a,a],[a+1,a+1],-x)*x^a
       +psi[0](a)*gamma(a)

Now the second derivative will work too:

(%i6) diff(%,a);
(%o6) -log(x)*((gamma(a)-gamma_incomplete(a,x))*log(x)
              -gamma(a)^2*hypergeometric_regularized([a,a],[a+1,a
+1],-x)*x^a)
       +gamma(a)^2*hypergeometric_regularized([a,a],[a+1,a
+1],-x)*x^a*log(x)
       +gamma(a)^2*'diff(hypergeometric_regularized([a,a],[a+1,a
+1],-x),a,1)
                  *x^a
       +2*psi[0](a)*gamma(a)^2*hypergeometric_regularized([a,a],[a+1,a
+1],-x)
         *x^a+psi[1](a)*gamma(a)+psi[0](a)^2*gamma(a)

The following gives a noun form too (the sign of x is not known):

(%i7) diff(gamma_incomplete(x,x),x);
(%o7) 'diff(gamma_incomplete(x,x),x,1)

I have tried the extension for other functions too. E. g. for the
jacobi_p function, with the following grad-property:

(putprop '$jacobi_p
         '((n a b x)
           nil
           nil
           nil
           ((mtimes)
            ((mexpt) ((mplus ) a b ((mtimes ) 2 n)) -1)
            ((mplus)
             ((mtimes) 2
              (($unit_step) n)
              ((mplus) a n) ((mplus) b n)
              (($jacobi_p) ((mplus) -1 n) a b x))
             ((mtimes) n (($jacobi_p) n a b x)
              ((mplus) a ((mtimes ) -1 b)
               ((mtimes)
                ((mplus) ((mtimes) -1 a) ((mtimes ) -1 b)
                 ((mtimes) -2 n)) x))))
            ((mexpt) ((mplus) 1 ((mtimes) -1 ((mexpt ) x 2))) -1)))
         'grad)

Now we get a noun form and not a Maxima Error for unknown derivatives:

(%i11) diff(jacobi_p(n,a,b,x),n);
(%o11) 'diff(jacobi_p(n,a,b,x),n,1)
(%i12) diff(jacobi_p(n,a,b,x),a);
(%o12) 'diff(jacobi_p(n,a,b,x),a,1)
(%i13) diff(jacobi_p(n,a,b,x),b);
(%o13) 'diff(jacobi_p(n,a,b,x),b,1)
(%i14) diff(jacobi_p(n,a,b,x),x);
(%o14) (n*jacobi_p(n,a,b,x)*((-2*n-b-a)*x-b+a)
        +2*(n+a)*(n+b)*jacobi_p(n-1,a,b,x)*unit_step(n))
        /((2*n+b+a)*(1-x^2))

Only two small extension to the function sdiffgrad are necessary to get
the above results.

Do you think that this is an interesting extension to Maxima?

Remark: Because we have new code to support the hypergeometric functions
(thanks to Barton), it is possible to improve the derivatives of the
gamma_incomplete and of other special functions further.

Dieter Kaiser