From: lph@paradigm.com
Newsgroups: sci.math.symbolic
Subject: Re: Gaussian moment integral (long Paramax soln)
Date: 23 Jul 92 13:03:29 GMT
Organization: Paradigm Associates Inc, Cambridge MA

Richard Waterman's query regarding the integral of x^a*exp(-b*x^2-c*x)
wrt x from 0 to inf can be handled several different ways. The general
case of the integral is not known to the definite integrator, but one
can use a=2 (or any non-negative integer) to get a closed-form answer
from the integrator. The Laplace transform in the general case is also
not known, but for any non-negative integer a the answer can be
obtained from differentiating the LT of exp(-b*x^2) from x into c wrt
c a total of a times. Also, a hypergeometric integrator can be called
for the case of c positive (the portion of the program for c negative
was not completed, apparently) yielding an answer in terms of the
Whittaker functions, and can be shown to give the same Taylor
series as the Laplace answer. Roger Alexander's results in Maple are
confirmed.

/* PARAMAX version 1 maintained by Paradigm Associates, Inc.
under VAX LISP version V2.2
on DEC VAX VUV2 running VMS
Written by LPH at Thu 7/23/92 12:05:22. */

/* can't get a closed-form but at least it asks the right questions. */
(C2) integrate(x^a*exp(-b*x^2-c*x),x,0,inf);
Is  A + 1  positive, negative, or zero?
P;
Is  B  positive, negative, or zero?
P;
Is  A  an integer?
YES;
			    INF
			   /		  2
			   [	 A   - B X  - C X
(D2) 			   I    X  %E		  dX
			   ]
			   /
			    0

/* remove the linear term in the exponential */
(C3) changevar(%,x=y-c/2/b,y,x);
					      2	 2    2
		      INF		   4 B  Y  - C
		     /			 - ------------
		     [		     A	       4 B
		     I	  (2 B Y - C)  %E	        dY
		     ]
		     /
		       C
		      ---
		      2 B
(D3) 		     -------------------------------------
				      A	 A
				     2  B

/* take the sub-case a=2 to allow closed-form computation */
(C4) %,a=2$

/* get the closed-form for a=2 */
(C5) %,integrate;

Is  B  positive or negative?
P;
Is  C  positive or negative?
N;
				      2
				     C
				     ---
			  2	     4 B	   C
      SQRT(%PI) SQRT(B) (C  + 2 B) %E    ERF(- ---------) - 2 B C
					       2 SQRT(B)
(D5) (-----------------------------------------------------------
				  2 B


								     2
								    C
								    ---
							 2	    4 B
					     SQRT(%PI) (C  + 2 B) %E	     2
					   + --------------------------)/(4 B )
						     2 SQRT(B)

/* confirms Alexander's result from Maple */
(C6) %,b=3/2,c=-1,numer;

(D6) 			       0.656798012459442

/* can't do the general case by Laplace, either */
(C7) laplace(x^a*exp(-b*x^2),x,c);

Is  A + 1  positive, negative, or zero?
P;
Is  B  positive, negative, or zero?
P;
Is  A  an integer?
YES;
					    2
				   A   - B X
(D7) 			  LAPLACE(X  %E	     , X, C)

/* express the definition of the Laplace transform from x to c */
(C8) %='integrate(x^a*exp(-b*x^2)*exp(-c*x),x,0,inf);

					  INF
			      2		 /	        2
		     A	 - B X		 [     A   - B X  - C X
(D8) 	    LAPLACE(X  %E      , X, C) = I    X  %E	        dX
					 ]
					 /
					  0

/* replace a with 0 */
(C9) %,a=0;

					  INF
			      2		 /	     2
			 - B X		 [      - B X  - C X
(D9) 	       LAPLACE(%E      , X, C) = I    %E	     dX
					 ]
					 /
					  0

/* now the Laplace transform can be carried out in closed form */
(C10) %,laplace;

Is  B  positive or negative?
P;
		      2
		     C
		     ---
		     4 B	      C		 INF
	 SQRT(%PI) %E    (1 - ERF(---------))   /	    2
				  2 SQRT(B)     [      - B X  - C X
(D10) 	 ------------------------------------ = I    %E		    dX
		      2 SQRT(B)		        ]
						/
						 0

/* the a=2 case is merely the second derivative wrt c */
(C11) diff(%,c,2);

		      2
		     C
		     ---
		 2   4 B	      C
      SQRT(%PI) C  %E    (1 - ERF(---------))
				  2 SQRT(B)
