Imaginary part of a ratio of polynomials

asked 2018-06-15 16:23:09 +0100

joaoff
220 ●1 ●6 ●21

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!

edit retag flag offensive close merge delete

add a comment

2 Answers

Sort by » oldest newest most voted

answered 2018-06-15 23:07:12 +0100

slelievre
17704 ●22 ●160 ●350 http://carva.org/samue...

The code in the question can be simplified in several ways.

Let us first provide simple code to plot the desired function, and then take the opportunity of this question to review the code in the question and offer some advice.

Plot imaginary part of rational function with complex coefficients

Here omega is renamed w and a few more names are changed for brevity and readability.

sage: R.<s> = CC[]
sage: w = polygen(CC, 'w')
sage: I = CC.0
sage: d = [0.902488008823108, 2.31898630138664, 3.15237117897600, 2.14081977623463, 1.0]
sage: s21 = 0.894292331092066/R(d)
sage: w21 = s21(I*w)
sage: w21_diff = w21.derivative()
sage: arg = w21_diff/w21
sage: f = lambda z: -arg(z).imag()
sage: plot(f, (-1.5, 1.5))

Voilà!

Avoid Sage's symbolic ring

The essential point is to avoid Sage's symbolic ring if at all possible.

What made the original code so slow was that imag(arg2) returned an expression living in Sage's symbolic ring, which is best avoided.

Code formatting, clarity, modularity

Since arg2 is just a more complicated expression of arg, it is not needed.

But let me still say how I would have reformatted its definition for readability.

After defining num, den, den_conj, one could define num2 and den2 as:

sage: 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)

sage: 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)

from there it appears that the goal is to suppress any real part or imaginary part that less than 1e-10, which can be made more readable as follows:

sage: epsilon = 1e-10
sage: truncate_x = lambda x: 0 if abs(x) < epsilon else x
sage: truncate_z = lambda z: (truncate_x(z[0]), truncate_x(z[1]))

sage: num2 = (num*den_conj).map_coefficients(truncate_z)
sage: den2 = (den*den_conj).map_coefficients(truncate_z)

This alternative formulation makes it a lot easier to change the truncation threshold; in case 1e-10 has to be changed to 1e-8 or 1e-12, it's better if it can be changed in only one place!

edit flag offensive delete link

Comments

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.

joaoff ( 2018-06-16 23:46:40 +0100 )edit

add a comment

answered 2018-06-19 08:17:48 +0100

Emmanuel Charpentier
7764 ●6 ●50 ●153

Your problem is that

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

thus prohibiting Sage to simplify expressions involving real_part, imag_part and others...

And this is attributable to

s = polygen(CC, "s")

omega = polygen(CC, "omega")

s and o no longer are symbolic expressions ; therefore, Sage can't assume anything about them.

Tell me what's wrong with the trivial solution :

var("s")
s21 = 0.894292331092066/(s^4 + 2.14081977623463*s^3 + 3.15237117897600*s^2 + 2.31898630138664*s + 0.902488008823108)
var("omega", domain="real")
o21=s21.subs(s=I*omega)
(o21.diff()/o21).imag_part().plot(omega,-3,3)

whose result you can meditate upon...

edit flag offensive delete link

Comments

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!

joaoff ( 2018-06-19 09:46:36 +0100 )edit

I am trying to avoid the symbolic ring

That is a useful heuristic. but it is only an heuristic. In the present case, there is, as far as I know, no way to express the fact that an indeterminate on the polynomial ring over CC is real. Therefore, you lose all the simplifications that can be done using that fact.

And your heuristic becomes a superstition.

In my (modified) example (I just added a timing command), the plot uses all of ... 164 ms. What is the time with Serge's solution ?

Emmanuel Charpentier ( 2018-06-19 17:41:25 +0100 )edit

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

164 ms against 13.2 ms

joaoff ( 2018-06-19 18:59:27 +0100 )edit

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

joaoff ( 2018-06-19 19:10:26 +0100 )edit

add a comment

Your Answer

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

Add Answer

Imaginary part of a ratio of polynomials

2 Answers

Plot imaginary part of rational function with complex coefficients

Avoid Sage's symbolic ring

Code formatting, clarity, modularity

Comments

Comments

Your Answer

Question Tools

Stats

Related questions

Imaginary part of a ratio of polynomials edit

2 Answers

Plot imaginary part of rational function with complex coefficients

Avoid Sage's symbolic ring

Code formatting, clarity, modularity

Comments

Comments

Your Answer

Question Tools

Stats

Related questions

Imaginary part of a ratio of polynomials