How to reorganize a mathematical expression into its polynomial form

Question

I'm trying to reorganize the expression polin:

polin=1+K*(1+taua*s)/((1+tau1*s)*(1+tau2*s))*(1-s*theta/2)/(1+s*theta/2)*Kc*(1+1/(s*tauI)+tauD*s)

to this form: (a)*s^4+ (b)*s^3+ (c)*s^2 + (d)*s + (e)

I need to find the values of a, b, c, d, and e and it's complicated to do it by hand.

This is the code I used

import numpy as np
from sympy import *
Kc, K, tauI, taua, theta, tauD, tau1,tau2, s = sympy.symbols('Kc K tauI taua theta tauD tau1 tau2 s')

polin=1+K*(1+taua*s)/((1+tau1*s)*(1+tau2*s))*(1-s*theta/2)/(1+s*theta/2)*Kc*(1+1/(s*tauI)+tauD*s)

pol=poly_from_expr(polin,gen=s)
  
print(pol)

This is the output

(Poly(-1/2*s**3*(1/(s**3*tau1*tau2*theta/2 + s**2*tau1*tau2 + s**2*tau1*theta/2 + s**2*tau2*theta/2 + s*tau1 + s*tau2 + s*theta/2 + 1))*K*Kc*tauD*taua*theta + s**2*(1/(s**3*tau1*tau2*theta/2 + s**2*tau1*tau2 + s**2*tau1*theta/2 + s**2*tau2*theta/2 + s*tau1 + s*tau2 + s*theta/2 + 1))*K*Kc*tauD*taua - 1/2*s**2*(1/(s**3*tau1*tau2*theta/2 + s**2*tau1*tau2 + s**2*tau1*theta/2 + s**2*tau2*theta/2 + s*tau1 + s*tau2 + s*theta/2 + 1))*K*Kc*tauD*theta - 1/2*s**2*(1/(s**3*tau1*tau2*theta/2 + s**2*tau1*tau2 + s**2*tau1*theta/2 + s**2*tau2*theta/2 + s*tau1 + s*tau2 + s*theta/2 + 1))*K*Kc*taua*theta - s*(1/(s**3*tau1*tau2*tauI*theta + 2*s**2*tau1*tau2*tauI + s**2*tau1*tauI*theta + s**2*tau2*tauI*theta + 2*s*tau1*tauI + 2*s*tau2*tauI + s*tauI*theta + 2*tauI))*K*Kc*taua*theta + s*(1/(s**3*tau1*tau2*theta/2 + s**2*tau1*tau2 + s**2*tau1*theta/2 + s**2*tau2*theta/2 + s*tau1 + s*tau2 + s*theta/2 + 1))*K*Kc*tauD + s*(1/(s**3*tau1*tau2*theta/2 + s**2*tau1*tau2 + s**2*tau1*theta/2 + s**2*tau2*theta/2 + s*tau1 + s*tau2 + s*theta/2 + 1))*K*Kc*taua - 1/2*s*(1/(s**3*tau1*tau2*theta/2 + s**2*tau1*tau2 + s**2*tau1*theta/2 + s**2*tau2*theta/2 + s*tau1 + s*tau2 + s*theta/2 + 1))*K*Kc*theta - (1/(s**3*tau1*tau2*tauI*theta + 2*s**2*tau1*tau2*tauI + s**2*tau1*tauI*theta + s**2*tau2*tauI*theta + 2*s*tau1*tauI + 2*s*tau2*tauI + s*tauI*theta + 2*tauI))*K*Kc*theta + (1/(s**3*tau1*tau2*tauI*theta/2 + s**2*tau1*tau2*tauI + s**2*tau1*tauI*theta/2 + s**2*tau2*tauI*theta/2 + s*tau1*tauI + s*tau2*tauI + s*tauI*theta/2 + tauI))*K*Kc*taua + (1/(s**3*tau1*tau2*theta/2 + s**2*tau1*tau2 + s**2*tau1*theta/2 + s**2*tau2*theta/2 + s*tau1 + s*tau2 + s*theta/2 + 1))*K*Kc + (1/(s**4*tau1*tau2*tauI*theta/2 + s**3*tau1*tau2*tauI + s**3*tau1*tauI*theta/2 + s**3*tau2*tauI*theta/2 + s**2*tau1*tauI + s**2*tau2*tauI + s**2*tauI*theta/2 + s*tauI))*K*Kc + 1, s, 1/(s**3*tau1*tau2*tauI*theta + 2*s**2*tau1*tau2*tauI + s**2*tau1*tauI*theta + s**2*tau2*tauI*theta + 2*s*tau1*tauI + 2*s*tau2*tauI + s*tauI*theta + 2*tauI), 1/(s**3*tau1*tau2*tauI*theta/2 + s**2*tau1*tau2*tauI + s**2*tau1*tauI*theta/2 + s**2*tau2*tauI*theta/2 + s*tau1*tauI + s*tau2*tauI + s*tauI*theta/2 + tauI), 1/(s**3*tau1*tau2*theta/2 + s**2*tau1*tau2 + s**2*tau1*theta/2 + s**2*tau2*theta/2 + s*tau1 + s*tau2 + s*theta/2 + 1), 1/(s**4*tau1*tau2*tauI*theta/2 + s**3*tau1*tau2*tauI + s**3*tau1*tauI*theta/2 + s**3*tau2*tauI*theta/2 + s**2*tau1*tauI + s**2*tau2*tauI + s**2*tauI*theta/2 + s*tauI), K, Kc, tauD, taua, theta, domain='QQ'), {'gen': s, 'gens': (s, 1/(s**3*tau1*tau2*tauI*theta + 2*s**2*tau1*tau2*tauI + s**2*tau1*tauI*theta + s**2*tau2*tauI*theta + 2*s*tau1*tauI + 2*s*tau2*tauI + s*tauI*theta + 2*tauI), 1/(s**3*tau1*tau2*tauI*theta/2 + s**2*tau1*tau2*tauI + s**2*tau1*tauI*theta/2 + s**2*tau2*tauI*theta/2 + s*tau1*tauI + s*tau2*tauI + s*tauI*theta/2 + tauI), 1/(s**3*tau1*tau2*theta/2 + s**2*tau1*tau2 + s**2*tau1*theta/2 + s**2*tau2*theta/2 + s*tau1 + s*tau2 + s*theta/2 + 1), 1/(s**4*tau1*tau2*tauI*theta/2 + s**3*tau1*tau2*tauI + s**3*tau1*tauI*theta/2 + s**3*tau2*tauI*theta/2 + s**2*tau1*tauI + s**2*tau2*tauI + s**2*tauI*theta/2 + s*tauI), K, Kc, tauD, taua, theta), 'domain': QQ, 'polys': False})

