Revision history [back]
 | 1 | initial version |
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!
