Correct sequence of commands for symbolic equation using sympy

Question

Note: I am brand new to sympy and trying to figure it out how it works.

What I have now: I do get the correct solutions but it takes 35 - 50 seconds.

Goal: To speed it up the calculations by defining symbolic equation once and then reusing it with different variables.

Set up: I need to calculate a polynomial G(t) (t = 6 roots) for every iteration of the loop. (220 iterations total) G(t) have 6 other variables, which are calculated and are known on every iterations. These variables are different on every iteration.

First try (slow): I simply put every into one python function, where I defined Gt symbolically and solved for t. It was running around 35 - 40 seconds. function_G is called on every iteration.

def function_G(f1, f2, a, b, c, d):
    t = sp.symbols('t')
    left = t * ((a * t + b)**2 + f2**2 * (c*t+d)**2)**2 
    right = (a*d-b*c) * (1+ f1**2 * t**2)**2 * (a*t+b) * (c*t+d)
    eq = sp.expand(left - right)
    roots = sp.solveset(Gt, t)
    return roots

• Then a person gave me a hint that:

You should only need to (symbolically) solve for the coefficients of the polynomial once, as a preprocessing step. After that, when processing each iterations, you simply calculate the polynomial coefficients, then solve for the roots.

• I asked for clarification the person added:

So I defined the function g(t) and then used sympy.expand to work out all parenthesis/exponents, and then sympy.collect to collect terms by powers of t. Finally I used .coeff on the output of collect to get the coefficients to feed into numpy.root.

Second try: To followed the advice, I defined a G(t) symbolically first and passed it to the function that runs the loop along with its symbolic parameters. Function constructGt() thus is called only once.

def constructGt():
    t, a, b, c, d, f1, f2 = sp.symbols('t a b c d f1 f2')
    left = t * ((a * t + b)**2 + f2**2 * (c*t+d)**2)**2 
    right = (a*d-b*c) * (1+ f1**2 * t**2)**2 * (a*t+b) * (c*t+d)
    gt = sp.Eq(left - right, 0)
    expanded = sp.expand(gt)
    expanded = sp.collect(expanded, t)

    g_vars = {
      "a": a,
      "b": b,
      "c": c,
      "d": d,
      "f1": f1,
      "f2": f2
    }

    return expanded, g_vars

then on every iteration I was passing the function and its parameters to get the roots:

#Variables values:
#a = 0.00011713490404073987 
#b = 0.00020253296124588926 
#c = 4.235688216068313e-07 
#d = 0.012262546040805029 
#f1= -0.012553203944721956 
#f2 = 0.018529776776949003

def function_G(f1_, f2_, a_, b_, c_, d_, Gt, v):
    Gt = Gt.subs([(v['a'], a_), (v['b'], b_), 
                  (v['c'], c_), (v['d'], d_),
                  (v['f1'], f1_), (v['f2'], f2_)])

    roots = sp.solveset(Gt, t)
    return roots 

But it got even slower around 56 seconds.

Question: I do not understand what Am I doing wrong? I also do not understand how this person is using .coeff() and then np.roots on the results.

Answer 1

Even if your f1 and f2 are linear in a variable, you are working with a quartic polynomial and the roots of that are very long. If this is univariate then it would be better to just use the expression, solve it at some value of constants where the solution is known and then use that value and new constants that are relatively close to the old ones and use nsolve to get the next root. If you are interested in more than one solution you may have to "follow" each root separately with nsolve...but I think you will be much happier with the overall performance. Using real_roots is another option, especially if the expression is simply a polynomial in some variable.

Given that you are working with a quartic you should keep this in mind: the general solution is so long and complicated (except for very special cases) that it is not efficient to work with the general solution and substitute in values as they are known. It is very easy to solve for numerical values, however, and it is much faster:

First create the symbolic expression into which values will be substituted; assumptionless "vanilla" symbols are used:

t, a, b, c, d, f1, f2 = symbols('t a b c d f1 f2')
left = t * ((a * t + b)**2 + f2**2 * (c*t+d)**2)**2 
right = (a*d-b*c) * (1+ f1**2 * t**2)**2 * (a*t+b) * (c*t+d)
eq = left - right

Next, define a dictionary of replacements to subtitute into the expression noting that dict(x=1) creates {'x': 1} and when this is used with subs a vanilla Symbol will be created for "x":

reps = dict(
a = 0.00011713490404073987 ,
b = 0.00020253296124588926 ,
c = 4.235688216068313e-07 ,
d = 0.012262546040805029 ,
f1= -0.012553203944721956 ,
f2 = 0.018529776776949003)

Evaluate the real roots of the expression:

from time import time
t=time();[i.n(3) for i in real_roots(eq.subs(reps))];'%s sec' % round(time()-t)
[-11.5, -1.73, 8.86, 1.06e+8]
'3 sec'

Find all 6 roots of the expression but take only the real parts:

>>> roots(eq.subs(reps))
{-11.4594523988215: 1, -1.73129179415963: 1, 8.85927293271708: 1, 106354884.4365
42: 1, -1.29328524826433 - 10.3034942999005*I: 1, -1.29328524826433 + 10.3034942
999005*I: 1}
>>> [re(i).n(3) for i in _]
[-11.5, -1.73, 8.86, 1.06e+8, -1.29, -1.29]

Change one or more values and do it again

reps.update(dict(a=2))
[i.n(3) for i in real_roots(eq.subs(reps))]
[-0.0784, -0.000101, 0.0782, 3.10e+16]

Update values in a loop:

>>> a = 1
>>> for i in range(3):
...     a += 1
...     reps.update(dict(a=a))
...     a, real_roots(eq.subs(reps))[0].n(3)
...
(2, -0.0784)
(3, -0.0640)
(4, -0.0554)

Note: when using roots, the real roots will come first in sorted order and then imaginary roots will come in conjugate pairs (but otherwise not in any given order).