home *** CD-ROM | disk | FTP | other *** search
/ NetNews Usenet Archive 1992 #16 / NN_1992_16.iso / spool / sci / math / symbolic / 2035 < prev    next >
Encoding:
Text File  |  1992-07-21  |  4.6 KB  |  131 lines
  1. Newsgroups: sci.math.symbolic
  2. Path: sparky!uunet!gatech!destroyer!news.iastate.edu!pv3448.vincent.iastate.edu!alex
  3. From: alex@iastate.edu (Roger Alexander)
  4. Subject: Gaussian moment integral Was: symbolic integral wanted
  5. Message-ID: <alex.711747852@pv3448.vincent.iastate.edu>
  6. Summary: Mma bug?
  7. Keywords: Mma, Gaussian moments, symbolic integration
  8. Sender: news@news.iastate.edu (USENET News System)
  9. Organization: Iowa State University, Ames IA
  10. Date: Tue, 21 Jul 1992 19:44:12 GMT
  11. Lines: 118
  12.  
  13. In a previous post Richard Waterman (rpw001@fisher.stat.psu.edu)
  14. describes trouble getting Mma to evaluate the integral
  15.  
  16. =>    ans1 = Integrate[x^a Exp[- b x^2 - c x],{x,0,Infinity}]
  17.  
  18. and reports
  19.  
  20. =>    Mma gives an "answer", no problem.
  21.  
  22. But there *is* a problem, as he recognizes:
  23.  
  24. =>   If I substitute in some numbers for {a,b,c}, say {2,1.5,-1},
  25. =>   into the symbolic answer I get a numerical result.
  26. =>   
  27. =>   The question is, why isn't the result the same as when I do:
  28. =>   
  29. =>   number2 = NIntegrate:x^2 Exp:- 1.5 x^2 + 1 x:,{x,0,Infinity}:
  30.  
  31. The answer is, Mma's "answer" to the first problem is wrong.  You 
  32. can check this by looking it up in Gradshteyn/Ryzhik, for example.
  33. Or let Mma display the inconsistency itself:
  34.  
  35. Mathematica 2.0 for DEC RISC
  36. Copyright 1988-91 Wolfram Research, Inc.
  37.  -- X11 windows graphics initialized -- 
  38.  
  39. In[1]:= ans1 = Integrate[x^a Exp[- b x^2 - c x],{x,0,Infinity}]
  40.  
  41.                                                             2
  42.          -1 - a/2                                    a  3  c
  43.         b         c Gamma[1 + a] HypergeometricU[1 + -, -, ---]
  44.                                                      2  2  4 b
  45. Out[1]= -------------------------------------------------------
  46.                                     a
  47.                                  4 2
  48.  
  49. (* Of course, the integral doesn't even converge unless certain
  50. conditions hold on a,b and c that Mma hasn't bothered to ask us about.
  51. An unsympathetic observer might feel that we deserve what happens
  52. next.  This post merely reports facts; the reader who wishes may
  53. assign blame as seems appropriate. *)
  54.  
  55. (* We assign values to the constants...*)
  56.  
  57. In[2]:= subs = {a -> 2, b -> 3/2, c -> -1}
  58.  
  59.                       3
  60. Out[2]= {a -> 2, b -> -, c -> -1}
  61.                       2
  62.  
  63. (* ... substitute into the above expression, evaluate numerically ... *)
  64.  
  65. In[3]:= %1 /. subs
  66.  
  67.                              3  1
  68.          -HypergeometricU[2, -, -]
  69.                              2  6
  70. Out[3]= -------------------------
  71.                     18
  72.  
  73. In[4]:= % //N
  74.  
  75. Out[4]= -0.103054
  76.  
  77. (* Thus Mma refutes its answer Out[1]: the integrand was positive.  If
  78. we substitute the constants before integrating we get a different
  79. (correct? Exercise!) answer. *)
  80.  
  81. In[5]:= Integrate[(x^a Exp[- b x^2 - c x])/.subs, {x,0,Infinity}]
  82.  
  83.                         1/6       1     1                         3  1
  84.         12 + Sqrt[384] E    Gamma[-, 0, -] + 3 HypergeometricU[2, -, -]
  85.                                   2     6                         2  6
  86. Out[5]= ---------------------------------------------------------------
  87.                                       54
  88.  
  89. In[6]:= % //N
  90.  
  91. Out[6]= 0.656798
  92.  
  93. =======================================
  94.  
  95. I tried Maple V out on this integral too:
  96.  
  97.     |\^/|      MAPLE V
  98. ._|\|   |/|_.  Copyright (c) 1981-1990 by the University of Waterloo.
  99.  \  MAPLE  /   All rights reserved.  MAPLE is a registered trademark of
  100.  <____ ____>   Waterloo Maple Software.
  101.       |        Type ? for help.
  102. > int(x^a * exp( -b*x^2 - c*x ),x=0..infinity);
  103.                         infinity
  104.                            /
  105.                           |       a          2
  106.                           |      x  exp(- b x  - c x) dx
  107.                           |
  108.                          /
  109.                          0
  110.  
  111. # Maple refuses to evaluate the integral without further information
  112. # on the constants a, b, and c.
  113.  
  114. > a := 2: b := 3/2: c := -1:
  115. ________________________________________________________________________
  116. > int(x^a * exp( -b*x^2 - c*x ),x=0..infinity);
  117.                    1/2           1/2  1/2
  118.             2/27 Pi    exp(1/6) 3    2    + 1/9
  119.  
  120.                           1/2           1/2  1/2          1/2  1/2
  121.                  + 2/27 Pi    exp(1/6) 3    2    erf(1/6 3    2   )
  122. ________________________________________________________________________
  123. > evalf(");
  124.                                   .6567979782
  125. ------------------------------------------------------------------------
  126. -- 
  127. *_________________*  "Faculty members encourage the free pursuit of       
  128.  Roger Alexander  | learning in their students and protect academic          
  129.  alex@iastate.edu | freedom of students." --ISU Faculty Handbook
  130. *_________________* [official statement of policies & procedures]
  131.