Menu
▾
▴
#1777 specint: wrong result for Parabolic Cylinder D function
A look at the implemented algorithm for the Laplace transform of a hypergeometric function shows that the integral for the Parabolic Cylinder D function should work only for an argument t^k with an exponent k < 1, e.g with an argument sqrt(t) (I have already introduced a new symbol in my sandbox for the Parabolic Cylinder D function):
(%i2) factor(ratsimp(specint(exp(-s*t)*parabolic_cylinder_d(0,sqrt(t)),t)));
(%o2) 4/(4*s+1)
This is a correct result.
We can verify the result using a specific expansion:
parabolic_cylinder_d(0,z) = exp(-z^2/4)
parabolic_cylinder_d(1,z) = z*exp(-z^2/4)
parabolic_cylinder_d(2,z) = (z^2-1)*exp(-z^2/4)
Now some more results with the first expansion for v=0:
(%i3) d0(z):=exp(-z^2/4)$
We can verify the result for the Parabolic Cylinder D function with an argument sqrt(t):
(%i6) factor(ratsimp(specint(exp(-s*t)*d0(sqrt(t)),t)));
(%o6) 4/(4*s+1)
Now we do the integration with an argument t using the expansion. We get a result which contains the Error function:
(%i7) factor(ratsimp(specint(exp(-s*t)*d0(t),t)));
(%o7) -sqrt(%pi)*%e^s^2*(erf(s)-1)
The following integral should not work, but it does and the answer is wrong:
(%i8) factor(ratsimp(specint(exp(-s*t)*parabolic_cylinder_d(0,t),t)));
(%o8) 2^(3/2)*gamma(3/4)/(4*s+1)^(3/4)
The routine f16p217test does not fail as expected and produces the wrong result.
Remarks:
1. $specint uses the Parabolic Cylinder D function and its Laplace transform for the Hermite and the Erfc function. For the Hermite function the algorithm does not work as expected. For the erfc function it is more simple to integrate (1-erf(z)).
2. The Parabolic Cylinder D function is transformed in terms of the Whittaker W function. In a second step the Whittaker W function is transformed to a hypergeometric representation. This representation is integrated.
But we already have a routine simpdtf which does the transformation to a hypergeometric representation in one step.
Dieter Kaiser
The error is in the routine f16p217test. The Laplace transform of the Whittaker W function is only valid for an argument a*t and not a*t^v. There is no check for v=1.
This is a possible correction of the involved part of the code in the routine f16p217test:
;; Ok, we satisfy the conditions. Now extract the arg.
;; The transformation is only valid for an argument a*t. We have
;; to specialize the pattern to make sure that we satisfy the condition.
(let ((l (m2-a*t a) ; more special pattern for the argument
; (m2 a
; '((mplus)
; ((coeffpt) (f hasvar) (a freevar))
; ((coeffpp) (c zerp)))
; nil)
))
If we correct this part of code the Laplace transform for the Parabolic Cylinder D function for an argument sqrt(t) will work, but will fail for an argument t as expected.
Furthermore, we have to change the algorithm of whittest. When f16p217test fails a transformation to a Whittaker M function is done. But the algorithm for the Laplace transform of the Whittaker W function is already the general Laplace transform of the corresponding representation in terms of Whittaker M. The problem is, that we get errors when trying this transformation. So it is the best to return a noun form at this point.
This could be the change to the routine whittest:
(defun whittest (r a i1 i2)
(cond ((f16p217test r a i1 i2))
(t
; The formula used in f16p217test already is already the general Laplace
; transformation. The transformation to a Whittaker M representation
; does not give any new information, but will fail with errors like undefined gamma
; functions or divisions by zero. Therefore, we return a noun form at this point.
; ;; Convert to M function and try again.
; (distrexecinit ($expand (mul (init r)
; (wtm a i1 i2))))
(setq *hyp-return-noun-flag* 'whittest-failed)
)))
With these changes the Laplace transform of the Parabolic Cylinder D function will work as expected. There are no problems with the testsuite.
Dieter Kaiser
I had again a look at the algorithm of whittest. It is not necessary to return in all cases a noun form. We only have to make sure that 1/2-u-k and 1/2+u-k are not zero or a negative integer. Only for these cases the transformation will fail.
This is the extension to check the indices and to call the transformation to the Whittaker M function:
(defun whittest (r a i1 i2)
(let (n m)
(cond ((f16p217test r a i1 i2))
((and
(not
(and (maxima-integerp (setq n (sub (sub '((rat simp) 1 2) i2) i1)))
(member ($sign n) '($zero $neg $nz))))
(not
(and (maxima-integerp (setq m (sub (add '((rat simp) 1 2) i2) i1)))
(member ($sign m) '($zero $neg $nz)))))
;; 1/2-u-k and 1/2+u-k are not zero or a negative integer
;; Transform to Whittaker M and try again.
(distrexecinit ($expand (mul (init r) (wtm a i1 i2)))))
(t
;; Both conditions fails, return a noun form.
(setq *hyp-return-noun-flag* 'whittest-failed)))))
Dieter Kaiser
Fixed in hypgeo.lisp revision 1.67,
Closing this bug report.
Dieter Kaiser