home *** CD-ROM | disk | FTP | other *** search
/ NetNews Usenet Archive 1992 #16 / NN_1992_16.iso / spool / sci / math / symbolic / 2061 < prev    next >
Encoding:
Internet Message Format  |  1992-07-23  |  11.4 KB
  1. Path: sparky!uunet!sun-barr!cs.utexas.edu!swrinde!mips!pacbell.com!decwrl!world!paradigm!lph
  2. From: lph@paradigm.com
  3. Newsgroups: sci.math.symbolic
  4. Subject: Re: Gaussian moment integral (long Paramax soln)
  5. Message-ID: <11001@paradigm.com>
  6. Date: 23 Jul 92 13:03:29 GMT
  7. Organization: Paradigm Associates Inc, Cambridge MA
  8. Lines: 449
  9.  
  10. Richard Waterman's query regarding the integral of x^a*exp(-b*x^2-c*x)
  11. wrt x from 0 to inf can be handled several different ways. The general
  12. case of the integral is not known to the definite integrator, but one
  13. can use a=2 (or any non-negative integer) to get a closed-form answer
  14. from the integrator. The Laplace transform in the general case is also
  15. not known, but for any non-negative integer a the answer can be
  16. obtained from differentiating the LT of exp(-b*x^2) from x into c wrt
  17. c a total of a times. Also, a hypergeometric integrator can be called
  18. for the case of c positive (the portion of the program for c negative
  19. was not completed, apparently) yielding an answer in terms of the
  20. Whittaker functions, and can be shown to give the same Taylor
  21. series as the Laplace answer. Roger Alexander's results in Maple are
  22. confirmed.
  23.  
  24. /* PARAMAX version 1 maintained by Paradigm Associates, Inc.
  25. under VAX LISP version V2.2
  26. on DEC VAX VUV2 running VMS
  27. Written by LPH at Thu 7/23/92 12:05:22. */
  28.  
  29. /* can't get a closed-form but at least it asks the right questions. */
  30. (C2) integrate(x^a*exp(-b*x^2-c*x),x,0,inf);
  31. Is  A + 1  positive, negative, or zero?
  32. P;
  33. Is  B  positive, negative, or zero?
  34. P;
  35. Is  A  an integer?
  36. YES;
  37.                 INF
  38.                /          2
  39.                [     A   - B X  - C X
  40. (D2)                I    X  %E          dX
  41.                ]
  42.                /
  43.                 0
  44.  
  45. /* remove the linear term in the exponential */
  46. (C3) changevar(%,x=y-c/2/b,y,x);
  47.                           2     2    2
  48.               INF           4 B  Y  - C
  49.              /             - ------------
  50.              [             A           4 B
  51.              I      (2 B Y - C)  %E            dY
  52.              ]
  53.              /
  54.                C
  55.               ---
  56.               2 B
  57. (D3)              -------------------------------------
  58.                       A     A
  59.                      2  B
  60.  
  61. /* take the sub-case a=2 to allow closed-form computation */
  62. (C4) %,a=2$
  63.  
  64. /* get the closed-form for a=2 */
  65. (C5) %,integrate;
  66.  
  67. Is  B  positive or negative?
  68. P;
  69. Is  C  positive or negative?
  70. N;
  71.                       2
  72.                      C
  73.                      ---
  74.               2         4 B       C
  75.       SQRT(%PI) SQRT(B) (C  + 2 B) %E    ERF(- ---------) - 2 B C
  76.                            2 SQRT(B)
  77. (D5) (-----------------------------------------------------------
  78.                   2 B
  79.  
  80.  
  81.                                      2
  82.                                     C
  83.                                     ---
  84.                              2        4 B
  85.                          SQRT(%PI) (C  + 2 B) %E         2
  86.                        + --------------------------)/(4 B )
  87.                              2 SQRT(B)
  88.  
  89. /* confirms Alexander's result from Maple */
  90. (C6) %,b=3/2,c=-1,numer;
  91.  
  92. (D6)                    0.656798012459442
  93.  
  94. /* can't do the general case by Laplace, either */
  95. (C7) laplace(x^a*exp(-b*x^2),x,c);
  96.  
  97. Is  A + 1  positive, negative, or zero?
  98. P;
  99. Is  B  positive, negative, or zero?
  100. P;
  101. Is  A  an integer?
  102. YES;
  103.                         2
  104.                    A   - B X
  105. (D7)               LAPLACE(X  %E         , X, C)
  106.  
  107. /* express the definition of the Laplace transform from x to c */
  108. (C8) %='integrate(x^a*exp(-b*x^2)*exp(-c*x),x,0,inf);
  109.  
  110.                       INF
  111.                   2         /            2
  112.              A     - B X         [     A   - B X  - C X
  113. (D8)         LAPLACE(X  %E      , X, C) = I    X  %E            dX
  114.                      ]
  115.                      /
  116.                       0
  117.  
  118. /* replace a with 0 */
  119. (C9) %,a=0;
  120.  
  121.                       INF
  122.                   2         /         2
  123.              - B X         [      - B X  - C X
  124. (D9)            LAPLACE(%E      , X, C) = I    %E         dX
  125.                      ]
  126.                      /
  127.                       0
  128.  
  129. /* now the Laplace transform can be carried out in closed form */
  130. (C10) %,laplace;
  131.  
  132. Is  B  positive or negative?
  133. P;
  134.               2
  135.              C
  136.              ---
  137.              4 B          C         INF
  138.      SQRT(%PI) %E    (1 - ERF(---------))   /        2
  139.                   2 SQRT(B)     [      - B X  - C X
  140. (D10)      ------------------------------------ = I    %E            dX
  141.               2 SQRT(B)                ]
  142.                         /
  143.                          0
  144.  
  145. /* the a=2 case is merely the second derivative wrt c */
  146. (C11) diff(%,c,2);
  147.  
  148.               2
  149.              C
  150.              ---
  151.          2   4 B          C
  152.       SQRT(%PI) C  %E    (1 - ERF(---------))
  153.                   2 SQRT(B)
  154. (D11) ---------------------------------------
  155.              5/2
  156.               8 B
  157.  
  158.  
  159.              2
  160.             C
  161.             ---
  162.             4 B             C               INF
  163.     SQRT(%PI) %E    (1 - ERF(---------))          /             2
  164.                  2 SQRT(B)      C     [        2   - B X  - C X
  165.       + ------------------------------------ - ---- = I       X  %E         dX
  166.               3/2              2   ]
  167.                4 B               4 B    /
  168.                                0
  169.  
  170. /* the a=2 case of the Laplace transform */
  171. (C12) d8,a=2;
  172.  
  173.                       INF
  174.                   2         /            2
  175.              2     - B X         [     2   - B X  - C X
  176. (D12)         LAPLACE(X  %E      , X, C) = I    X  %E            dX
  177.                      ]
  178.                      /
  179.                       0
  180.  
  181. /* two things equal to a third are equal to each other */
  182. (C13) subst(reverse(d11),%);
  183.  
  184.                            2
  185.                           C
  186.                           ---
  187.                           2      4 B           C
  188.             2       SQRT(%PI) C  %E    (1 - ERF(---------))
  189.            2   - B X                       2 SQRT(B)
  190. (D13) LAPLACE(X  %E     , X, C) = ---------------------------------------
  191.                               5/2
  192.                            8 B
  193.  
  194.  
  195.                          2
  196.                         C
  197.                         ---
  198.                         4 B         C
  199.                     SQRT(%PI) %E    (1 - ERF(---------))
  200.                                  2 SQRT(B)        C
  201.                   + ------------------------------------ - ----
  202.                               3/2              2
  203.                            4 B               4 B
  204.  
  205. /* get the right hand side */
  206. (C14) rhs(%)$
  207.  
  208. /* factor the first two terms together */
  209. (C15) factor(part(%,[1,2]))+rest(%,2);
  210.  
  211. Is  B  positive, negative, or zero?
  212. P;
  213.  
  214.                      2
  215.                     C
  216.                     ---
  217.              2        4 B         C
  218.          SQRT(%PI) (C  + 2 B) %E    (ERF(---------) - 1)
  219.                          2 SQRT(B)            C
  220. (D15)        - ----------------------------------------------- - ----
  221.                     5/2                  2
  222.                  8 B                   4 B
  223.  
  224. /* the numerical values, also confirms Alexander's Maple result */
  225. (C16) %,b=3/2,c=-1;
  226.  
  227.                  1/6           SQRT(2)
  228.          2 SQRT(2) %E    SQRT(%PI) (- ERF(---------) - 1)
  229.          1                      2 SQRT(3)
  230. (D16)          - - ------------------------------------------------
  231.          9                9 SQRT(3)
  232.  
  233. /* same as before */
  234. (C17) %,numer;
  235.  
  236. (D17)                    0.656798012459442
  237.  
  238. /* now use the hypergeometric integrator */
  239. (C18) specint(x^a*exp(-b*x^2-c*x),x);
  240.  
  241. Is  C  positive, negative, or zero?
  242. N;
  243.  
  244. /* doesn't handle this case as of yet */
  245. (D18)                 OTHER-DEFINT-TO-FOLLOW-NEGTEST
  246.  
  247. /* ok, get an answer for positive c */
  248. (C19) ''c18;
  249.  
  250. Is  C  positive, negative, or zero?
  251. P;
  252.                  2
  253.                 C
  254.                 ---
  255.        A + 1            8 B
  256. (D19) 2         GAMMA(A + 1) %E
  257.  
  258.  
  259.          - A - 1
  260.          ------- + 1/2                    2
  261.         2        1/4                   C
  262.   SQRT(%PI) 2           B    %M              (---)
  263.                   - A - 1           4 B
  264.                   ------- + 1/4, - 1/4
  265.                      2
  266.  (---------------------------------------------------------
  267.                1   - A - 1
  268.          GAMMA(- - -------) SQRT(C)
  269.                2      2
  270.  
  271.  
  272.         - A - 2
  273.         ------- + 1                   2
  274.            2         1/4              C
  275.    2 SQRT(%PI) 2        B    %M             (---)
  276.                    - A - 1          4 B     A + 1  A + 1
  277.                    ------- + 1/4, 1/4          -----  -----
  278.                       2                    2      2
  279.  - -------------------------------------------------------)/(8        B      )
  280.               - A - 1
  281.           GAMMA(- -------) SQRT(C)
  282.                  2
  283.  
  284. /* use values of a and b */
  285. (C20) %,a=2,b=3/2$
  286.  
  287. /* the answer in c in terms of Whittaker functions */
  288. (C21) factor(%);
  289.  
  290.            2
  291.           C
  292.           --          2                 2
  293.           12         C                    C
  294.         %E   (4 %M            (--) - SQRT(%PI) %M           (--))
  295.               - 5/4, 1/4 6           - 5/4, - 1/4 6
  296. (D21)       - --------------------------------------------------------
  297.                  1/4  1/4
  298.                   3 2    3      SQRT(C)
  299.  
  300. /* use the Abramowitz & Stegun formula 13.1.32 */
  301. (C22) %,%m[a,b](z):=exp(-z/2)*z^(1/2+b)*m(1/2+b-a,1+2*b,z)$
  302.  
  303. /* c is positive now */
  304. (C23) %,assume_pos:true;
  305.  
  306.                  2                        2
  307.                 C                       C
  308.             2          - --                2         - --     2
  309.         3  C    3/2     12         3  1  C           12    C
  310.      4 M(2, -, --) C    %E         SQRT(%PI) M(-, -, --) SQRT(C) %E         --
  311.         2  6                 2  2  6             12
  312.     (------------------------- - ------------------------------------) %E
  313.             3/4                      1/4
  314.            6                     6
  315. (D23) - -----------------------------------------------------------------------
  316.                      1/4  1/4
  317.                   3 2    3    SQRT(C)
  318.  
  319. /* use A&S formula 13.1.2 for the Kummer function */
  320. (C24) %,m(a,b,z):=sum(gamma(a+n)*gamma(b)*z^n/(gamma(a)*gamma(b+n)*n!),n,0,inf)$
  321.  
  322. /* use c positive */
  323. (C25) %,assume_pos:true;
  324.  
  325.                     2
  326.                    C  INF
  327.                  - -- ====   2 N
  328.                3/2       12 \        C    GAMMA(N + 2)
  329.        2  2 SQRT(%PI) C    %E      >    ------------------
  330.       C                  /         N            3
  331.       --                  ====  6  N! GAMMA(N + -)
  332.       12                  N = 0            2
  333. (D25) - %E   (------------------------------------------------
  334.                      3/4
  335.                     6
  336.  
  337.  
  338.                   2
  339.                  C  INF    2 N         3
  340.                - -- ====  C       GAMMA(N + -)
  341.                  12 \             2
  342.      2 SQRT(%PI) SQRT(C) %E      >    ------------------
  343.                 /      N          1
  344.                 ====  6  N! GAMMA(N + -)
  345.                 N = 0              2          1/4  1/4
  346.    - ---------------------------------------------------)/(3 2      3    SQRT(C))
  347.                  1/4
  348.                 6
  349.  
  350. /* factor the expression */
  351. (C26) factor(%)$
  352.  
  353. /* put terms inside of sums */
  354. (C27) intosum(multthru(%))$
  355.  
  356. /* combine into one sum, note there are odd and even terms in c */
  357. (C28) sumcontract(%);
  358.  
  359.       INF          2 N  - N - 1/4       2 N + 3
  360.       ====   2 SQRT(%PI) C    6             GAMMA(-------)
  361.       \                          2
  362. (D28)  >    (------------------------------------------
  363.       /              1/4  1/4            2 N + 1
  364.       ====       3 2      3    N! GAMMA(-------)
  365.       N = 0                   2
  366.  
  367.  
  368.                           2 N + 1  - N - 3/4
  369.                  2 SQRT(%PI) C          6             GAMMA(N + 2)
  370.                    - --------------------------------------------)
  371.                        1/4  1/4         2 N + 3
  372.                     3 2    3    N! GAMMA(-------)
  373.                                 2
  374.  
  375. /* factor the whole thing */
  376. (C29) factor(%);
  377.  
  378.            INF
  379.            ====
  380.            \      2 N  - N - 3/4
  381. (D29) 2 SQRT(%PI) ( >     C    6
  382.            /
  383.            ====
  384.            N = 0
  385.  
  386.  
  387.            2 2 N + 3             2 N + 1
  388.  (SQRT(6) GAMMA (-------) - C GAMMA(N + 2) GAMMA(-------))
  389.             2                    2
  390.  
  391.  
  392.        2 N + 1      2 N + 3     1/4  1/4
  393. /(N! GAMMA(-------) GAMMA(-------)))/(3 2    3     )
  394.           2             2
  395.  
  396. /* get the first several terms of the taylor series around c=0 */
  397. (C30) taylor(%,c,0,4);
  398.  
  399.           1/4   1/4 3   1/4 3            1/4      1/4 3      1/4 3
  400.      SQRT(6) 6    (2   )  (3   )  SQRT(%PI)      (6    (2   )  (3   ) ) C
  401. (D30)/T/ -------------------------------------- - ------------------------
  402.               108                     27
  403.  
  404.  
  405.          1/4   1/4 3   1/4 3         2         1/4   1/4 3   1/4 3   3
  406.    (SQRT(6) 6    (2   )  (3   )  SQRT(%PI)) C     (2 6    (2   )  (3   ) ) C
  407.  + ------------------------------------------- - ---------------------------
  408.                216                     243
  409.  
  410.  
  411.            1/4   1/4 3   1/4 3           4
  412.    (5 SQRT(6) 6       (2   )  (3   )  SQRT(%PI)) C
  413.  + --------------------------------------------- + . . .
  414.                7776
  415.  
  416.  
  417. /* simplify each term with radcan because of the radicals */
  418. (C31) map('radcan,%);
  419.  
  420.                    4      3                 2
  421.       5 SQRT(2) SQRT(3) SQRT(%PI) C    4 C    SQRT(2) SQRT(3) SQRT(%PI) C
  422. (D31) ------------------------------ - ---- + ----------------------------
  423.            1296                81           36
  424.  
  425.  
  426.                         2 C   SQRT(2) SQRT(3) SQRT(%PI)
  427.                           - --- + -------------------------
  428.                          9         18
  429.  
  430. /* the direct laplace transform in question */
  431. (C32) laplace(x^2*exp(-3/2*x^2),x,c)$
  432.  
  433. /* taylor around c=0 to 4 terms */
  434. (C33) taylor(%,c,0,4);
  435.  
  436.                                     2
  437.      SQRT(2) SQRT(3) SQRT(%PI)   2 C   (SQRT(2) SQRT(3) SQRT(%PI)) C
  438. (D33)/T/ ------------------------- - --- + ------------------------------
  439.             18              9             36
  440.  
  441.  
  442.                    3                      4
  443.                 4 C    (5 SQRT(2) SQRT(3) SQRT(%PI)) C
  444.                   - ---- + -------------------------------- + . . .
  445.                  81             1296
  446.  
  447. /* laplace transform answer is identical to the hypergeometric integrator answer
  448. through c^4, at least */
  449. (C34) %-d31;
  450.  
  451. (D34)/T/                0 + . . .
  452.  
  453. Leo Harten
  454. Paradigm Associates, Inc.
  455. 29 Putnam Avenue, Suite 6
  456. Cambridge, MA 02139/USA
  457. LPH@PARADIGM.COM
  458. +1(617)492-6079
  459.