ASKSAGE: Sage Q&A Forum - Latest question feedhttps://ask.sagemath.org/questions/Q&A Forum for SageenCopyright Sage, 2010. Some rights reserved under creative commons license.Tue, 09 Jun 2020 04:15:48 -0500Monomial with power modulo nhttps://ask.sagemath.org/question/51855/monomial-with-power-modulo-n/ I need to implement the following formal structure
$$ax^\gamma, \gamma \in \mathbb{Z}/n\mathbb{Z}, a \in \mathbb{R}$$ and $x$ is a formal variable.
I tried
```
ZZ6 = Integers(6)
x = var('x')
x^ZZ6(9) + x^ZZ6(11)
```
And get value `x^8`, whereas I need `x^ZZ6(8)`.
How can I do this in sage?
only1saleTue, 09 Jun 2020 04:15:48 -0500https://ask.sagemath.org/question/51855/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)joaoffMon, 26 Mar 2018 04:45:25 -0500https://ask.sagemath.org/question/41760/A question on symbolic Matrices - unexpected Decimals in algebraic entryhttps://ask.sagemath.org/question/8942/a-question-on-symbolic-matrices-unexpected-decimals-in-algebraic-entry/Firstly, I'd just like to say that I am starting out with Sage adn Python at the moment, really impressed by the power of Sage and its potential! I am coming at this from an engineering and C++ background, and am currently trying to look at an engineering problem using some matrix linear algebra. I am setting up a from of stiffness matrix and am trying to extract the symbolic eigenvalues and eigenvectors from this. The stiffness matrix below is generated through a few different matrix multiplication operations, based on node and connectivity matrices. Despite trying to define these as being "SR" rings, I seem to end up with entries like (b + 1.0g + 1.0s) rather than the (b + g + s) I was expecting. I suspect this is due to me not understanding rings properly and how to apply them. I have tried the explanations in the Sage documentation, can anybody recommend some further reading or tutorials on this, because it does seem to be a vital topic in getting to grips with Sage?
I enclose my code below, in case anybody can spot a glaringly stupid thing that I've done...
many thanks,
Brian
N = matrix(QQ,[[1,0,-1,0],[0,1,0,-1]])
C = matrix(QQ,[[-1,0,1,0],[0,-1,0,1],[-1,1,0,0],[0,-1,1,0],[0,0,-1,1],[1,0,0,-1]])
D = C.transpose()
space = N.nrows()
M = N*D
M
II = identity_matrix(QQ,space)
b,s,g = var('b s g')
numMembers = M.ncols()
LVector = II
#loop through members to create LVector
for i in range(numMembers):
m1 = M.matrix_from_columns([i])
mt = m1.transpose()
mag = m1.norm()
Coeff = m1*mt/(mag^2)
size = m1.norm()
if i<2:
#for bars
L = -g*(II - Coeff) + b*Coeff
else:
L = g*(II - Coeff) + s*Coeff
L.factor()
if i == 0:
LVector = L
else:
LVector = LVector.block_sum(L)
print L.str()
print LVector.str()
Pre = D.tensor_product(II)
Post = C.tensor_product(II)
K = Pre*LVector
K = K*Post
Ks = K.change_ring(SR)
print Ks.str()
Eigen = K.eigenvectors_left()
print Eigen.str()BrianLoudonTue, 01 May 2012 00:22:55 -0500https://ask.sagemath.org/question/8942/