ASKSAGE: Sage Q&A Forum - RSS feedhttps://ask.sagemath.org/questions/Q&A Forum for SageenCopyright Sage, 2010. Some rights reserved under creative commons license.Mon, 26 Mar 2018 16:34:27 +0200- What is the best way to work with ratios of polynomials?https://ask.sagemath.org/question/41760/what-is-the-best-way-to-work-with-ratios-of-polynomials/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)Mon, 26 Mar 2018 11:45:25 +0200https://ask.sagemath.org/question/41760/what-is-the-best-way-to-work-with-ratios-of-polynomials/
- Comment by B r u n o for <p>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.</p>
<p>I perceived that making</p>
<pre><code># symbolic_expression is some symbolic polynomial with real or complex coefficients
symbolic_expression.polynomial(RR)
symbolic_expression.polynomial(CC)
</code></pre>
<p>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.</p>
<p>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?</p>
<p>Cordially</p>
<p><strong><em></em></strong><em> EDIT <strong></strong></em></p>
<p>I complement my question to answer to <a href="/users/8630/b-r-u-n-o/">@B r u n o</a> comment.</p>
<p>First, I would like to thank <a href="/users/8630/b-r-u-n-o/">@B r u n o</a> and <a href="/users/1092/slelievre/">@slelievre</a> for the tips! I start to better understand things now.</p>
<p>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.</p>
<p>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.</p>
<p>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.</p>
<pre><code>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)
</code></pre>
https://ask.sagemath.org/question/41760/what-is-the-best-way-to-work-with-ratios-of-polynomials/?comment=41762#post-id-41762The 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'))`.)Mon, 26 Mar 2018 13:24:44 +0200https://ask.sagemath.org/question/41760/what-is-the-best-way-to-work-with-ratios-of-polynomials/?comment=41762#post-id-41762
- Answer by slelievre for <p>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.</p>
<p>I perceived that making</p>
<pre><code># symbolic_expression is some symbolic polynomial with real or complex coefficients
symbolic_expression.polynomial(RR)
symbolic_expression.polynomial(CC)
</code></pre>
<p>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.</p>
<p>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?</p>
<p>Cordially</p>
<p><strong><em></em></strong><em> EDIT <strong></strong></em></p>
<p>I complement my question to answer to <a href="/users/8630/b-r-u-n-o/">@B r u n o</a> comment.</p>
<p>First, I would like to thank <a href="/users/8630/b-r-u-n-o/">@B r u n o</a> and <a href="/users/1092/slelievre/">@slelievre</a> for the tips! I start to better understand things now.</p>
<p>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.</p>
<p>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.</p>
<p>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.</p>
<pre><code>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)
</code></pre>
https://ask.sagemath.org/question/41760/what-is-the-best-way-to-work-with-ratios-of-polynomials/?answer=41774#post-id-41774I highly recommend chapters 1, 2, 5, 7, 9, 11 of the book
[Calcul mathÃ©matique avec Sage](http://sagebook.gforge.inria.fr/),
also available in an
[English version](https://members.loria.fr/PZimmermann/sagebook/english.html).
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:
- [SageMath reference manual: Isolate Real Roots of Polynomials](https://doc.sagemath.org/html/en/reference/polynomial_rings/sage/rings/polynomial/real_roots.html)
- [SageMath reference manual: Isolate Complex Roots of Polynomials](https://doc.sagemath.org/html/en/reference/polynomial_rings/sage/rings/polynomial/complex_roots.html)
and you can also search the web for "root finding algorithms" and "root isolation".Mon, 26 Mar 2018 16:34:27 +0200https://ask.sagemath.org/question/41760/what-is-the-best-way-to-work-with-ratios-of-polynomials/?answer=41774#post-id-41774