Extension of the integrator - special functions

To support integrals for special functions I have done four modifications of
sinint. The testsuite has no problems with these modifications. To get more
results it is also necessary to have more integrals with power functions.

At first as an example some integrals for the function expintegral_ei:

(%i17) integrate(expintegral_ei(z),z);
(%o17) z*expintegral_ei(z)-%e^z
(%i18) integrate(z*expintegral_ei(z),z);
(%o18) z^2*expintegral_ei(z)/2-(z-1)*%e^z/2
(%i19) integrate(z^(-2)*expintegral_ei(z),z);
(%o19) gamma_incomplete(-1,-z)-expintegral_ei(z)/z

Short description of the necessary modifications:

1. Put the integral of the elementary function on the property list.
2. Extend the routine integrallookups with an algorithm to look up 
   the integral from the property list.
3. Write a more general routine to do a partial integration and call this
   routine in the routine intform.
4. Generalize diffdiv to support the integration of functions with more
   than one argument (to support gamma_incomplete and expintegral_e).

The code below is a first step. Maxima can solve only simple integrals with
special functions. I have tested it for the Exponential Integrals and the
Incomplete Gamma function. It is also possible to cut off the special code for
the logarithm function and to use the more general approach.

Further generalizizations and extensions are possible.

Remark:

The code of sinint is not easy to read and it is not always easy to see what is
going on. One problem is the use of a lot of global variables in combination
with a highly recursive algorithm.

The new routine partial-integration and the routine ratlog do almost the same.
If you compare the routines you can see that no global variables are necessary.

The only thing to do to give support for new functions is to put the elmentary
integral on the property list. This is an example for the Exponential Integral
Ei:

;;; Integral of Exponential Ei

(defprop %expintegral_ei
  ((x)
   ((mplus) 
      ((mtimes) -1 ((mexpt) $%e x))
      ((mtimes) x ((%expintegral_ei) x))))
  integral)

Do you think we should start with this extension to support more integrals with
special functions?

Dieter Kaiser


This is a diff with the changes to support the special functions:

Index: sin.lisp
===================================================================
RCS file: /cvsroot/maxima/maxima/src/sin.lisp,v
retrieving revision 1.33
diff -u -r1.33 sin.lisp
--- sin.lisp   18 Nov 2008 13:02:32 -0000   1.33
+++ sin.lisp   22 Nov 2008 15:42:42 -0000
(defmacro op (frob)
   `(get ,frob 'operators))
 
@@ -135,6 +138,22 @@
            (cdr (sassq 'c y 'nill))
            (cdr (sassq 'd y 'nill)))))
        (t (return nil)))))))
+
+      ;; We have a special function with an integral on the property list.
+      ;; We try to do a partial integration.
+
+      ((and (not (atom (car expres)))
+            (get (caar expres) 'integral))
+       (when *debug-integrate* 
+         (format t "~&INTFORM with Integral on property list~%"))
+       (cond
+         ((setq arg
+            (m2 exp 
+             `((mtimes) ((,(caar expres)) (b rat8)) ((coefftt) (c rat8prime)))
+              nil))
+          (partial-integration (cons (cons 'a expres) arg) var))
+    (t nil)))
+
    ((optrig (caar expres))
     (cond ((not (setq w (m2 (cadr expres) *c* nil)))
       (intform (cadr expres)))

@@ -498,6 +529,38 @@
              '((mplus)
                ((mtimes) x ((%log) x))
                ((mtimes) -1 x))))
+
+        ;; The integral of the Log function is directly implemented in the
+        ;; algorithm. This can be generalized to a lookup algorithm for any
+        ;; special function. The integral is put on the property list.
+        ;; In a first step we support functions with one argument and add a
+        ;; lookup algorithm for the functions gamma_incomplete and 
+        ;; expintegral_e with two arguments.
+
+        ((and (not (atom (car exp)))
+              (setq form (get (caar exp) 'integral))
+              (= (length (car form)) 1))
+         ;; Found an integral on the property list with one argument.
+         ($substitute (cadr exp) (caar form) (cadr form)))
+
+        ((and (not (atom (car exp)))
+              (setq form (get (caar exp) 'integral))
+              (= (length (car form)) 2))
+         ;; Found an integral on the property list with two arguments.
+         (cond
+           ((and (freevar (cadr exp)) ; Integration only wrt the 2. argument
+                 (member (caar exp) '(%gamma_incomplete %expintegral_e)))
+            ;; We have no general algorithm for functions with more than
+            ;; one argument. This is a first workaround for the Incomplete
+            ;; Gamma function and the Exponential Integral E
+            (let ((arg1 (cadr exp))
+                  (arg2 (caddr  exp)))
+              ($substitute
+                arg2
+                (cadar form)
+                ($substitute arg1 (caar form) (caddr form)))))
+           (t nil)))
+
    ((eq (caar exp) 'mplus)
     (muln (list '((rat simp) 1 2) exp exp) nil))
    ((eq (caar exp) 'mexpt)

@@ -1500,6 +1567,26 @@
               (list '(mtimes) y *d*)
               (list '(mtimes) -1 z))))))
 
+;;; partial-integration is an extension of the algorithm of ratlog to support 
+;;; the technique of partial integration for more cases. The integrand
+;;; is like g(x)*f'(x) and the result is g(x)*f(x)-integrate(g'(x)*f(x),x).
+
+(defun partial-integration (form var)
+  (let ((g  (cdr (assoc 'a form)))   ; part g(x)
+        (df (cdr (assoc 'c form)))   ; part f'(x)
+        (f  nil))
+
+    (setq f (integrator df var))     ; integrate f'(x) wrt var
+
+    (cond 
+      ((or (isinop f '%integrate)    ; no result or
+           (isinop f (caar g)))      ; g in result
+       nil)                          ; we return nil
+      (t
+       ;; Build the result: g(x)*f(x)-integrate(g'(x)*f(x))
+       (add (mul f g)
+            (mul -1 (integrator (mul f (sdiff g var)) var)))))))
+
 ;; returns t if argument of every trig operation in y matches arg
 (defun every-trigarg-alike (y arg)
   (cond ((atom y) t)

@@ -1595,7 +1682,19 @@
           (cond ((freevar (cadr y)) (caddr y))
            ((freevar (caddr y)) (cadr y))
            (t 0)))
-         (t (cadr y))))
+                          ;; At this point the algorithm assume that we have a
+                          ;; function with one argument. This argument is taken
+                          ;; with (cadr y). The algorithm has to be extended
+                          ;; to support functions with more than one argument.
+                          ;; In a first step we only support the functions
+                          ;; gamma_incomplete and expintegral_e.
+                          ((member (caar y) '(%gamma_incomplete
%expintegral_e))
+                           ;; We assume the integration wrt the second argument
+                           ;; of the function.
+                           (caddr y))
+                          (t
+                           ;; Take the argument of a function with one
variable.
+                           (cadr y))))