(D11) ---------------------------------------
			 5/2
		      8 B


		     2
		    C
		    ---
		    4 B		     C		       INF
	SQRT(%PI) %E    (1 - ERF(---------))	      /		     2
				 2 SQRT(B)      C     [	    2   - B X  - C X
      + ------------------------------------ - ---- = I	   X  %E	     dX
			  3/2			  2   ]
		       4 B		       4 B    /
						       0

/* the a=2 case of the Laplace transform */
(C12) d8,a=2;

					  INF
			      2		 /	        2
		     2	 - B X		 [     2   - B X  - C X
(D12) 	    LAPLACE(X  %E      , X, C) = I    X  %E	        dX
					 ]
					 /
					  0

/* two things equal to a third are equal to each other */
(C13) subst(reverse(d11),%);

						   2
						  C
						  ---
					      2	  4 B		   C
			2	   SQRT(%PI) C  %E    (1 - ERF(---------))
	       2   - B X				       2 SQRT(B)
(D13) LAPLACE(X  %E	 , X, C) = ---------------------------------------
						      5/2
						   8 B


						 2
						C
						---
						4 B		 C
				    SQRT(%PI) %E    (1 - ERF(---------))
							     2 SQRT(B)	    C
				  + ------------------------------------ - ----
						      3/2		      2
						   4 B			   4 B

/* get the right hand side */
(C14) rhs(%)$

/* factor the first two terms together */
(C15) factor(part(%,[1,2]))+rest(%,2);

Is  B  positive, negative, or zero?
P;

				     2
				    C
				    ---
			 2	    4 B		 C
	     SQRT(%PI) (C  + 2 B) %E    (ERF(---------) - 1)
					     2 SQRT(B)	        C
(D15) 	   - ----------------------------------------------- - ----
				    5/2				  2
				 8 B			       4 B

/* the numerical values, also confirms Alexander's Maple result */
(C16) %,b=3/2,c=-1;

			     1/6		   SQRT(2)
		 2 SQRT(2) %E    SQRT(%PI) (- ERF(---------) - 1)
	     1					  2 SQRT(3)
(D16) 	     - - ------------------------------------------------
	     9			    9 SQRT(3)

/* same as before */
(C17) %,numer;

(D17) 			       0.656798012459442

/* now use the hypergeometric integrator */
(C18) specint(x^a*exp(-b*x^2-c*x),x);

Is  C  positive, negative, or zero?
N;

/* doesn't handle this case as of yet */
(D18) 		        OTHER-DEFINT-TO-FOLLOW-NEGTEST

/* ok, get an answer for positive c */
(C19) ''c18;

Is  C  positive, negative, or zero?
P;
			     2
			    C
			    ---
       A + 1		    8 B
(D19) 2	     GAMMA(A + 1) %E


	     - A - 1
	     ------- + 1/2			        2
		2	    1/4			       C
  SQRT(%PI) 2		   B    %M		      (---)
				  - A - 1	       4 B
				  ------- + 1/4, - 1/4
				     2
 (---------------------------------------------------------
		       1   - A - 1
		 GAMMA(- - -------) SQRT(C)
		       2      2


		- A - 2
		------- + 1			       2
		   2	     1/4		      C
   2 SQRT(%PI) 2	    B    %M		     (---)
				   - A - 1	      4 B     A + 1  A + 1
				   ------- + 1/4, 1/4	      -----  -----
				      2			        2      2
 - -------------------------------------------------------)/(8	    B	  )
			  - A - 1
		  GAMMA(- -------) SQRT(C)
			     2

/* use values of a and b */
(C20) %,a=2,b=3/2$

/* the answer in c in terms of Whittaker functions */
(C21) factor(%);

	       2
	      C
	      --		  2				 2
	      12		 C			        C
	    %E   (4 %M	        (--) - SQRT(%PI) %M	       (--))
		      - 5/4, 1/4 6		   - 5/4, - 1/4 6
(D21) 	  - --------------------------------------------------------
				 1/4  1/4
			      3 2    3	  SQRT(C)

/* use the Abramowitz & Stegun formula 13.1.32 */
(C22) %,%m[a,b](z):=exp(-z/2)*z^(1/2+b)*m(1/2+b-a,1+2*b,z)$

/* c is positive now */
(C23) %,assume_pos:true;

				 2				        2
				C				       C
		    2	      - --		        2	     - --     2
		3  C    3/2     12		 3  1  C	       12    C
	 4 M(2, -, --) C    %E	     SQRT(%PI) M(-, -, --) SQRT(C) %E	     --
		2  6				 2  2  6		     12
	(------------------------- - ------------------------------------) %E
		    3/4				      1/4
		   6				     6
(D23) - -----------------------------------------------------------------------
				     1/4  1/4
				  3 2    3    SQRT(C)

