maxima-bugs

[Maxima-bugs] [ maxima-Bugs-2867727 ] specint: wrong result for Parabolic Cylinder D function

[SourceForge.net <no...@so...>]

Bugs item #2867727, was opened at 2009-09-27 01:30
Message generated for change (Tracker Item Submitted) made by crategus
You can respond by visiting: 
https://sourceforge.net/tracker/?func=detail&atid=104933&aid=2867727&group_id=4933

Please note that this message will contain a full copy of the comment thread,
including the initial issue submission, for this request,
not just the latest update.
Category: Lisp Core - Integration
Group: None
Status: Open
Resolution: None
Priority: 5
Private: No
Submitted By: Dieter Kaiser (crategus)
Assigned to: Nobody/Anonymous (nobody)
Summary: specint: wrong result for Parabolic Cylinder D function

Initial Comment:
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

----------------------------------------------------------------------

You can respond by visiting: 
https://sourceforge.net/tracker/?func=detail&atid=104933&aid=2867727&group_id=4933

[Maxima-bugs] [ maxima-Bugs-2867727 ] specint: wrong result for Parabolic Cylinder D function

[SourceForge.net <no...@so...>]

Bugs item #2867727, was opened at 2009-09-27 01:30
Message generated for change (Comment added) made by crategus
You can respond by visiting: 
https://sourceforge.net/tracker/?func=detail&atid=104933&aid=2867727&group_id=4933

Please note that this message will contain a full copy of the comment thread,
including the initial issue submission, for this request,
not just the latest update.
Category: Lisp Core - Integration
Group: None
Status: Open
Resolution: None
Priority: 5
Private: No
Submitted By: Dieter Kaiser (crategus)
Assigned to: Nobody/Anonymous (nobody)
Summary: specint: wrong result for Parabolic Cylinder D function

Initial Comment:
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

----------------------------------------------------------------------

>Comment By: Dieter Kaiser (crategus)
Date: 2009-09-27 22:37

Message:
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

----------------------------------------------------------------------

You can respond by visiting: 
https://sourceforge.net/tracker/?func=detail&atid=104933&aid=2867727&group_id=4933

[Maxima-bugs] [ maxima-Bugs-2867727 ] specint: wrong result for Parabolic Cylinder D function

[SourceForge.net <no...@so...>]

Bugs item #2867727, was opened at 2009-09-27 01:30
Message generated for change (Comment added) made by crategus
You can respond by visiting: 
https://sourceforge.net/tracker/?func=detail&atid=104933&aid=2867727&group_id=4933

Please note that this message will contain a full copy of the comment thread,
including the initial issue submission, for this request,
not just the latest update.
Category: Lisp Core - Integration
Group: None
Status: Open
Resolution: None
Priority: 5
Private: No
Submitted By: Dieter Kaiser (crategus)
Assigned to: Nobody/Anonymous (nobody)
Summary: specint: wrong result for Parabolic Cylinder D function

Initial Comment:
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

----------------------------------------------------------------------

>Comment By: Dieter Kaiser (crategus)
Date: 2009-09-28 00:14

Message:
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

----------------------------------------------------------------------

Comment By: Dieter Kaiser (crategus)
Date: 2009-09-27 22:37

Message:
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

----------------------------------------------------------------------

You can respond by visiting: 
https://sourceforge.net/tracker/?func=detail&atid=104933&aid=2867727&group_id=4933

[Maxima-bugs] [ maxima-Bugs-2867727 ] specint: wrong result for Parabolic Cylinder D function

[SourceForge.net <no...@so...>]

Bugs item #2867727, was opened at 2009-09-27 01:30
Message generated for change (Settings changed) made by crategus
You can respond by visiting: 
https://sourceforge.net/tracker/?func=detail&atid=104933&aid=2867727&group_id=4933

Please note that this message will contain a full copy of the comment thread,
including the initial issue submission, for this request,
not just the latest update.
Category: Lisp Core - Integration
Group: None
>Status: Closed
>Resolution: Fixed
Priority: 5
Private: No
Submitted By: Dieter Kaiser (crategus)
Assigned to: Nobody/Anonymous (nobody)
Summary: specint: wrong result for Parabolic Cylinder D function

Initial Comment:
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

----------------------------------------------------------------------

>Comment By: Dieter Kaiser (crategus)
Date: 2009-09-28 23:41

Message:
Fixed in hypgeo.lisp revision 1.67,
Closing this bug report.

Dieter Kaiser

----------------------------------------------------------------------

Comment By: Dieter Kaiser (crategus)
Date: 2009-09-28 00:14

Message:
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

----------------------------------------------------------------------

Comment By: Dieter Kaiser (crategus)
Date: 2009-09-27 22:37

Message:
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

----------------------------------------------------------------------

You can respond by visiting: 
https://sourceforge.net/tracker/?func=detail&atid=104933&aid=2867727&group_id=4933