Python Integration Using Sympy

Question

I am currently using Sympy to help me do math calculations. Right now, I'm trying to do numerical integration, but still get an error message when running 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'm sure there is probably something wrong with actual physics / math, but right now I just want it to integrate. Here is the result that I get when running 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: Okay, so I pulled out the numbers the program used and put them in alpha-tungsten. It turns out that the integral does not converge, therefore, an error. I assume this is just a math error.

Answer 1

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

First, find out how variables work in python:

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

After this operation, ig no longer a symbol, it is just the return value of flux .

Define all your characters, and then make an expression out of them. The integral is simple enough to handle it.

Finally, integral, as simple as const*x/(x-const) , as in your case, should be executed manually, and not disappear from the program.

[EDIT]: I rewrote it cleanly, and yet sympy does not integrate correctly due to an error. You can report it on the mailing list or a problem and they will try to fix it. Moreover, the expression is so simple that it can be integrated manually.

[EDIT2]:

 In [5]: integrate(a*x/(b*x+c), x) Out[5]: ⎛ ⎛ 2 ⎞⎞ ⎜xc⋅log⎝b ⋅x + b⋅c⎠⎟ a⋅⎜─ - ─────────────────⎟ ⎜b 2 ⎟ ⎝ b ⎠