Revision history [back]

I am not aware on any "ready-made" solution. But deriving one is (relatively) easy :

Let's create the right field. Since you want numerical solutions, and did not specify any need for "high" precision, let's use the field of rational fractions of polynomials with complex coefficients, represented in CDF :

sage: P.<t>=PolynomialRing(CDF);P
Univariate Polynomial Ring in t over Complex Double Field
sage: FF=FractionField(P);FF
Fraction Field of Univariate Polynomial Ring in t over Complex Double Field

Create some random coefficients :

sage: alpha=[CDF(complex(2*(random()-0.5),2*(random()-0.5))) for p in range(6)];
....: alpha
[0.38301071001798026 - 0.8607425517858578*I,
 0.8424253929971652 + 0.01200809058669372*I,
 0.8708778554509589 - 0.36419026643134567*I,
 -0.27318942080730335 - 0.5784963036839359*I,
 -0.360559101520443 - 0.27656553577286846*I,
 0.6003797314834083 + 0.9787486618223218*I]

Create the rational fraction of interest :

sage: F=alpha[1]+alpha[2]*t+alpha[3]*t^2+(alpha[4]+alpha[5]*t^2)/(t^2-alpha[0]);
....: F
((-0.27318942080730335 - 0.5784963036839359*I)*t^4 + (0.8708778554509589 - 0.36419026643134567*I)*t^3 + (2.045375983144977 + 0.9771812732391981*I)*t^2 + (-0.020081486491519862 + 0.8890904001221609*I)*t - 0.6935529239631297 + 0.44394661938314584*I)/(t^2 - 0.38301071001798026 + 0.8607425517858578*I)

... and factorize it :

sage: fF=F.factor();fF
(-0.27318942080730335 - 0.5784963036839359*I) * (t - 0.8139787066102013 + 0.5287254720522134*I)^-1 * (t + 0.8139787066102013 - 0.5287254720522134*I)^-1 * (t - 1.6006015285417514 + 1.231139790778643*I) * (t - 0.49228617795296975 + 0.45260231434788073*I) * (t + 0.6737076471983565 - 0.20663751771064562*I) * (t + 1.3526456579041384 - 0.0031083744050989758*I)

The potential roots are the roots of each of these factors :

sage: PR=reduce(union, [q[0].roots(multiplicities=False) for q in fF if q[1]>0])
....: ; PR
[1.6006015285417514 - 1.231139790778643*I,
 -0.6737076471983565 + 0.20663751771064562*I,
 0.49228617795296975 - 0.45260231434788073*I,
 -1.3526456579041384 + 0.0031083744050989758*I]

But we must be careful of the potential poles (shit happens, at least numerically) :

sage: PP=reduce(union, [q[0].roots(multiplicities=False) for q in fF if q[1]<0])
....: ; PP
[-0.8139787066102013 + 0.5287254720522134*I,
 0.8139787066102013 - 0.5287254720522134*I]

By difference, get the effective roots :

sage: ER=set(PR)-set(PP);ER
{-1.3526456579041384 + 0.0031083744050989758*I,
 -0.6737076471983565 + 0.20663751771064562*I,
 0.49228617795296975 - 0.45260231434788073*I,
 1.6006015285417514 - 1.231139790778643*I}

Numerical check :

sage: [F.subs(t=q).abs() for q in ER]
[6.984039901725846e-15,
 4.258227406219383e-16,
 2.047869171424458e-15,
 2.4068798976194395e-16]

Close enough for government work...

Programming that in a "nice" function, with ability to pass options (such as the base ring and its precision, etc...) is left as an exercise for the reader... If well done, it might be added to Sage, if you think that the problem arises often enough to warrant such an inclusion.

By the way, Sage's symbolic solve is also able to get your roots directly :

sage: foo=SR(repr(F).replace("t","x")).solve(x)
sage: [q.rhs().n() for q in foo]
[-1.35264565790414 + 0.00310837440509715*I,
 -0.673707647198353 + 0.206637517710645*I,
 0.492286177952966 - 0.452602314347879*I,
 1.60060152854175 - 1.23113979077864*I]

But it seems to be much slower. I'll also leave you the delicate pleasure of discovering the form of the "exact" expression Sage uses to express those roots (hint to print them, plan for a couple of bedsheets...).

HTH,