Ask Your Question
1

What is the best way to work with ratios of polynomials?

asked 2018-03-26 04:45:25 -0500

joaoff gravatar image

updated 2018-03-27 03:26:10 -0500

What is the best way to work with ratios of polynomials with (real and complex) coefficients that can vary from very small to very big values? As I have to keep precision, I am working with symbolic expressions (type 'sage.symbolic.expression.Expression'), but as I operate with the expressions, they become bigger and bigger, and it is more and more difficult to plot them or find its roots, for example.

I perceived that making

# symbolic_expression is some symbolic polynomial with real or complex coefficients
symbolic_expression.polynomial(RR)
symbolic_expression.polynomial(CC)

simplifies the polynomials and eases computations but, as I understand, this converts the symbolic expressions to the real (RR) and complex (CC) rings, rounding the coefficients of the polynomials to machine precision and, in this manner, numerical errors would propagate in my computations. In this scenario, I would need to increase this conversion precision to an arbitrary value that guarantees correct results at the end of all computations.

I've read something about polynomial rings, symbolic ring, etc. I am studying this, but I still don't know what a ring is. Is symbolic ring the same thing as symbolic expressions? Does polynomial rings enable better performance and arbitrary precision?

Cordially

EDIT

I complement my question to answer to @B r u n o comment.

First, I would like to thank @B r u n o and @slelievre for the tips! I start to better understand things now.

Basically, I generate three polynomials. These polynomials are the basis of my work. I then operate on these polynomials generating other functions, e.g., as the ratio of these polynomials, and continue operating on the new functions, transforming them, plotting, finding the roots of the numerator and denominator polynomials, etc.

I'm slowly improving the code. I started with a plain script, then I created some functions to organize the code, in a procedural approach, and now I am moving to an object-oriented paradigm. This helps me better visualize and understand the problem. I even created a graphical interface to ease the interaction with the program, so as I can easily experiment with different input variables.

The code snippet below reflects the current state of the code it illustrates some inputs to test the code, how I generate my three basic polynomials and some derived functions. Sorry if the code seems silly :-/, I'm not used to develop software, I am still studying.

class FilterPrototype:

    def __init(self):
        pass

    def generate_polynomials(self):
        pass

    def generate_s21(self, e_poly, p_poly):
        return p_poly.polynomial(CC)/e_poly.polynomial(CC)

    def generate_s11(self, e_poly, f_poly):
        return f_poly.polynomial(CC)/e_poly.polynomial(CC)

    def generate_group_delay(self, s21):
        return -(1/s21(s=i*omega)*s21(s=i*omega).diff(omega)).imag()

    def plot(self, symb_express, bounds):
        pl = plot(symb_express, bounds[0], bounds[1])
        pl.show()
        return list(pl[0])

class ApproxArbitraryPhase(FilterPrototype):

    def __init__(self, spec_type, spec, e_order, p_order):
        self.spec_type = spec_type # type of specification
        self.spec = spec           # filter specification
        self.e_order = e_order     # order of the denominator polynomial
        self.p_order = p_order     # order of the numerator polynomial

    def generate_polynomials(self):
        listOmega = [] # List of frequencies (in radians per second)
        for a in range(self.e_order + 1):
            listOmega.append(a/self.e_order)

        if self.spec_type == "group delay":
            group_delay_spec(omega) = sage_eval(self.spec, locals={'omega':omega})
        elif self.spec_type == "phase":
            root.destroy()
            sys.exit("You cannot specify a phase. This is yet not implemented.")
        else:
            root.destroy()
            sys.exit("""The spec_type should be "group delay" or "phase".""")

        phase(omega) = group_delay_spec(omega).integral(omega)
        listPhases = [0]
        for a in range(1, self.e_order+1):
            listPhases.append(phase(listOmega[a]))

        alpha0 = listOmega[1]/tan(listPhases[1])
        listAlpha = [alpha0]
        for a in range(1, self.e_order):
            firstline = listOmega[a+1]^2 - listOmega[a]^2
            lastline = alpha0 - listOmega[a+1]/tan(listPhases[a+1])
            ic = 1
            while ic < a:
                lastline = listAlpha[ic] - ((listOmega[a+1]^2 - listOmega[ic]^2)/lastline)
                ic+=1
            listAlpha.append(firstline/lastline)

        ePolynomials = [1, s + listAlpha[0]]
        for a in range(2, self.e_order+1):
            ePolynomials.append(listAlpha[a-1]*ePolynomials[a-1] + (s^2 + listOmega[a-1]^2)*ePolynomials[a-2])

        ePoly = ePolynomials[self.e_order]
        pMod = sqrt(ePoly(s=i*omega)*ePoly(s=-i*omega))
        integrand = []
        for a in range(self.p_order+1):
            integrand.append(pMod*chebyshev_T(a, omega)/sqrt(1-omega^2))

        c = []
        for a in range(self.p_order+1):
            if a == 0:
                c.append((1/pi)*Rational(integrand[a].nintegral(omega, -1, 1)[0]))
            else:
                c.append((2/pi)*Rational(integrand[a].nintegral(omega, -1, 1)[0]))

        up = 0
        for a in range(self.p_order+1):
            up = up + c[a]*chebyshev_T(a, omega)

        eMod(omega) = (ePoly(s=i*omega)).polynomial(CC)
        den = lambda omega: up/abs(eMod(omega))
        pl = plot(den(omega), -1, 1) # to improve the precision of norm_factor, the plot_points=200 parameter of the plot function should be augmented
        plPoints = []
        for pair in list(pl[0]):
            plPoints.append(pair[1])
        normFactor = Rational(max(plPoints))

        pPoly(omega) = up/normFactor
        pPoly(s) = pPoly(omega=s/I)

        fMod2 = ePoly(s=s)*ePoly(s=-s) - pPoly(s=s)*pPoly(s=-s)

        solutions = fMod2.roots(ring=CC)
        fPoly = 0
        for solution in solutions:
            if solution[0].real() < 0 or (solution[0].real() == 0 and solution[0].imag() < 0):
                if fPoly == 0:
                    fPoly = (s - solution[0])
                else:
                    fPoly = fPoly*(s - solution[0])

        return (ePoly, fPoly, pPoly)


