How do I calculate integral of exp(a*x**2+b*x+c)?

Viewed 270

So i want to calculate integral of exp(-ax^2-bx-c) assuming that a>0, b>0, c>0 using sympy.

So i done this:

from sympy import *

x = symbols('x')
a, b, c = symbols('a, b, c')
assumptions = Q.real(a), Q.positive(a), \
              Q.real(b), Q.positive(b), \
              Q.real(c), Q.positive(c)

with assuming(*assumptions):
    expression = exp(-a*x**2 - b*x - c)
    ans = integrate(expression, (x, -oo, oo))
    ans = simplify(ans)

print(str(ans).replace(',', ',\n\n'))

And i get this:

Piecewise((sqrt(pi)*(erf(b/(2*sqrt(a))) + erfc(b/(2*sqrt(a))) + 1)*exp(-c + b**2/(4*a))/(2*sqrt(a)),

 (Abs(arg(a)) < pi/2) | ((Abs(arg(a)) <= pi/2) & (2*Abs(arg(b)) < pi) & (Abs(2*arg(b) + 2*pi) < pi))),

 (Integral(exp(-a*x**2 - b*x - c),

 (x,

 -oo,

 oo)),

 True))

But i know, that answer must be:

sqrt(pi/a) * exp(b**2/(4*a) - c)

As i think, assuming isnt working, so i tried this:

from sympy import *

a = symbols('a')
assumptions = [ Q.positive(a) ]
with assuming(*assumptions):
    print(ask(Q.is_true(a > 0)))

And the answer is:

None

And this is very strange. What am i doing wrong?

1 Answers

The sympy documentation is confusing about this. You are trying to use the "new" assumptions system which is still experimental and is not actually used internally by sympy. You need to use the core assumptions system in which assumptions are set when creating symbols:

In [50]: from sympy import * 
    ...:  
    ...: x = symbols('x') 
    ...: a, b, c = symbols('a, b, c', positive=True)                                                                                           

In [51]: integrate(exp(-a*x**2 - b*x - c), (x, -oo, oo))                                                                                       
Out[51]: 
  ⎛         2                     ⎞                             
  ⎜        b                     2⎟                             
  ⎜       ───                   b ⎟                 2           
  ⎜       4⋅a    ⎛ b  ⎞        ───⎟                b            
  ⎜  π⋅b⋅ℯ   ⋅erf⎜────⎟        4⋅a⎟               ───           
  ⎜              ⎝2⋅√a⎠   π⋅b⋅ℯ   ⎟  -c       -c  4⋅a     ⎛ b  ⎞
  ⎜- ────────────────── - ────────⎟⋅ℯ     √π⋅ℯ  ⋅ℯ   ⋅erfc⎜────⎟
  ⎝         2⋅√a            2⋅√a  ⎠                       ⎝2⋅√a⎠
- ───────────────────────────────────── + ──────────────────────
                   √π⋅b                            2⋅√a         

That result can be simplified. It looks like simplify doesn't handle erf plus erfc but we can give it a little help:

In [52]: integrate(exp(-a*x**2 - b*x - c), (x, -oo, oo)).simplify()                                                                            
Out[52]: 
                                        2
                                       b 
                                 -c + ───
   ⎛   ⎛ b  ⎞       ⎛ b  ⎞    ⎞       4⋅a
√π⋅⎜erf⎜────⎟ + erfc⎜────⎟ + 1⎟⋅ℯ        
   ⎝   ⎝2⋅√a⎠       ⎝2⋅√a⎠    ⎠          
─────────────────────────────────────────
                   2⋅√a                  

In [53]: integrate(exp(-a*x**2 - b*x - c), (x, -oo, oo)).simplify().rewrite(erfc)                                                              
Out[53]: 
           2
          b 
    -c + ───
         4⋅a
√π⋅ℯ        
────────────
     √a 
Related