#833 specint(gammaincomplete(v,a*t)*exp(-p*t),t) seems wrong

[Raymond Toy]

(%i12) assume(v>0,p>0);

(%o12) [v>0,p>0]
(%i13) specint(gammaincomplete(v,a*t)*exp(-p*t),t);

SIMP2F1-WILL-CONTINUE-IN
(%o13)
a^v*(p+a)^(-v-1)*gamma(v+1)*%f[2,1]([v+1,3/2],[2],p/(p+a))

Compare this with

(%i14) specint(gammagreek(v,a*t)*exp(-p*t),t);
(%o14) a^v*p^(-v-1)*gamma(v+1)/((a/p+1)^v*v)

This matches formula 34, p 179 in Tables of Transforms.
Considering gammaincomplete = gamma(n)-gammagreek, the
expression for gammaincomplete seems wrong.

It might still be right if the hypergeometric function
simplifies, but maxima can't, and I can't think of any
way to simplify it either.

[Raymond Toy]

Logged In: YES
user_id=28849

The fix for transforming %w causes this return the result in
terms of the associated Legendre function Q. Somewhat
better, but I do not know if this is equivalent.

[Dieter Kaiser]

Logged In: YES
user_id=2039760
Originator: NO

I would like to suggest the following algorithm to get an evaluated integral for the Gammaincomplete function:

(defun gammaincomplete-to-hypergeometric (a x)
(if (eq (asksign a) '$zero)

;; The representation as Whittaker W function for a=0

(mul*
(power x (div (sub a 1) 2))
(power '$%e (div x -2))
($whittaker_w (div (sub a 1) 2) (div a 2) x)) ; sign of (div a 2) correct?

;; In all other cases the representation as Hypergeometric 1F1 function

(sub
(list '(%gamma) a)
(mul
(power x a)
(inv a)
($hypergeometric
(list '(mlist) a)
(list '(mlist) (add a 1))
(mul -1 x))))))

This routine has to be replaced for GAMMAINCOMPLETETW and called in the routine FRACTEST2. The idea is to use the Hypergeometric Function 1F1 when the parameter A is not equal to zero. That is equivalent to use gamma(a) - gammagreek(a,x) for the evaluation. Now we use the transformation to Whittaker W only for a=0, but for this case we get a nice result too.

Here some examples with the above routine:

(%i1) assume(s>0,a>0);
(%o1) [s > 0,a > 0]

We get the expected simple answer for gammaincomplete:

(%i2) radcan(specint(%e^(-s*t)*gammaincomplete(a,t),t));
(%o2) (a*gamma(a)*(s+1)^a-gamma(a+1))/(a*s*(s+1)^a)

If we take a=0 we get a nice and correct result too:

(%i3) radcan(specint(%e^(-s*t)*gammaincomplete(0,t),t));
(%o3) log(s+1)/s

The function %ei uses the code of gammaincomplete with a=0:

(%i4) radcan(specint(%e^(-s*t)*%ei(t),t));
(%o4) -log(1-s)/s

For special values like a=1/2 or a=1 both algorithm get simple expressions. The results calulated with 1F1 are identical to the results of the orginal code which use the transformation to Whittaker W:

(%i5) radcan(specint(%e^(-s*t)*gammaincomplete(1/2,t),t));
(%o5) (sqrt(%pi)*sqrt(s+1)-sqrt(%pi))/(s*sqrt(s+1))

(%i6) radcan(specint(%e^(-s*t)*gammaincomplete(1,t),t));
(%o6) 1/(s+1)

In the above code I have used the routines $hypergeometric and $whittaker_w. I have introduced this functions to make the code more expressive. It is easy to use the alternavtives %f and %w. Furthermore $specint has to be extended to integrate %f.

The testsuite has no problems with this alogorithm.

Dieter Kaiser

[Dieter Kaiser]

Logged In: YES
user_id=2039760
Originator: NO

I have suggested to use a representation as a Hypergeometric 1F1 for a<>0 for the Gammaincomplete function to get more simple results for the Laplace transform of the Gammaincomplete function.

Because $specint has no code to do the Laplace transform for the Hypergeometric function itself ($specint knows only an intern representation for the Hypergeometric function) it is more easy to transform to the Gammagreek function. The Gammagreek function will be transformed to the Hypergeometric 1F1 function.

The following changes produce the correct results shown in this thread:

===================================================================
RCS file: /cvsroot/maxima/maxima/src/hypgeo.lisp,v
retrieving revision 1.40
diff -u -r1.40 hypgeo.lisp
--- hypgeo.lisp 4 Jul 2008 14:12:52 -0000 1.40
+++ hypgeo.lisp 18 Jul 2008 22:03:29 -0000
@@ -2301,7 +2302,7 @@
((eq flg 'kbateman)
(kbatemantw i1 a1))
((eq flg 'gammaincomplete)
- (gammaincompletetw a1 i1))
+ (gammaincomplete-to-gammagreek a1 i1))
((eq flg 'kti)
(kti i1 a1))
((eq flg 'erfc)
@@ -2582,6 +2583,28 @@
(power '$%e (div x -2))
(wwhit x (div (sub a 1) 2)(div a 2))))

+;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
+;;; Only for a=0 we use the general representation as a Whittaker W function.
+;;; In all other cases we transform to a Gammagreek function. The Gammagreek
+;;; function will be further transformed to a Hypergeometric 1F1 representation.
+;;; With this change we get more simple and correct results for the Laplace
+;;; transform of the Gammaincomplete function.
+;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
+
+(defun gammaincomplete-to-gammagreek (a x)
+ (if (eq (asksign a) '$zero)
+ ;; The representation as a Whittaker W function for a=0
+ (mul
+ (power x (div (sub a 1) 2))
+ (power '$%e (div x -2))
+ (wwhit x (div (sub a 1) 2) (div a 2)))
+ ;; In all other cases the representation as a Gammagreek function
+ (sub
+ (list '(%gamma) a)
+ (list '($gammagreek) a x))))
+
+;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
+
(defun distrexecinit (fun)
(cond ((and (consp fun)
(consp (car fun))

===================================================================

The testsuite has no problems with this algorithm. With this changes this bug report could be closed. Should I commit this change of the code?

Dieter Kaiser

[Dieter Kaiser]

Logged In: YES
user_id=2039760
Originator: NO

The code is changed as suggested.
Closing the bug report.

Dieter Kaiser