The result seems to be a list of dicts so you should be able to use the standard Python operations c_solutions = [sol[c] for sol in solutions] ...

Answering my own question: in the latest version of sympy 0.7.6 (was using a slightly older version previously), indexedBase objects are made to be commutative by default. So just download the newer version of sympy and i can use the Mod operator. https://github.com/sympy/sympy/pull/7355...

matrix,sympy,symbolic-math,polynomial-math

Since collect(expand(V), ...) already shows V as a linear polynomial in the variables v_1, theta_1, v_2, theta_2, instead of using V.match(form), perhaps an easier, more direct way to get the coefficients is to use the V.coeff method: N = sy.Matrix([V.coeff(v) for v in (v_1, theta_1, v_2, theta_2)]).transpose() import sympy as...

TL; DR It is a bug. If you don't believe it, try this: In [1]: from sympy import factor, Symbol In [2]: factor(1e-20*Symbol('t')-7.292115e-5) Out[2]: -2785579325.00000 Two years ago, the default value for the parameter tol in RealField.__init__ was changed from None to False in commit polys: Disabled automatic reduction to...

You normally get first the parameter and then evaluate to obtain the function value. For example: from sympy import * x = Symbol('x', real=True) # parameter f = -2 * x**2 + 4*x # function derivative = f.diff(x) # -4*x + 4 solve(derivative, x) # -4*x + 4 = 0...

You can use the first construct with splatting: X = subs(X,zip(toin,toout)...) ...

