How to calculate definite integral and get back a numerical value with python?

Question

I am doing a little complicated integral with python and I would need the result to be a numerical value, for example 2003708.58843. Now my python program gives some really weird result containing some "hyper" function etc. If someone can help me with this I would appreciate it really much.

Current program is here:

import sympy as sp

x = sp.symbols('x') 
f1 = 5000 + ((-(-7000+x)*(-3000+x))**(1/2))
f2 = 5000 - ((-(-7000+x)*(-3000+x))**(1/2))
minimum1 = 3000
maximum1 = (500/3)*(35-2*((61)**(1/2)))
int1 = sp.integrate(f1, (x,minimum1,maximum1))
int2 = sp.integrate(f2, (x,minimum1,maximum1))
minimum2 = (500/3)*(35-2*((61)**(1/2)))
maximum2 = 4500
f3 = 5000 + (1/2)*((-(-4500+x)*(-500+x))**(1/2))
f4 = 5000 - (1/2)*((-(-4500+x)*(-500+x))**(1/2))
int3 = sp.integrate(f3, (x,minimum2,maximum2))
int4 = sp.integrate(f4, (x,minimum2,maximum2))
I1 = int1-int2
I2 = int3-int4
total = I1 + I2
print(total)

Answer 1

When using sympy you should be careful about using floats e.g. 1/2 and instead use exact rationals e.g. Rational(1, 2). Changing that I get:

In [2]: import sympy as sp

In [3]: x = sp.symbols('x')
   ...: f1 = 5000 + ((-(-7000+x)*(-3000+x))**(Rational(1, 2))) 
   ...: f2 = 5000 - ((-(-7000+x)*(-3000+x))**(Rational(1, 2))) 
   ...: minimum1 = 3000 
   ...: maximum1 = (Rational(500, 3))*(35-2*((61)**(Rational(1, 2)))) 
   ...: int1 = sp.integrate(f1, (x,minimum1,maximum1))

In [4]: int1
Out[4]: 
                            _________________                      5/2                          3/2               ⎛        _________________⎞           
                           ╱ 8500   1000⋅√61      ⎛8500   1000⋅√61⎞            ⎛8500   1000⋅√61⎞                  ⎜       ╱ 8500   1000⋅√61 ⎟           
                4000000⋅  ╱  ──── - ────────      ⎜──── - ────────⎟       3000⋅⎜──── - ────────⎟                  ⎜√10⋅  ╱  ──── - ──────── ⎟           
  5000000⋅√61           ╲╱    3        3          ⎝ 3        3    ⎠            ⎝ 3        3    ⎠                  ⎜    ╲╱    3        3     ⎟   42500000
- ─────────── - ───────────────────────────── - ─────────────────────── + ───────────────────────── + 4000000⋅asin⎜─────────────────────────⎟ + ────────
       3                _________________             _________________         _________________                 ⎝           200           ⎠      3    
                       ╱ 3500   1000⋅√61             ╱ 3500   1000⋅√61         ╱ 3500   1000⋅√61                                                        
                      ╱  ──── + ────────        2⋅  ╱  ──── + ────────        ╱  ──── + ────────                                                        
                    ╲╱    3        3              ╲╱    3        3          ╲╱    3        3                                                            

If you want to evaluate that result in floating point to 15 decimal digits then you can do:

In [5]: int1.evalf()
Out[5]: 1294014.90427420

To compute the approximate numeric result directly using numerical quadrature (analogous to the answer from snatchysquid) you can do

In [6]: sp.Integral(f1, (x, minimum1, maximum1)).evalf()
Out[6]: 1294014.90427420

Answer 2

Try this:

from scipy.integrate import quad

def integrand(x):
    return x**2

ans, err = quad(integrand, 0, 1)
print ans

Simply change the integrand function to whatever you want and 1, 0 to the limits of integration you desire.

code taken from here.