Enhanced Laplace transforms and desolve for Maxima



On 25/01/2011 7:49 AM, Barton Willis wrote:
> There is a symbol property mechanism (I think) for telling Maxima the 
> antiderivative of a function.

Yes.   A simple example - in gamma.lisp.  nil is returned as the 
integral wrt the first argument is unknown.

    ;;; Integral of the Incomplete Gamma function

    (defprop %gamma_incomplete
       ((a z)
        nil
        ((mplus)
           ((mtimes) -1 ((%gamma_incomplete) ((mplus) 1 a) z))
           ((mtimes) ((%gamma_incomplete) a z) z)))
       integral)

You can also use a function if you use putprop to avoid quoting issues - 
as in simp.lisp -

    (defun abs-integral (x)
       (mul (div 1 2) x (take '(mabs) x)))

    (putprop 'mabs `((x) ,#'abs-integral) 'integral)

or a lambda expression  - as in bessel.lisp.  I should rewrite this.  I 
couldn't get the quoting right when I tried to use a function.  Obvious 
now :-)

    ;; Integral of the Bessel function wrt z
    (putprop '%bessel_j
       `((v z)
        nil
       ,(lambda (v z)
         (case v
               (0
                ;; integrate(bessel_j(0,z)
                ;; = (1/2)*z*(%pi*bessel_j(1,z)*struve_h(0,z)
                ;;            +bessel_j(0,z)*(2-%pi*struve_h(1,z)))
                '((mtimes) ((rat) 1 2) z
                  ((mplus)
                   ((mtimes) $%pi
                    ((%bessel_j) 1 z)
                    ((%struve_h) 0 z))
                   ((mtimes)
                    ((%bessel_j) 0 z)
                    ((mplus) 2 ((mtimes) -1 $%pi
                                ((%struve_h) 1 z)))))))
               (1
                ;; integrate(bessel_j(1,z) = -bessel_j(0,z)
                '((mtimes) -1 ((%bessel_j) 0 z)))
               (otherwise
                ;; http://functions.wolfram.com/03.01.21.0002.01
                ;; integrate(bessel_j(v,z)
                ;;  = 2^(-v-1)*z^(v+1)*gamma(v/2+1/2)
                ;;   *
    hypergeometric_regularized([v/2+1/2],[v+1,v/2+3/2],-z^2/4)
                ;;  =
    2^(-v-1)*z^(v+1)*hypergeometric([v/2+1/2],[v+1,v/2+3/2],-z^2/4)
                ;;   / ((v/2+1/2)*gamma(v+1))
                '((mtimes)
                   (($hypergeometric)
                     ((mlist)
                       ((mplus) ((rat) 1 2) ((mtimes) ((rat) 1 2) v)))
                     ((mlist)
                       ((mplus) ((rat) 3 2) ((mtimes) ((rat) 1 2) v))
                       ((mplus) 1 v))
                     ((mtimes) ((rat) -1 4) ((mexpt) z 2)))
                   ((mexpt) ((mplus) ((rat) 1 2) ((mtimes) ((rat) 1 2)
    v)) -1)
                   ((mexpt) 2 ((mplus) -1 ((mtimes) -1 v)))
                   ((mexpt) ((%gamma) ((mplus) 1 v)) -1)
                   ((mexpt) z ((mplus ) 1 v)))))))
       'integral)