Integrating in Python using Sympy

Question

I am currently using Sympy to help me perform mathematical calculations. Right now, I am trying to perform a numerical integration, but keep getting an error any time I run the script. Here is the script:

from sympy import *
cst = { 'qe':1.60217646*10**-19, 'm0':N(1.25663706*10**-6) }

d = 3.6*10**-2
l = 20.3*10**-2
n = 217.0
I = 10.2

# Circum of loops
circ = l/n;
# Radius
r = N( circ/(2*pi) )
# Flux through a ring a distance R from the ceter
def flux(rad, I, loopRad):
    distFromWire = loopRad - rad
    bPoint = cst['m0']*I/(2*pi*distFromWire)
    return ( bPoint*2*pi*rad )

# Integrate from r=0 to r=wireRad
x = Symbol('x')
ig = Symbol('ig')
ig = flux(x, I, r)
print(ig)
integrate(ig*x,x)

I am sure there is probably something wrong with the actual physics/math, but right now I just want it to integrate. Here is the output I get when I run the script:

8.05359718208634e-5*x/(-6.28318530717959*x + 0.000935483870967742)
Traceback (most recent call last):
  File "script.py", line 34, in <module>
    integrate(ig*x,x)
  File "C:\Python27\lib\site-packages\sympy\utilities\decorator.py", line 24, in threaded_func
    return func(expr, *args, **kwargs)
  File "C:\Python27\lib\site-packages\sympy\integrals\integrals.py", line 847, in integrate
    return integral.doit(deep = False)
  File "C:\Python27\lib\site-packages\sympy\integrals\integrals.py", line 364, in doit
    antideriv = self._eval_integral(function, xab[0])
  File "C:\Python27\lib\site-packages\sympy\integrals\integrals.py", line 577, in _eval_integral
    parts.append(coeff * ratint(g, x))
  File "C:\Python27\lib\site-packages\sympy\integrals\rationaltools.py", line 42, in ratint
    g, h = ratint_ratpart(p, q, x)
  File "C:\Python27\lib\site-packages\sympy\integrals\rationaltools.py", line 124, in ratint_ratpart
    H = f - A.diff()*v + A*(u.diff()*v).quo(u) - B*u
  File "C:\Python27\lib\site-packages\sympy\core\decorators.py", line 75, in __sympifyit_wrapper
    return func(a, sympify(b, strict=True))
  File "C:\Python27\lib\site-packages\sympy\polys\polytools.py", line 3360, in __mul__
    return f.mul(g)
  File "C:\Python27\lib\site-packages\sympy\polys\polytools.py", line 1295, in mul
    _, per, F, G = f._unify(g)
  File "C:\Python27\lib\site-packages\sympy\polys\polytools.py", line 377, in _unify
    F = f.rep.convert(dom)
  File "C:\Python27\lib\site-packages\sympy\polys\polyclasses.py", line 277, in convert
    return DMP(dmp_convert(f.rep, f.lev, f.dom, dom), dom, f.lev)
  File "C:\Python27\lib\site-packages\sympy\polys\densebasic.py", line 530, in dmp_convert
    return dup_convert(f, K0, K1)
  File "C:\Python27\lib\site-packages\sympy\polys\densebasic.py", line 506, in dup_convert
    return dup_strip([ K1.convert(c, K0) for c in f ])
  File "C:\Python27\lib\site-packages\sympy\polys\domains\domain.py", line 85, in convert
    raise CoercionFailed("can't convert %s of type %s to %s" % (a, K0, K1))
sympy.polys.polyerrors.CoercionFailed: can't convert DMP([1, 0], ZZ) of type ZZ[_b1] to RR
[Finished in 0.3s with exit code 1]

EDIT: Alright, so I pulled out the numbers that the program was using and put it in wolfram alpha. Turns out that the integral doesn't converge, hence the error. I guess it WAS just a math error.

Answer 1

It is a very bad idea to use a symbolic library for numerical work. Just use scipy/numpy. That being said, for such a simple integral you could have used sympy. However you should actually use sympy expressions, not dump everything in opaque functions.

Firstly, learn how do variables in python work:

ig = Symbol('ig')
ig = flux(x, I, r)

After this operation ig is not a Symbol any more, it is just the return value of flux.

Define all your symbols and then make an expression out of them. The integral is sufficiently simple for sympy to handle it.

Finally, integral as simple as const*x/(x-const) as in your case should be done by hand, not wasted on software.

[EDIT]: I have rewritten it cleanly, and still sympy does not integrate correctly because of a bug. You could report it on the mailing list or issue tracker and they will try to correct it. That being said, the expression is so simple that it can be integrated by hand.

[EDIT2]:

In [5]: integrate(a*x/(b*x+c), x)
Out[5]: 

  ⎛         ⎛ 2        ⎞⎞
  ⎜x   c⋅log⎝b ⋅x + b⋅c⎠⎟
a⋅⎜─ - ─────────────────⎟
  ⎜b            2       ⎟
  ⎝            b        ⎠