if __name__ == "__main__":
    # execute only if run as a script
    s = var("s")
    omega = var("omega", domain=RR)

    filtro = ApproxArbitraryPhase("group delay", "omega + 5/2", 4, 4)
#     filtro = ApproxArbitraryPhase("group delay", "omega + 9/2", 6, 6)
#     filtro = ApproxArbitraryPhase("group delay", "-omega + 27/5", 6, 6)
#     filtro = ApproxArbitraryPhase("group delay", "piecewise([([0,0.5], 2*omega + 143/20), ((0.5,1), 163/20), ([1,1], 163/20)])", 8, 8)
#     filtro = ApproxArbitraryPhase("group delay", "3/5*omega + 106/25", 6, 4)
    e,f,p = filtro.generate_polynomials()
    print(e)
    print(f)
    print(p(s))
    s21 = filtro.generate_s21(e, p(s))
    print(s21)
    s11 = filtro.generate_s11(e, f) # possui polinomio do numerador com coeficientes imaginarios
    print(s11)
    gd = filtro.generate_group_delay(s21)
    print(gd)
    points_s21 = filtro.plot(20*log(abs(s21(s=i*omega)),10), [-1.5, 1.5])
#     print(points_s21)
    points_s11 = filtro.plot(20*log(abs(s11(s=i*omega)),10), [-1.5, 1.5])
#     print(points_s11)
    points_gd = filtro.plot(gd, [-1.5, 1.5])
#     print(points_gd)
edit retag flag offensive close merge delete

Comments

The symbolic ring is the parent of symbolic expressions ("the place where they live"). And polynomial rings are specialized classes to handle polynomials specifically (rather than generic expressions).

How are your polynomials defined?

To work with more precision, you can either work with a real (or complex) field with more bits of precision (RealField(n) to get n bits of precision), or with interval or ball arithmetic to get guarantees (RealBallField(n) or RealIntervalField(n)), or with exactly represented algebraic numbers (using AlgebraicField() (complex) or AlgebraicRealField()). But the choice depends on your exact needs! (Finally note that do have quotients of polynomials over some ring R, you can use FractionField(PolynomialRing(R,'x')).)

B r u n o gravatar imageB r u n o ( 2018-03-26 06:24:44 -0500 )edit

1 answer

Sort by » oldest newest most voted
1

answered 2018-03-26 09:34:27 -0500

slelievre gravatar image

I highly recommend chapters 1, 2, 5, 7, 9, 11 of the book Calcul mathématique avec Sage, also available in an English version.

Note that to find roots of a quotient of polynomials, you can use the numerator only, and forget about the denominator.

There are well-known methods to count, find, isolate roots of polynomials. Some of them are implemented in Sage.

See these pages in SageMath's reference manual:

and you can also search the web for "root finding algorithms" and "root isolation".

edit flag offensive delete link more

Your Answer

Please start posting anonymously - your entry will be published after you log in or create a new account.

Add Answer

Question Tools

1 follower

Stats

Asked: 2018-03-26 04:45:25 -0500

Seen: 60 times

Last updated: Mar 27