Use the rational=False flag: >>> print filldedent(solve([eq1,eq2,eq3],[b,n,k], rational=False)) [(-3.44827586206897*log((-0.000335859591913345*p + 0.00110833665331404)/(2.24927535168052e-5*p - 0.000139455071804192)) + 2, 0.689655172413793*p*log((-0.000335859591913345*p + 0.00110833665331404)/(2.24927535168052e-5*p - 0.000139455071804192)), (0.000335859591913345*((-0.000335859591913345*p + 0.00110833665331404)/(2.24927535168052e-5*p - 0.000139455071804192))**1.13793103448276 +...

Use the 'ordered' function: >>> list(ordered((g(x),f(x)))) [f(x), g(x)] It returns a generator -- that's why the result is wrapped with list....

If you never want the denominator pulled out, set the long_frac_ratio to oo. The ee defined below has the same symptoms as your expression; notice that setting the ratio forces it into simple over/under format: >>> ee (-2*x*y*z + x)/(x - 1) >>> latex(ee) '\\frac{1}{x - 1} \\left(- 2 x...

When you're doing advanced expression manipulation like this (okay, your example is still simple), the replace method is very useful: test.replace(f, lambda arg: 2*arg) # returns: 2*x + 2*y From the docs: [replace] Traverses an expression tree and performs replacement of matching subexpressions from the bottom to the top of...

You can't use mpmath functions with symbolic SymPy objects. You need to define sinc symbolically. One way is to just define it as a Python function that returns the equivalent (here sin is sympy.sin): def sinc(x): return sin(x)/x Or, you could write your own SymPy class for it. class sinc(Function):...

python,python-3.x,sympy,polynomial-math,gmpy

There are a few obstacles: gmpy.mpfr has no ._sympy_ method, so it will be converted to float in an intermediate step. sympy.Poly, by default, uses sympy.polys.domain.RR for floating point coefficients, and will use it to convert. RR, when it loads, uses 53 as its precision, and ignores mpmath.mp.precision. When converting,...

For simple display replacement you could use: def format_math(string): return (string.replace("**", "^")).replace("*", "") Then you could use it versus user input to compare their input answer versus yours. x = format_math("-25*x**(3/5)/3 + 6*x**(4/3) - 5*x**6/3 + x**2/2 - 4*x") # -25x^(3/5)/3 + 6x^(4/3) - 5x^6/3 + x^2/2 - 4x user_input...

Except for some shortcuts in notation, subtracting the terms you've already put into A is exactly what you should, and must, do. Here's a slightly shorter version: import sympy # Define constants R1=R3=R5 = 1e3 R2=R4=R6 = 2e3 C1 = 1e-6 C2 = 0.5e-6 x_plus = 3 x_minus = 0...

Make a distinction between the symbols and the expressions: import sympy as sy a, b, c = sy.symbols('a:c') a_expr = sy.sympify("b + 5") b_expr = sy.sympify("c + 5") c_expr = sy.sympify("5") sy.solve([sy.Eq(a,a_expr), sy.Eq(b,b_expr), sy.Eq(c,c_expr)]) yields {c: 5, b: 10, a: 15} Above, a is a Symbol, but a_expr is an...

The problem is operator precedence in Python/SymPy: you need to surround the inequalities with parentheses otherwise 1 | 2 is evaluated first and an error is raised by SymPy. You can write: >>> sympify('(2>1) | (2<1)') True # SymPy bool Of course, you don't really need SymPy's power for logical...

The function you are looking for to generate python code is python. Although it generates python code, that code will need some tweaking to remove dependence on SymPy objects as Oliver W pointed out. >>> print python( Matrix([[x**2,exp(y) + x]]).jacobian([x, y]) ) x = Symbol('x') y = Symbol('y') e =...

python,numpy,sympy,symbolic-math,polynomial-math

In SymPy this would be simply: from sympy.abc import x f = x**2 + 3*x + 2 g = f.subs({x:2/x}) Resulting in: print(g) #2 + 6/x + 4/x**2 ...

You can use ordered to put the indices in order: >>> from sympy import * >>> i, j = symbols('i j', cls=Wild) >>> x = IndexedBase('x') >>> e = x[1]*x[3] + x[2]*x[1] + x[3]**2 >>> def new(o, x): ... if o.is_Mul: ... i,j=list(ordered([i.args[1] for i in o.args])) ... elif o.is_Pow:...

python,matrix,sympy,symbolic-math

My first answer used phi.match(form) to find the coefficients, but that doesn't seem to work so well when matching many Wild symbols. So instead, I think a better approach is to use phi = collect(expand(...)) and then use phi.coeff to find the coefficients: import sympy as sy phi_1, phi_2, x,...

Unpacking your tuple will help you with the first step: >>> v=var('x:10') >>> f=Function('f') >>> f(*v) f(x0, x1, x2, x3, x4, x5, x6, x7, x8, x9) ...

This appears to be a bug. You can get around it by calling diff directly on the Poly instance. Ideally calling the function diff from the top level sympy module should yield the same result as calling the method diff. In [1]: from sympy import * In [2]: from sympy.abc...

python,sympy,polynomials,coefficients

I don't know why sympy truncates small coefficients when constructing a polynomial over reals, but it doesn't do this over rationals. So as a workaround, you can construct a polynomial with domain='QQ', extract coefficients, and then convert back to floats. Example using your polynomial: import sympy z_s = symbols('z_s') f...

Found out where the weird part is coming from: stringpict.py line 417: if num.binding==prettyForm.NEG: num = num.right(" ")[0] This is being done ONLY for the numerator. It adds a space after the numerator if the numerator is negative… Weird! I'm not sure if there can be a fixed other than...

The only thing close to this that I found in SymPy is is_zero_dimensional, which checks if an Ideal generated by a Groebner is zero dimensional. It's possible that I missed something, though, particularly in the AGCA module.

Yes, sympy 0.7.6 is wrong in this case, and in some other such cases. Generally, I don't know any symbolic math package that I would trust to do calculus with piecewise defined functions. Note that although inner = Piecewise((0, (x>=y)|(x<0)), (a,True)) throws a TypeError at integration time, a logically equivalent...

The following will give you what you ask for >>> s=R(u(x,y)).diff(x) >>> s.replace(Subs, lambda a,b,c: a.xreplace(dict(zip(b,c)))) (It will revert to a Subs instance if you apply the doit method.)...

IIUC, you might be able to combine as_coeff_exponent and as_ordered_terms, and get away with [[float(x) for x in term.as_coeff_exponent(z)] for term in G_F.as_ordered_terms()] which gives me >>> terms = [[float(x) for x in term.as_coeff_exponent(z)] for term in G_F.as_ordered_terms()] >>> terms.sort(key=lambda x: x[1]) >>> terms [[0.001953125, -1.75], [0.01171875, -1.5], [0.013671875, -1.25],...

python,sympy,units-of-measurement,quantum

The documentation says that the energy is returned in the units of hw (that is h.omega), so you use whichever units you want (they are just rescalings). That is, hw is J if you take h(bar) = 1.055e-34 J.s and w in s-1 and use SI units for everything else...

SymPy expressions are not mutable. If you want to change them, you'll need to use a function like subs to create a new expression (like z2 = z.subs(x[i + 1], 2)).

ipython-notebook,sympy,mathjax

I won't save with the notebook (you'll probably need to execute some Javascript in a cell to do that), but if you right click on an equation there is an option to scale it up or down (this is a MathJax feature, which has nothing to do with the IPython...

python,math,sympy,equation-solving

The docstring of nsolve warns that the numerator of your expression is used. In your case this is catastrophic and it is advisable to just use mpmath.findroot directly: >>> m = 10**4; c = 4 >>> f=lambda n,m,c:( 2**(n - 2)*n*(n - 1)*binomial(m, n)*factorial(m - S(3)/2)* ... factorial(m - n)/(binomial(2*m,...

python,scipy,sympy,symbolic-math

subs calls sympify on it arguments. See https://github.com/sympy/sympy/blob/master/sympy/core/basic.py#L845 for example. This means that strings will be converted to symbols but sympify has no way of knowing that you want your strings to be converted to symbols with non-default assumptions. The only way to do that is to actually pass in...

You can use: f.free_symbols which will return a set of all free symbols. Example: >>> import sympy >>> x, y, z = sympy.symbols('x:z') >>> f = sympy.exp(x + y) - sympy.sqrt(z) >>> f.free_symbols set([x, z, y]) ...

I solved with sympy: from sympy import * x, y = symbols("x y") f = (x ** 2 + y ** 2) res = integrate(f, (y, 20, x-2), (x, 22, 30)) Basically sympy.integrate is able to deal with multiple integrations, also with variable boundaries....

python,performance,numpy,lambda,sympy

It's hard to make any concrete timings without the actual equations, but there are a few recommendations regarding your code. First, let's talk about the equations: If there's always a column of zeros and ones, why bother evaluating? There seems to be symmetry in the equations: your symEqns[0][0] == symEqns[1][3]....

python,numpy,scipy,linear-algebra,sympy

Your equation represents a system of equations, where each element of v0 is expressed as a sum of the respective elements in the arrays v1,v2,v3,v4,v5. This is a perfectly determined case, i.e. the number of unknowns a1,a2,a3,s4,s5 equals the number of equations, which is the length of the vectors v1,v2,v3,v4,v5....

This is happening because of less information. I think you want to do this: In [7]: x = Symbol('x', real=True) In [8]: (log(exp(exp(x)))).simplify() Out[8]: exp(x) ...

The following should get you headed in the right direction: >>> from sympy.tensor import IndexedBase, Indexed >>> from sympy import sift >>> x = IndexedBase('x') >>> y = IndexedBase('y') >>> e = x[1]* y[2] + x[5]* y[10] >>> e.atoms(IndexedBase) set([y, x]) >>> e.atoms(Indexed) set([x[5], y[10], x[1], y[2]]) >>> sifted =...

>>> for k,v in z.items(): ... print '%s, %s' % (k.n(3, chop=True), v) ... 0.178 + 1.13*I, 1 0.178 - 1.13*I, 1 2.64, 1 ...

python,python-2.7,numpy,scipy,sympy

The problem is that your integral has no (or has a hard one) analytical solution, and therefore SymPy is returning the unevaluated integral expression. If you want the numerical value as an answer, why not use scipy.integrate.quad, for example: from scipy.integrate import quad from numpy import * from sympy import...

Constructing a sy.Poly would take a lot less time: using_loop : 43.56 using_P2 : 12.64 using_poly : 0.03 from timeit import default_timer as tc import numpy as np import sympy as sy from numpy.polynomial.polynomial import polyval2d as P2 def using_poly(coefficient_matrix, S=sy.S): order = coefficient_matrix.shape[0] x = sy.Symbol('x') y = sy.Symbol('y')...

python,math,numpy,sympy,callable-object

Looks to me like you're simply missing a *: for a,b,c,d,e,f,g,h in zip(uAFOURIERL,IAFOURIERL,IBFOURIERL,ICFOURIERL,uAFOURIERR,IAFOURIERR,IBFOURIERR,ICFOURIERR): m.append(solve(a-(x*(Rl+Xl*1J)*b+x*(Rr+Xr*1J)*c+x*(Rr+Xr*1J)*d)- (e-((1-x)*(Rl+Xl*1J)*f+(1-x)*(Rr+Xr*1J)*g+(1-x)*(Rr+Xr*1J)*h))*math.e**(alpha*1J))) ^ and turning a multiplication of (1-x) by (Rl+Xl*1J) into a function call....

Are you sure that you a running the same script that you've given us here? I say this because I can copy and paste your example verbatim and it works with absolutely no issue. Once you've checked, could you post which version of Python, SymPy, NumPy and Matplotlib you're using...

plot,ipython,sympy,polynomials

Welcome to SO! The subs() method is not meant to be used with numpy arrays. lambdify() does what you want. Try: import numpy as np import matplotlib.pyplot as plt import sympy as sy sy.init_printing() # nice formula rendering in IPython x = sy.symbols("x", real=True) # the sample polynomial: pp =...

Using sympy this can be done as follows... x = symbols('x') // user_input and original_eq are strings expr1 = sympify(user_input) // user input expr2 = sympify(original_eq) // Answer print expr1==expr2 Double equals signs (==) are used to test equality. However, this tests expressions exactly, not symbolically http://docs.sympy.org/dev/gotchas.html#double-equals-signs But this returns...

First a comment about your code: since you define q to be the sum a+b+c+d you will never see the change, even if it did work (but it doesn't). Something else that does work is the following: >>> p = a + b + 2*c + d >>> q =...

Regarding implementing it in SymPy, there is a guide on how to implement special functions here. We would love pull requests for well-known special functions. For numerical routines, they are implemented in mpmath. It looks like it just uses the definition directly. ...

Your a and b does not represent similar objects, actually a is a 1x3 "matrix" (one row, 3 columns), namely a vector, while b is a 3x1 matrix (3 rows, one column). >>> a array([1, 2, 3]) >>> b Matrix([ [1], [2], [3]]) The numpy equivalent would be numpy.array([[1], [2],...

algorithm,sympy,polynomial-math,polynomials

This is pretty well studied, as mcdowella's comment indicates. Here is how the Cantor-Zassenhaus random algorithm works for the case where you want to find the roots of a polynomial, instead of the more general factorization. Note that in the ring of polynomials with coefficients mod p, the product x(x-1)(x-2)...(x-p+1)...

If you have SSH access to the server you can install it with pip install sympy or you can use Anaconda to install it. SymPy can be used with Django, Flask or any python web framework. You can also use cgi to run python scripts.

python,python-2.7,ipython,ipython-notebook,sympy

Given that you are new to Python I would advise that you install a distribution that already includes the complete scientific python stack such as WinPython or Anaconda. If it is specifically sympy you are after you can play around online at Sympy live. If you want to stick to...

You don't need the parentheses here, .rhs is an attribute. You can write: rhs = sol.rhs sol.rhs will return a SymPy Float object....

python,python-2.x,sympy,integral

You may invoke the doit method: >>> f1.doit() f(t1) I believe SymPy is reluctant to perform these operations automatically since they may be arbitrarily expensive, and there is no universal system of predicting how expensive they will be. But maybe it would be wise to add some heuristics for integrals...

This is likely due to floating point round-off. Simplifying gives: In [10]: sympy.simplify(Lagrange(Lx,Ly)) Out[10]: X*(1.85037170770859e-17*X**2 + 1.0*X - 1.11022302462516e-16) Which is basically X**2. Try getting rid of those float casts: def Lagrange (Lx, Ly): X=sympy.symbols('X') if len(Lx)!= len(Ly): print "ERROR" return 1 y=0 for k in range ( len(Lx) ):...

python,matplotlib,plot,3d,sympy

I am turning my comment into an answer. I suggest to use mayavi and contour3d for this task. You can always rewrite your implicit function to be f(x,y,z)=0. For a sphere we have x^2 + y^2 + z^2 = r^2, that can be rewritten as f(x,y,z) = x^2 + y^2...

Your code wasn't working for me, but I found using old_poly_ring instead of poly_ring did the trick: from sympy import * class AlgebraicSet(object): def __init__(self, polynomials,field): self.polynomials = polynomials self.field = field def get_symbols(self): symbols = set() for f in self.polynomials: symbols = set(symbols | f.atoms(Symbol)) return symbols def get_coordinate_ring(self):...

python,sympy,nonlinear-functions

Although you had the right sign of the difference of r in f1 and f3, if you write in terms of squares (where the sign no longer matters) as is already done in f2 you obtain 2 real answers: >>> f1=(x0-x1)**2+(y0-y1)**2-(r0 - r1)**2 >>> f2=(zqua-x0)**2+(yqua-y0)**2-r0**2 >>> f3=(x0-x3)**2+(y0-rm)**2-(r3 - r0)**2 >>>...

You can manually adjust the long_frac_ratio long_frac_ratio: The allowed ratio of the width of the numerator to the width of the denominator before we start breaking off long fractions. The default value is 2. >>> latex(e,long_frac_ratio=3) '\\frac{2 + 3}{7}' ...

Same @smichr here with basically the same solution: def cSum(s): con, dep = factor_terms(s.function.as_independent(*s.variables)) return con*Sum(dep, *s.args[1:]) var('a b x N i') eq = Sum(2*a*x(i)**2, (i, 1, N)) + Sum(2*b*x(i), (i, 1, N)) + \ Sum(-2*x(i)**2, (i, 1, N)) >>> pprint(eq.replace(lambda s: isinstance(s, Sum), lambda s: cSum(s))) N N N...

The problem comes from using the from pylab import * line right after from sympy import * You've created a namespace collision, because solve exists both under pylab (where it is in fact an alias for numpy.linalg.solve) as well as sympy. You should try to avoid such generic imports. For...

What about this: isinstance(z, tuple(sympy.core.all_classes)) ...

python,scipy,sympy,numerical-methods,differential-equations

Here and here are some examples. As for your problem, you can write your equation like: y' + p(t)y - q(t) = 0 and then use dsolve(). import sympy t = sympy.Symbol('t') y = sympy.Function('y')(t) p = sympy.Function('p')(t) q = sympy.Function('q')(t) y_ = sympy.Derivative(y, t) # y' + p(t)y -...

If you give your symbols a specific order, as in the code below, you could convert the expression to a polynomial and get its coefficients: >>> from sympy import * >>> x, y, z, t = symbols('x y z t') >>> a1, a2, a3, a4 = symbols('a[1], a[2], a[3], a[4]')...

This behaviour is deeply engraved into sympy and you won't be able to override it by making your own class of Integral. From the Sympy docs: Finally, one last note. You may have noticed that the order we entered our expression and the order that it came out from srepr...

You're approaching it wrong. From the documentation The coefficients of Legendre polynomials can be recovered using degree-n Taylor expansion: >>> for n in range(5): ... nprint(chop(taylor(lambda x: legendre(n, x), 0, n))) ... [1.0] [0.0, 1.0] [-0.5, 0.0, 1.5] [0.0, -1.5, 0.0, 2.5] [0.375, 0.0, -3.75, 0.0, 4.375] ...

The full "fu" method tries many different combinations of transformations to find "the best" result. The individual transforms used in the Fu-routines can be used to do targeted transformations. You will have to read the documentation to learn what the different functions do, but just running through the functions of...

I'm assuming you really mean -(-1)**Rational(1, 3) + (-1)**Rational(2, 3), as literally -(-1)**(1/3) + (-1)**(2/3) is all Python (no SymPy), and evaluates numerically. Most SymPy objects do not do any kind of nontrivial simplification automatically. The reason is that sometimes you might want to represent -(-1)**(1/3) + (-1)**(2/3) without it...

python,geometry,networkx,sympy

I found a solution that sped the process up by about 13x (for a polygon with 35 points (like the data listed above), the old method from the code in the question took about 4hours to find all line segments inside the polygon. This new method took 18 minutes instead.)...

I guess those 900MB are not dependencies, but recommendations. $ apt-cache show python-sympy Package: python-sympy Priority: optional Section: universe/python Installed-Size: 14889 Maintainer: Ubuntu Developers <[email protected]> Original-Maintainer: Georges Khaznadar <[email protected]> Architecture: all Source: sympy Version: 0.7.4.1-1 Depends: python (>= 2.7), python (<< 2.8), python:any (>= 2.7.1-0ubuntu2) Recommends: python-imaging, python-ctypes, ipython, python-numpy,...

In version 3 (and earlier) I get the following >>> import sympy as sp >>> M = sp.Matrix(3,3,lambda i,j: i+j) >>> V = sp.Matrix.zeros(3, 1) >>> M.col_insert(1,V) Matrix([ [0, 0, 1, 2], [1, 0, 2, 3], [2, 0, 3, 4]]) >>> print(M) Matrix([[0, 1, 2], [1, 2, 3], [2, 3,...

Thanks to @smichr for reminding me! This question has been answered on the GitHub issue page.

I think you should use sp.log instead of np.log. I run your code but it seems like the equation is too complex and no algorithms are implemented to solve it.

>>> x in a.free_symbols, y in a.free_symbols (True, True) >>> x in b.free_symbols, y in b.free_symbols (False, True) ...

Use re.search instead of re.match because match tries to match from the beginning of the input string. >>> import re >>> s = '3(18(7x-3)+x) = 90' >>> re.search(r'(?<=\()[^()]*(?=\))', s).group() '7x-3' ...

This is a nice example of using an intelligent initial guess. If you just provide a tolerance you can find the additional solution: >>> len([nsolve(cos(x)*cosh(x)-1,x,3.14/2+3.14*n,tol=1e-12) for n in range(10)]) 10 Note, however, that the function is very steep in the region of the roots and it is unlikely that you...

You asigned values to all the variables that are on the rhs of x so when you show x you see the value that it took on with the variables that you defined. Rather than input values, why not try solve the ode symbolically if possible? >>> from sympy import...

javascript,python,ajax,django,sympy

You can include the SymPy library into your views.py file and then create a new url that acts like an API. So when you pass the data to the API it runs through SymPy then prints out the value to the page. You can use jQuery to send a POST...

Try real = True and positive = True (lower case): import sympy as sp d, R, c = sp.symbols('d R c', positive = True, real = True) dt = sp.symbols('\Delta_t', real = True) dt = (1/c**2)*(-R*c+sp.sqrt(c**2*(R+d)**2)) print sp.simplify(dt) Output: d/c ...

python,sympy,computer-algebra-systems

It seems that function composition works as you would expect in sympy: import sympy h = sympy.cos('x') g = sympy.sin(h) g Out[245]: sin(cos(x)) Or if you prefer from sympy.abc import x,y g = sympy.sin('y') f = g.subs({'y':h}) Then you can just call diff to get your derivative. g.diff() Out[246]: -sin(x)*cos(cos(x))...

Use the star expressions In [13]: g(*v) Out[13]: g(x1, x2, x3, x4) In [14]: g(*v).diff(x1) Out[14]: Derivative(g(x1, x2, x3, x4), x1) ...

Turns out, this question has a very obvious answer: MatrixSymbols in sympy can be indexed like a matrix, i.e.: X[i,j] gives the element-wise equations. If one wants to subset more than one element, the MatrixSymbol must first be converted to a sympy.Matrix class: X = sympy.Matrix(X) X # lists all...

Somewhat surprisingly, I couldn't find a way to do this "out of the box". You can use the method in this question to make any one variable the sole subject of the left hand side (LHS) of the inequality but I can't seem to make the constant term the subject....

python,numpy,combinations,sympy,itertools

Assuming I understand your aim (not sure what you might want to happen in cases of duplicate elements, namely whether you want to consider them distinct or not), you could use itertools.product: import itertools def everywhere(seq, width): for locs in itertools.product(range(width), repeat=len(seq)): output = [[] for _ in range(width)] for...

I suspect that may not be a feature in sympy for these reasons: First, if theta is not algebraic over the integers, then adjoining theta to a polynomial ring over the integers, is isomorphic. For example, pi is not algebraic over the integers, because there are no integer coefficients that,...

This document contains the items that are identified as LaTeX entities.

expr.lhs and expr.rhs will give you the left- and right-hand sides of the equation.

Very basic indeed. Multiplication in python is achieved with the * character. Hence, you should set string as follows: string = "10*x+4=7" If you'd like to auto-insert the * character in front of every x in your equations, you can define a function like this one to do the job:...

Maybe the definition of x1 and x2 had something wrong. I get the desired result using what was shown: >>> a = (k*x_2 + m)/(x_2 + 1) + (k*x_1 + m)/(x_1 + 1) >>> b = collect(collect(cancel(a), m), k) >>> b.subs({x_1 + x_2: -(2*k*m - 8)/k**2, x_1*x_2: m**2/k**2}).simplify() 8*(k +...

my previous answer was based on long time ago experiance this is the correct proccess (that worked for me based on the latest anaconda on a windows machine) I have updated the meta.yaml only with the following changes: changing to get the source from git including mpmath in build run...

A more targeted approach would use replace: Target only Numbers >>> eq.replace( ... lambda x: x.is_Symbol and S(x.name).is_Number, ... lambda x: S(x.name)) x + 1 >>> _.atoms(Number) set([1]) Target any number expression >>> Symbol('1+2').replace( ... lambda x: x.is_Symbol and S(x.name).is_Number, ... lambda x: S(x.name)) 3 Target any expression >>> Symbol('1+x').replace(...

The symbol x_i is not related to the symbol i. It prints with a subscript, but that is just a printing convention. A better choice is IndexedBase. There's a fix from https://github.com/sympy/sympy/pull/8065 that makes this easier, but for now, you can do something like this: In [13]: x = IndexedBase('x')...