/* use A&S formula 13.1.2 for the Kummer function */
(C24) %,m(a,b,z):=sum(gamma(a+n)*gamma(b)*z^n/(gamma(a)*gamma(b+n)*n!),n,0,inf)$

/* use c positive */
(C25) %,assume_pos:true;

				    2
				   C  INF
				 - -- ====   2 N
			   3/2	   12 \	    C    GAMMA(N + 2)
	   2  2 SQRT(%PI) C    %E      >    ------------------
	  C			      /	     N		    3
	  --			      ====  6  N! GAMMA(N + -)
	  12			      N = 0		    2
(D25) - %E   (------------------------------------------------
				     3/4
				    6


			      2
			     C  INF    2 N	     3
			   - -- ====  C	   GAMMA(N + -)
			     12 \		     2
     2 SQRT(%PI) SQRT(C) %E      >    ------------------
				/      N	      1
				====  6  N! GAMMA(N + -)
				N = 0		      2	      1/4  1/4
   - ---------------------------------------------------)/(3 2	  3    SQRT(C))
			     1/4
			    6

/* factor the expression */
(C26) factor(%)$

/* put terms inside of sums */
(C27) intosum(multthru(%))$

/* combine into one sum, note there are odd and even terms in c */
(C28) sumcontract(%);

      INF		  2 N  - N - 1/4       2 N + 3
      ====   2 SQRT(%PI) C    6	         GAMMA(-------)
      \						  2
(D28)  >    (------------------------------------------
      /		      1/4  1/4	        2 N + 1
      ====	   3 2	  3    N! GAMMA(-------)
      N = 0				   2


					      2 N + 1  - N - 3/4
				 2 SQRT(%PI) C	      6	         GAMMA(N + 2)
			       - --------------------------------------------)
					   1/4  1/4	     2 N + 3
					3 2    3    N! GAMMA(-------)
								2

/* factor the whole thing */
(C29) factor(%);

		   INF
		   ====
		   \	  2 N  - N - 3/4
(D29) 2 SQRT(%PI) ( >	 C    6
		   /
		   ====
		   N = 0


	       2 2 N + 3			 2 N + 1
 (SQRT(6) GAMMA (-------) - C GAMMA(N + 2) GAMMA(-------))
		    2				    2


	   2 N + 1	  2 N + 3	 1/4  1/4
/(N! GAMMA(-------) GAMMA(-------)))/(3 2    3	 )
	      2		     2

/* get the first several terms of the taylor series around c=0 */
(C30) taylor(%,c,0,4);

		  1/4   1/4 3   1/4 3		    1/4	  1/4 3	  1/4 3
	 SQRT(6) 6    (2   )  (3   )  SQRT(%PI)	  (6    (2   )  (3   ) ) C
(D30)/T/ -------------------------------------- - ------------------------
			  108				     27


	     1/4   1/4 3   1/4 3	     2	     1/4   1/4 3   1/4 3   3
   (SQRT(6) 6    (2   )  (3   )  SQRT(%PI)) C	 (2 6    (2   )  (3   ) ) C
 + ------------------------------------------- - ---------------------------
		       216				     243


	       1/4   1/4 3   1/4 3	       4
   (5 SQRT(6) 6	   (2   )  (3   )  SQRT(%PI)) C
 + --------------------------------------------- + . . .
		       7776


/* simplify each term with radcan because of the radicals */
(C31) map('radcan,%);

				   4	  3				 2
      5 SQRT(2) SQRT(3) SQRT(%PI) C    4 C    SQRT(2) SQRT(3) SQRT(%PI) C
(D31) ------------------------------ - ---- + ----------------------------
		   1296		        81		   36


						2 C   SQRT(2) SQRT(3) SQRT(%PI)
					      - --- + -------------------------
						 9		 18

/* the direct laplace transform in question */
(C32) laplace(x^2*exp(-3/2*x^2),x,c)$

/* taylor around c=0 to 4 terms */
(C33) taylor(%,c,0,4);

									2
	 SQRT(2) SQRT(3) SQRT(%PI)   2 C   (SQRT(2) SQRT(3) SQRT(%PI)) C
(D33)/T/ ------------------------- - --- + ------------------------------
		    18		      9			 36


				   3				      4
				4 C    (5 SQRT(2) SQRT(3) SQRT(%PI)) C
			      - ---- + -------------------------------- + . . .
				 81		     1296

/* laplace transform answer is identical to the hypergeometric integrator answer
through c^4, at least */
(C34) %-d31;

(D34)/T/ 			   0 + . . .

Leo Harten
Paradigm Associates, Inc.
29 Putnam Avenue, Suite 6
Cambridge, MA 02139/USA
LPH@PARADIGM.COM
+1(617)492-6079