joaoff's profile - activity

I'm confused about the equality test of symbolic expressions. I saw the recently asked question https://ask.sagemath.org/question/417... and, coincidently, I read something similar on the first page of the second chapter of this book http://dl.lateralis.org/public/sagebo.... It says that Sage evaluates the following expressions as below,

sage: bool(arctan(1+abs(x)) == pi/2 - arctan(1/(1+abs(x))))
False
sage: bool(arctan(sqrt(2)) == pi/2 - arctan(1/sqrt(2)))
True

In Wikipedia (https://en.wikipedia.org/wiki/List_of...), we can verify the corresponding trigonometric identity.

$ arctan\left(x\right) + arctan\left(\frac{1}{x}\right) = \frac{\pi}{2}, \; if \; x>0$

and

$ arctan\left(x\right) + arctan\left(\frac{1}{x}\right) = -\frac{\pi}{2}, \; if \; x<0$

The above evaluation works as is in the version I am using.

sage: version()
'SageMath version 7.5.1, Release Date: 2017-01-15'

However, in Sage Math Cell or Cocalc ('SageMath version 8.1, Release Date: 2017-12-07') both tests evaluate to False.

Mathematica is not able to evaluate any of the expressions and return the entire equality test without evaluation.

Questions:

1) I understand that both tests should evaluate to True, am I right? If so, and as Sage 7.5.1 can verify the second equality, why it can't verify the first?

2) Why the newer version of Sage cannot verify the second equation anymore?

3) Differently from Mathematica, that returns True if both expressions are numerically equal, False otherwise and returns the test unevaluated if it cannot establish the equality, Sage returns True if it can prove that the difference between both expressions is zero and False otherwise. Is this correct?

Take a look at the "Tcl/Tk" subsection of the "Additional software" section of the installation guide. This should help you. http://doc.sagemath.org/html/en/insta...

Have you tried to use Tkinter? If it is not installed, you should be able to install it easily.

If you could try the fresnel_sin() command, tell me if it works. I’m still using Sagemath 7.5 :/

Hey! It seems that the support to Fresnel integrals starts with Sagemath 8.3 beta4.

https://trac.sagemath.org/ticket/24212

I didn’t read all the ticket, but if it is correct, it is very nice!

Ok, thanks for explaining!

However, this doesn't solve the problem. I cannot discretize the interval. This is what I do in the first workaround.

As you can see, the correct result is

(1.0053729181238868, -0.32744025779801017)

and you have obtained

(1.00331600582447, 0.23606797749978958)

I can obtain a more precise result with the first workaround mentioned above.

With your answer, I discovered that I need to upgrade my Sage system. I am using version 7.5.1 and in this version, your solution doesn't work, it returns

RuntimeError: ECL says: THROW: The catch MACSYMA-QUIT is undefined.

but it works in Sage Cell.

Hello!

