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: 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        ⎠
    