Thank you

score 1 · Accepted Answer · answered May 26 '21 at 09:08

Your expression looks like the following:

The order of terms in the polynomial in the numerator goes from -1 to 3 while the denominator goes from 0 to 3. This means that it is bottom-heavy and you will get a 1/s term. This means that there does not exist constants a, b, c, d, e such that polin == (a)*s^4+ (b)*s^3+ (c)*s^2 + (d)*s + (e).

The best you can do is an approximation using a series expansion. It is important to note that one of the terms is 1/s and so you will get something of the form (a)*s^4+ (b)*s^3+ (c)*s^2 + (d)*s + (e) + (f)/s.

To solve this, you could manually get the coefficients or multiply polin by s so that all powers are nonnegative.

Method 1:

from sympy import *

Kc, K, tauI, taua, theta, tauD, tau1, tau2, s = symbols('K_c K tau_I tau_a theta tau_D tau_1 tau_2 s', real=True)

polin = 1 + K * (1 + taua * s) / ((1 + tau1 * s) * (1 + tau2 * s)) * (1 - s * theta / 2) / (1 + s * theta / 2) * Kc * (
            1 + 1 / (s * tauI) + tauD * s)

n=5  # highest power of polin is s**(n-1)
x = series(polin, s, n=n).removeO()

coefficients = [x.coeff(s**(i)) for i in range(-1, n)]
for i, coef in enumerate(coefficients):
    print(f"Coefficient of s**{i-1} is {coef}")

Method 2:

from sympy import *

Kc, K, tauI, taua, theta, tauD, tau1, tau2, s = symbols('K_c K tau_I tau_a theta tau_D tau_1 tau_2 s', real=True)

polin = 1 + K * (1 + taua * s) / ((1 + tau1 * s) * (1 + tau2 * s)) * (1 - s * theta / 2) / (1 + s * theta / 2) * Kc * (
            1 + 1 / (s * tauI) + tauD * s)

n=6  # highest power of polin is s**(n-2)
x = series(s*polin, s, n=n).removeO()

sp: Poly = Poly(x, s)  # s multiplied by polin
coefficients = sp.all_coeffs()
for i, coef in enumerate(coefficients):
    print(f"Coefficient of s**{n-i-2} is {coef}")

probabily not a good idea to use `from ... import *` (in general): https://stackoverflow.com/a/3615206/10465165 I think in this case it would be better to explicitly import functionalities needed from sympy — bk_, May 26 '21 at 09:24
@bk_ considering the nature of symbolic mathematics, SymPy's use is very interactive. You generally do not know which tools you will need (for example OP probably did not think to use `series`). Waiting an extra second to import everything is much less time consuming than having to rewrite `from sympy import ...` in the REPL. Half the issues for SymPy on stack overflow and github import * and it is "okay" according to Python's [documentation](https://docs.python.org/3/tutorial/modules.html#more-on-modules). — Chris du Plessis, May 27 '21 at 11:23
In that case I would recommend to either `import sympy` or use an alias: `import sympy as sy` (I do not know whether there is convention like numpy - np). And then use everything out of the package with `sy.function`. That way its clear out of which package `function` originated. If you import everything from a second package (e.g. `from numpy import *`, you're giving the reader of your script a hard time determining out of which package the functions you are using originate (are they inbuilt? are they from sympy? are they from numpy?) — bk_, May 27 '21 at 11:29

How to reorganize a mathematical expression into its polynomial form

1 Answers1