I need to find the local maximum of an expression in an interval. It seems that the find_local_maximum() function has a problem when the argument expression contains the abs() function. This problem was reported to the sage-support group and two tickets were opened (https://trac.sagemath.org/ticket/24536 and https://trac.sagemath.org/ticket/24537), but it seems that it will take a while to be solved.

As a workaround, I can plot the function in that interval and look for the point with the maximum value, but this is not precise and produces a slightly different value each time the code is run.

Another workaround is to interface with Mathematica, but this is something I would like to avoid.

So, does anyone knows an alternative way to find the local maximum in an interval for an expression containing abs() functions? A piece of code is shown below as an example.

s = polygen(CC, "s")
omega = polygen(CC, "omega")
pi = CC.pi()
I = CC.0

p_order = 6
ePoly = s^6 + 2.52349407705767*s^5 + 4.57149144245730*s^4 + 4.95095921397024*s^3 + 3.70209076479638*s^2 + 1.64521458323101*s + 0.361999139300727

##### Generate P(s) polynomial #####
pMod(omega) = sqrt(ePoly(s=I*omega)*ePoly(s=-I*omega))
integrand = []
for a in range(p_order+1):
    integrand.append(pMod*chebyshev_T(a, omega)/sqrt(1-omega^2))

c = []
for a in range(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(p_order+1):
    up = up + c[a]*chebyshev_T(a, omega)

eMod = ePoly(s=I*omega)

### First workaround: Find the plot point with maximum value

den = lambda omega: up/abs(eMod)
pl = plot(den(omega), -1, 1)
plPoints = []
for pair in list(pl[0]):
    plPoints.append(pair[1])
# normFactor is multiplied by a constant to ensure that P(s) lies below the abscissa axis
normFactor = 1.00001*Rational(max(plPoints))

### Second workaround: Interfacing with Mathematica to calculate normFactor

# den = up/abs(eMod)
# mden = mathematica(den)
# mNormFactor = mden.NMaximize(-1 < omega < 1, omega)[1]
# normFactor = mNormFactor.sage()

pPoly = up/normFactor
pPoly = pPoly(omega=s/I)
pPoly = pPoly.polynomial(CC)
pPoly = pPoly.subs(x=s)

print(pPoly)

If I declare my expression as below, it works

den = lambda omega: (0.141162317194057*omega^6 + 0.173051768044343*omega^4 + 0.140157830888798*omega^2 + 0.359575336011182)/abs(-1.00000000000000*omega^6 + 2.52349407705767*I*omega^5 + 4.57149144245730*omega^4 - 4.95095921397024*I*omega^3 - 3.70209076479638*omega^2 + 1.64521458323101*I*omega + 0.361999139300727)
print(find_local_maximum(den,-1,1))

the result is

(1.0053729181238868, -0.32744025779801017)

If I declare my expression as

den = lambda omega: up/abs(eMod)
print(find_local_maximum(den,-1,1))

it returns

((0.141162317194057*omega^6 + 0.173051768044343*omega^4 + 0.140157830888798*omega^2 + 0.359575336011182)/abs(-1.00000000000000*omega^6 + 2.52349407705767*I*omega^5 + 4.57149144245730*omega^4 - 4.95095921397024*I*omega^3 - 3.70209076479638*omega^2 + 1.64521458323101*I*omega + 0.361999139300727), 0.99999996297566163)

Is that all you want?

One way to do that.

for k in range(9):
    if k != 3:
        if k != 4:
            print(k)
        else:
            pass
    else:
        print(20)
        # i want to skip the next iteration

Does this help? this is just an idea. You can improve it, of course.

var('x')
anti=integrate((4*x - sqrt(-x^2 + 1))/(sqrt(-x^2 + 1) + 5), x)
#anti=integrate(1/(sqrt(x + 1)*sqrt(-x + 1) + 5), x)
if str(anti).find("integrate") != -1:
    print("Result contains unevaluated integral!")
anti

Output:

Result contains unevaluated integral!
-x - 4*sqrt(-x^2 + 1) + 5*integrate(1/(sqrt(x + 1)*sqrt(-x + 1) + 5), x) + 20*log(sqrt(-x^2 + 1) + 5)

If you use the timeit command to profile your code, adjusting the command options (e.g., preparse, number, repeat, precision) adequately, you should obtain more or less the same timing every time you measure your code execution time, either in sagemath cloud (cocalc) or running it through the terminal.

Is there really no way to apply an assumption on the indeterminate when working on a polynomial ring? this is intriguing.

Your example with timeit and @slelievre's example with timeit.

164 ms against 13.2 ms

Well, there is nothing wrong with your solution but, as @slelievre suggested, I am trying to avoid the symbolic ring. In fact, I started working in the symbolic ring, but I began to face several difficulties (probably due to my inexperience). Due to this, I changed the code to use the polynomial ring over the complex field and, despite losing precision, found easier to work this way and I could improve code performance. Maybe later I will try again to work in the symbolic ring, or work in polynomial rings over rationals, to regain an infinite precision.

By the way, thanks for your response!

Great answer, thanks a lot!

I'm changing my code to use polynomial rings instead of the symbolic ring. The problem is that sometimes the response changes to the symbolic ring without me realizing or knowing why. Your answers helped me better understand what is happening in my code. Concerning the formatting, clarity and modularity of the code, this is something that I have great difficulty. Generally, I tend to use big descriptive variable names to assure that I will immediately understand the code when returning to it after not working on it for a while.

When working with polynomial rings with complex coefficients, how to pass the assumption that the unknown variable belongs to the Reals set?

I am trying to recover and plot the imaginary part of the the ratio of polynomials arg2, but as sage assumes the independent variable is complex, it takes a long time and returns a huge answer. But I think the answer could be much simpler if it assumes that the independent variable belongs to the Reals set. The assume() function doesn't seems to work in this case.

arg has complex terms in the numerator and denominator. arg2 has complex terms only on the numerator. So, it should be easy to obtain the imaginary part of arg2.

s = polygen(CC, "s")
omega = polygen(CC, "omega")
I = CC(I)

# assume(omega, "real") # this does not work

s21 = 0.894292331092066/(s^4 + 2.14081977623463*s^3 + 3.15237117897600*s^2 + 2.31898630138664*s + 0.902488008823108)

s21_omega = s21.subs(s=I*omega)
s21_omega_diff = s21_omega.derivative()
arg = 1/s21_omega*s21_omega_diff

num = arg.numerator()
den = arg.denominator()
den_conj = den.map_coefficients(lambda z: z.conjugate())

num2 = (num*den_conj).map_coefficients(lambda z: 0 if abs(z.real())<1e-10 and abs(z.imag())<1e-10 else z.imag()*I if abs(z.real())<1e-10 and abs(z.imag())>=1e-10 else z.real() if abs(z.real())>=1e-10 and abs(z.imag())<1e-10 else z)
den2 = (den*den_conj).map_coefficients(lambda z: 0 if abs(z.real())<1e-10 and abs(z.imag())<1e-10 else z.imag()*I if abs(z.real())<1e-10 and abs(z.imag())>=1e-10 else z.real() if abs(z.real())>=1e-10 and abs(z.imag())<1e-10 else z)
arg2 = num2/den2
print(arg2)
print(imag(arg2))
# plot(-imag(arg2), (-1.5, 1.5)).show()

Output:

(-3.19903509380033*omega^15 - 1.71213939841908*I*omega^14 + 9.63823791915898*omega^13 + 3.11426578979509*I*omega^12 - 15.8129908994689*omega^11 - 4.60245136409720*I*omega^10 + 13.7326241083798*omega^9 + 1.24770230131999*I*omega^8 - 9.58497293293397*omega^7 - 0.760732566992820*I*omega^6 + 4.72291963120453*omega^5 - 2.52135706735033*I*omega^4 - 2.44038907746050*omega^3 - 0.463630243025939*I*omega^2 + 0.203401406783201*omega - 1.36326886734867*I)/(0.799758773450081*omega^16 - 2.75378226261685*omega^14 + 5.27099696648963*omega^12 - 5.49304964335183*omega^10 + 4.79248646646703*omega^8 - 3.14861308746971*omega^6 + 2.44038907746050*omega^4 - 0.406802813566401*omega^2 + 0.530548112702672)

2.55845638284151*imag_part(omega)^31/(0.639614095710378*imag_part(omega)^32 + 10.2338255313661*imag_part(omega)^30*real_part(omega)^2 + 76.7536914852462*imag_part(omega)^28*real_part(omega)^4 + 358.183893597743*imag_part(omega)^26*real_part(omega)^6 + 1164.09765419317*imag_part(omega)^24*real_part(omega)^8 + 2793.83437005430*imag_part(omega)^22*real_part(omega)^10 + 5122.02967846394*imag_part(omega)^20*real_part(omega)^12 + 7317.18525490165*imag_part(omega)^18*real_part(omega)^14 + 8231.83341181278*imag_part(omega)^16*real_part(omega)^16 + 7317.18525493145*imag_part(omega)^14*real_part(omega)^18 + 5122.02967846394*imag_part(omega)^12*real_part(omega)^20 + 2793.83437005430*imag_part(omega)^10*real_part(omega)^22 + 1164.09765419317*imag_part(omega)^8*real_part(omega)^24 + 358.183893597743*imag_part(omega)^6*real_part(omega)^26 + ... (huge answer)

Thank you!

Thank you, I was never going to find this out.

Does anyone knows why after the substitution below, the scaled_s21 variable becomes of type 'sage.symbolic.expression.Expression' instead of continuing as 'sage.rings.fraction_field_element.FractionFieldElement_1poly_field'?

I'm confused. I just wanted to use the numerator() and denominator() methods, but with the type changing, I'm unable.

Thnak you!

s = polygen(CC, "s")
p = 0.894289785676221
e = s^4 + 2.14081977623463*s^3 + 3.15237117897600*s^2 + 2.31898630138664*s + 0.902488008823108
s21 = p/e
print(s21)
print(type(s21))
print(s21.parent())
omega0 = sqrt(2*pi*750e6*2*pi*1250e6)
scaled_s21 = s21.subs(s=omega0/(2*pi*1250e6-2*pi*750e6)*(s/omega0 + omega0/s))
print(scaled_s21)
print(type(scaled_s21))
print(scaled_s21.parent())

Output:

0.894289785676221/(s^4 + 2.14081977623463*s^3 + 3.15237117897600*s^2 + 2.31898630138664*s + 0.902488008823108)
<class 'sage.rings.fraction_field_element.FractionFieldElement_1poly_field'>
Fraction Field of Univariate Polynomial Ring in s over Complex Field with 53 bits of precision
0.894289785676221/(1.00000000000000*(3.75000000000000e9*pi/s + (1.00000000000000e-9)*s/pi)^4 + 2.14081977623463*(3.75000000000000e9*pi/s + (1.00000000000000e-9)*s/pi)^3 + 3.15237117897600*(3.75000000000000e9*pi/s + (1.00000000000000e-9)*s/pi)^2 + 8.69619863019990e9*pi/s + (2.31898630138664e-9)*s/pi + 0.902488008823108)
<type 'sage.symbolic.expression.Expression'>
Symbolic Ring

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)

Does anyone know a function in Sage equivalent to Mathematica _Chop()_ function ( http://reference.wolfram.com/language... )? I couldn't find any.

In fact, I would like to drop the imaginary parts that are close to zero of complex polynomial coefficients.

This shouldn't be difficult to implement but is something that I use a lot, it is strange not find an implementation.

Cordially

Very interesting! I can't see the applicability of the mentioned function in this particular problem, but it probably can help with my other problem, plotting big symbolic expressions (https://ask.sagemath.org/question/416...).

In Mathematica, when plotting a very big and complex symbolic expression, I generally use the Evaluate[] function (e.g. http://reference.wolfram.com/language...). I was looking for a similar function in Sage and, apparently, fast_callable is this function. I will give it a try.

Thanks a lot!

Does anyone know how to plot a ratio of polynomials where the coefficients of the polynomials vary from extremely big to extremely small values? e.g.

var("omega")

numerator = 1/434425210923785365644729072814805003095647594571477081143688112179*omega^82 - 1/282168800363107583713154019070571550253203986968116550714*omega^80 + 1/376705284799486258761328410091603792725536189789*omega^78 - 1/775984649289050778795273868963329671283*omega^76 + 1/2194552498783346688198647992959*omega^74 - 1/7997055354595006175132*omega^72 + 37/1335393102730856*omega^70 - 455327/89425453412*omega^68 + 6451507166/8156449*omega^66 - 41921205646451/398*omega^64 + 12148488559833861120*omega^62 - 1222797178046083079207911425*omega^60 + 107973546113555797789307958819028992*omega^58 - 8389589284510832786372039497137019252899840*omega^56 + 574047573454150250336814935212839009181133950681089*omega^54 - 34513862440990797385476079081710948473169710438951140982784*omega^52 + 1811420748484519385045717526620321156600093813148380753880114266113*omega^50 - 81787566131583729215998508705538459221433180523307219172710272077796474880*omega^48 + 3074808395164144704086261320345404352244943366045903809072779420014923038874140672*omega^46 - 88158128188124577566928559996352418914343887707417923271162886096982268067866251309875201*omega^44 + 1276001990602982330745542293770646665590092983251428019005953641290986839752143036376864694206464*omega^42 + 50369501995070658072835654251040011409949189464195973029581067274496453807222446814723919181288107409408*omega^40 - 5422646273721796408605794927447900588248745194267196182599017580889204367989136820721458419743182469561086640129*omega^38 + 294713197016069575830983289812182224150427565530453919249784823279683263369786981460098799522297673811884617147560755200*omega^36 - 12215223571528656750794949703498245908596636883090964706566516953450365851729599110977560326879247539894799131557796203328438273*omega^34 + 421566491928721680419007829094184826541933844529415542073312144440199053433647609288595755917957129357084944852204800077564039202340864*omega^32 - 12516212267147404429229647731387721025123388989961149553475986047447463987490660724352358293390873442017645025300898526431390960031953078714369*omega^30 + 324384318879686142437160140587334223782579714788774469706503221987046551849501681825766358071172441792448781390801119733389149415560400234273231077377*omega^28 - 7387292185075462418862140718878310563299525368703143556507183065676300148539804750553952220382528962008497259504339924871164528693646839499632967763533758464*omega^26 + 148147576549670480176424742570505113382496592810105724265276479687408049874080107426078378616519348224098107322066761943275443363035890238181101932475433564088827904*omega^24 - 2614352983825249603402079949177727066634757002976206933577141540374477238146660413804057344751020839726456625121445634542631825250862617604413287932677361615994705332404224*omega^22 + 40472930950614252427322391199345022952821847319764623404423487895464586861905495460361611937362853552737369167685282072188371812469915846107436301498874685121753362553033067593728*omega^20 - 546796479239276734969817343152988884071617215095039987059486720343626501580823236262129279129106406271648455070886215057036604477324486364130961502454958245222249278481585821602157690881*omega^18 + 6398345482963010405651688543654226237482870669667751595384527999380305228718755002303831885277585496704965845755258790817010267468004294456270161640983423358987394850223629580900096084482719744*omega^16 - 64180410520272272525286908376473230431734138441676083665323637745196749056460187115412218416398327877779227771800677121061088325929088938579331150090894803951427168305021791509246069685553904083795968*omega^14 + 544208653291391185309536180524159975422536018745994149714204580484188377200450473903397939463023201727188280315535659623140071365231759008260639937606202290659210954860546302434710221600085853152285286203393*omega^12 - 3827148869236131645057753748864256011706470094795063609207751917613788002030507913841732628668035491051373157906167394821657884872636999755793875871427355836181313542579221345353770527185046940717229548007726252032*omega^10 + 21731619560104735623412248166557020474247938972624063082336680903151973764168130050518970368891893871902596734207676031052134545024607407417007348957518735674156914150659282214614029610623838986362932572240851152089907201*omega^8 - 95767367248890998859080840076211432090249689089635000943797128834097631391862494890486062534305250474327151750740126931966007985728703285377423534457675359002872120478001359450002524654691233467116529177706330479483962206978048*omega^6 + 307398380835349522251122973685605867590098291648854481500785164487397373217556109794667977004956100249926660979665856367573287860824766046271883357679682890511797573992243503579120101379656959151097962039219471247899783414197876948992*omega^4 - 639479912680222769699574758361470606739916016079184564965894672578774414095139913044892342897844958391466586738992899235402376325166165958853889423782954279875630911656947991004581824214687332076761443029188086281653700846286211570814943233*omega^2 + 647222323629001353412618991139040623627834152189367973808816854381809885796462656684768497116303794227646539449954073178816871620703525349577186183599165579353490438287333256004323407793268413622639886629815730954228627165293634675091625356558336

denominator = 1/137761281603135848379924667497137087691869145590694743389457616256325*omega^84 - 1/83086834766162983807691746766150178643427527433316275522827*omega^82 + 1/102667767044170877623975342730336352676592625615621*omega^80 - 1/195052973429231983176184037372982355474872*omega^78 + 1/506764970405918852128607114944932*omega^76 - 1/1689091026941888449775972*omega^74 + 1/6938475729210507*omega^72 - 5933/202766817954*omega^70 + 210652963/41684674*omega^68 - 3785514892410/5023*omega^66 + 98176715539323913*omega^64 - 11274496021929013098315776*omega^62 + 1149760802785547342617432310153216*omega^60 - 104740316526399710014129038296835497656321*omega^58 + 8564675098852218158330202009453158486665628483585*omega^56 - 631107412233357996092481381992488681322377443105002487809*omega^54 + 42040826550209669791480148216500925719160340847286518196157284352*omega^52 - 2538156657414023370036118838246125316621096397799706645672301006268923904*omega^50 + 139157701207028251044471300427394350671149356489494088577874141239880561868144640*omega^48 - 6938806439383168134484769650927631938020954059414471869084090238768620526122595210231809*omega^46 + 314993119151231424855218501630778642353808903970689645637119382608673180637105584990275091038209*omega^44 - 13026311952157375634499439761001799003899804334363328258128277741022032909434938363115295924632479072256*omega^42 + 490832913386158964971261411043934330175709099657034165101518631974970770766908644352744356524979896994217590784*omega^40 - 16848063598754041432195740457548489087523338761447189453678864215130767690468257118676488451717823178276208643556769792*omega^38 + 526508123680113455403918388001376421311837933201943096641237980895995131957184498765311703601193482744568171270729639711997953*omega^36 - 14964040137529112789336002405814742574969978116283981052553508948071227427385569689635300763344329290377090567018871249668589630783489*omega^34 + 386219588613305226787305963970336593407234250843368709064728339440768265801980913517208370322485059415982349821971022988796755219950914764801*omega^32 - 9034392239290768500236248643616567747019665018495660272110199314250181111897744491977335480100194146015876816462940326122858266556363974484342865921*omega^30 + 191046534415310850754936618207854871774036324685618690113988679137516259702053348085287489422888909156122197081115795038125721730271097600962795307182063616*omega^28 - 3640615079002674868268803209292339284671635529598000125401708718259709095257942353116876576750774280418512074696423958832428090501175351634961736513120496483565568*omega^26 + 62273172374254986415765227778440797173097273754339193891373421481657115467189617254708606757292829726836650836046370485734292092085898801412029454712950542649407572017153*omega^24 - 951531652629767939573304547715201546281156310628865493286337971797183608774930095570943761071381252056585979625290081231116086193892640991438658207371842803101206913382977896449*omega^22 + 12911217909000424455313473021628455155346806719690986449907489732516140526838114148296344376049550653353703159891539620541282581173653621727840150941148476570112554868892393249471201280*omega^20 - 154437621072641585160192364210764638535150742299453909620719685765830187029475768884884341432582322941793749721471314800162159259791070934517719047186939224330486171808759549778279712277135361*omega^18 + 1613671375800221558063007173458743660146871072783243100013792246050123478339713407322883872334136082537546559536378103103080560054853188334529829307580805332330325585126522606938881345808027181121537*omega^16 - 14559082441673569008977537112513409409617326807467488077449500018188832706180435122319535717730339986042614839268533169311805347564162239812492671767653840792592722234547429968326169980798063761096993406976*omega^14 + 111744222840860258586904336704874254832742992254602728119554334607230626003726011474866600712329730965417565026295984240724114851734359742329870679353469856198907573647755676148914916780030398787721120780818317313*omega^12 - 715267855855396257045094528442779490633060111715529941383146691210931540467322126692565592976853593736823767984638431457105145557082338533900534849060590000205707733135840156253287305930305666272715830450207103255576577*omega^10 + 3714905633019119124314145982488666373006597207416225659385733815586350823207206537470310618757702145455955559354122274345149874879416613750368567586179362486479162207672443597112579110941483842098201059445524219729129391521792*omega^8 - 15039534985098642489127495070878079816478699521964932372348795943530923555616039497846344032659763466508345351341047837858908515818598488134274908478846308392419344474084655780759561686415523039911831647254805083163507636676075716609*omega^6 + 44523108397807453960903723278683589626375425344092590311786839673418629587601620660013463353874988948445949661864697766839069544771099155966119222987528998667273718783412135043865598080899844899439778708053188505321146926956138446907768832*omega^4 - 85727466500775256683870209515231431180659141489998596286522465643370098508698223577125979654943088313884931697767434884070649253899301725049084279835311596901096918650207437145768941287635219222105620194647347139209573654378699720526506978443265*omega^2 + 80567082948458938117841918763525355285225446281873107502306332907427746631501737918800317997386306466900211456674973441765747419882727073257243150087944844059450542039822520151117935401484171322088826759710569301983550807117385762337391493721136037889

#numerator = 1/434425210923785365644729072814805003095647594571477081143688112179*omega^82 + 647222323629001353412618991139040623627834152189367973808816854381809885796462656684768497116303794227646539449954073178816871620703525349577186183599165579353490438287333256004323407793268413622639886629815730954228627165293634675091625356558336

#denominator = 1/137761281603135848379924667497137087691869145590694743389457616256325*omega^84 + 80567082948458938117841918763525355285225446281873107502306332907427746631501737918800317997386306466900211456674973441765747419882727073257243150087944844059450542039822520151117935401484171322088826759710569301983550807117385762337391493721136037889

gd = numerator/denominator
plt3 = plot(gd, (2*pi*0.98*10^3, 2*pi*1.02*10^3))
plt3.show(gridlines=True)

As you can perceive, I tried to rationalize the coefficients, even though, I get the following message, with an empty plot

verbose 0 (3757: plot.py, generate_plot_points) WARNING: When plotting, failed to evaluate function at 200 points.
verbose 0 (3757: plot.py, generate_plot_points) Last error message: '(34, 'Numerical result out of range')'

The plot should resemble an "M" letter centered at 2pi10^3.