Speed up SymPy equation solver

Question

I am trying to solve a set of equations using the following python code (using SymPy of course):

def Solve(kp1, kp2):
    a, b, d, e, f = S('a b d e f'.split())
    equations = [
      Eq(a+b, 2.6),
      Eq(2*a + b + d + 2*f, 7),
      Eq(d + e, 2),
      Eq(a*e,kp2*b*d),
      Eq( ((b * (f**0.5))/a)*((a+b+d+e+f+13.16)**-0.5), kp1)
    ]
    return solve(equations)

The code solves the equations successfully but after approximately 35 seconds. The Solve() function is used inside an iterative block(about 2000 iterations) in another file so speed really matters to me.

Is there any way to speed up the solver? If not can you recommend another way to solve system of equations using python?

Answer 1

You only need to solve the equation one time. After that you will have equations of the form:

• a = f_1(kp1, kp2)
• b = f_2(kp1, kp2)
• ...

so you can simply compute a, ..., e in dependency of kp1 and kp2. For example solving the first, third and fourth equations gives:

• b: -a + 2.6
• e: 2.0*kp2*(5.0*a - 13.0)/(5.0*a*kp2 - 5.0*a - 13.0*kp2),
• d: 10.0*a/(-5.0*a*kp2 + 5.0*a + 13.0*kp2)

Solving all five equations on my pc is too slow, but if it gives you an expression you just have to plug in (substitute) the values for kp1 and kp2, you don't have to solve the equations again. For substitution take a look at the sympy documentation.

So your loop should look like:

solutions = sympy.solve(eqs, exclude=[kp1, kp2])
for data_kp1, data_kp2 in data:
    for key, eq in solutions:
        solution = eq.subs([(kp1, data_kp1), (kp2, data_kp2)])
        final_solutions.append("{key}={solution}".format(key=key, solution=solution))

Answer 2

Solve the 1st 4 linear equations for a,b,d,e. This generates two solutions. Substitute each of these into the 5th equation (in this form: Eq(b**2*f, d**2*kp1**2*(a + b + 2*d + f + 13.16))) to give two equations (e51 and e52). These two equations are nonlinear in f and, if you use unrad on them you will end up with 2 sixth-order polynomials in f -- which is likely why you get no solution since only quartics are exactly solvable in general. Call those two equations e51u, and e52u.

If you are interested in the real roots, you can use real_roots to give the roots to these polynomials after substituting in the desired values for kp1 and kp2 to leave only f as the unknown. For example,

>>> ans = solve(equations[:-1],a,b,d,e)
>>> l=equations[-1]  # modified as suggested
>>> l=l.lhs-l.rhs
>>> l
b**2*f - d**2*kp1**2*(a + b + 2*d + f + 13.16)
>>> e51=l.subs(ans[0]); e51u=unrad(e51)[0]
>>> e52=l.subs(ans[1]); e52u=unrad(e52)[0]
>>> import random
>>> for r1,r2 in [[random.random() for i in range(2)] for j in range(3)]:
...     print [i.n(2) for i in real_roots(e51u.subs(dict(kp1=r1,kp2=r2)))]
...     print [i.n(2) for i in real_roots(e52u.subs(dict(kp1=r1,kp2=r2)))]
...     print '^_r1,r2=',r1,r2
...     
[1.7, 2.9, 3.0, 8.2]
[1.7, 2.9, 3.0, 8.2]
^_r1,r2= 0.937748743197 0.134640776315
[1.3, 2.3, 4.7, 7.4]
[1.3, 2.3, 4.7, 7.4]
^_r1,r2= 0.490002815309 0.324553144174
[1.1, 2.1]
[1.1, 2.1]
^_r1,r2= 0.308803300429 0.595356213169

It appears that e51u and e52u both give the same solutions, so perhaps you only need to use one of them. And you should check the answers in the original equations to see which one(s) are true solutions:

>>> r1,r2
(0.30880330042869408, 0.59535621316941589)
>>> [e51.subs(dict(kp1=r1,kp2=r2,f=i)).n(2) for i in real_roots(e51u.subs(dict(kp1=r1,kp2=r2)))]
[1.0e-12, 13.]

So here you see that only the first solution (seen to be 1.1 from above) is actually a solution; 2.1 was a spurious solution.