CODEWITHSUNDEEP
pythonnumpyscipywolfram-mathematicasolver
Alexander McFarlane
Alexander McFarlane

Reputation: 11293

Numerically Representing Mathematica's Root Object in Open-Source Language

Question

I would like to reproduce a Root[] object, ideally in a python function.

Is there any particular library that would be suited for this process?

Attempts

If I understand the Root[] function properly, it is simply finding the nth degree root of a polynomial and so I and taking a stab that numpy.roots would suffice by taking #1 as the argument x in the numpy docs.

Background

I have a number of 5th order polynomial root which cannot be reduced with ToRadicals due to the order that were obtained from a particularly nasty Inverse Laplace Transform.

Minimal Mathematica example

r = 1/τ;
ct = Cos[θ];
r2 = r^2;
r3 = r^3;
r4 = r2^2;
ct2 = ct^2;
ct3 = ct^3;
ct4 = ct2^2;
ct5 = ct^5;
ct6 = ct2^3;
p2 = ϕ^2;

fn := Root[2*ρ*p2*r^5 + 2*ρ*p2*r^5*ct - 2*ρ^2*p2*r^5*ct - 2*ρ*p2*r^5*ct2 - 2*ρ*p2*r^5*ct3 + 2*ρ^2*p2*r^5*ct3 + (r4 + 4*p2*r4 + 4*ρ*p2*r4 + r4*ct - 2*ρ*r4*ct + 4*p2*r4*ct - 2*ρ*p2*r4*ct - 2*ρ^2*p2*r4*ct - r4*ct2 - 4*p2*r4*ct2 - r4*ct3 + 2*ρ*r4*ct3 - 4*p2*r4*ct3 + 6*ρ*p2*r4*ct3 -  2*ρ^2*p2*r4*ct3)*#1 + (4*r3 + 8*p2*r3 + 2*ρ*p2*r3 +  3*r3*ct - 6*ρ*r3*ct + 4*p2*r3*ct - 4*ρ*p2*r3*ct - 3*r3*ct2 - 4*p2*r3*ct2 + 2*ρ*p2*r3*ct2 - 2*r3*ct3 + 4*ρ*r3*ct3)*#1^2 + (6*r2 + 4*p2*r2 + 3*r2*ct -  6*ρ*r2*ct - 3*r2*ct2 - r2*ct3 +  2*ρ*r2*ct3)*#1^3 + (4*r + r*ct - 2*ρ*r*ct -  r*ct2)*#1^4 + #1^5 &, 5]

Upvotes: 0

Views: 167

Answers (1)

Warren Weckesser
Warren Weckesser

Reputation: 114811

If you are interested in symbolic calculations, you can use SymPy. In particular, SymPy has polynomial objects and the classes RootOf and CRootOf to represent the roots of polynomials.

For example,

In [103]: from sympy import Symbol, poly

In [104]: x = Symbol('x')

In [105]: p = poly(x**4 - 3*x**2 + x - 1)

In [106]: p
Out[106]: Poly(x**4 - 3*x**2 + x - 1, x, domain='ZZ')

In [107]: p.root(0)
Out[107]: CRootOf(x**4 - 3*x**2 + x - 1, 0)

CRootOf(poly, k) is a placeholder for the kth root of the polynomial. To find its numerical value, use the .evalf() method:

In [109]: p.root(0).evalf()
Out[109]: -1.94397243715073

Here are the numerical values of all the roots:

In [110]: [p.root(k).evalf() for k in range(p.degree())]
Out[110]: 
[-1.94397243715073,
 1.66143946800762,
 0.141266484571554 - 0.538201812325831*I,
 0.141266484571554 + 0.538201812325831*I]

Upvotes: 2

Related Questions