Getting all (complex) solutions of a non polynomial equation

asked 2017-12-26 15:08:51 +0200

yoda789
13 ●1 ●1 ●5

updated 2017-12-26 21:12:50 +0200

Hi !

I was used to solve the following equation with Mathematica. \begin{equation} \alpha_1 + \alpha_2x + \alpha_3x^2 + x^4 + \frac{\alpha_4}{x^2-\alpha_0} + \frac{\alpha_5 x^2}{x^2-\alpha_0}=0 \end{equation} where $\alpha_i$ are constants.

The Mathematica function "Solve" gives me all the numerical roots of this non polynomial equation very easily. These roots can be real or complex.

I'm a very beginner at Sage. I have tried several methods to solve this equation automatically but it seems that all methods I've used work only for polynomial equations. Here they are :

alpha0 = 0.25
alpha1 = -2.5
alpha2 = 6.9282
alpha3 = -5.5
alpha4 = 0.5
alpha5 = -0.5
x = var('x')
eq = alpha1 + alpha2*x + alpha3*x**2 + x**4 + alpha4/(x**2 - alpha0) + alpha5*x**2/(x**2 - alpha0) == 0.

# test 1
# solve(eq, x, ring=CC)
# ==> [0 == 20000*x^6 - 115000*x^4 + 138564*x^3 - 32500*x^2 - 34641*x + 22500]

# test 2
# solve(eq, x, ring=CC, solution_dict=True)
# ==> [{0: 20000*x^6 - 115000*x^4 + 138564*x^3 - 32500*x^2 - 34641*x + 22500}]

# test 3
# eq.roots(x, ring=CC,multiplicities=False)
# ==> TypeError: fraction must have unit denominator

Do you have any idea of a method to get the approximated roots of the equation ?

Thanks in advance

EDIT : correction of an error in the equation ; add few tests

edit retag flag offensive close merge delete

add a comment

3 Answers

Sort by » oldest newest most voted

answered 2017-12-26 23:08:09 +0200

Emmanuel Charpentier
7057 ●6 ●46 ●137

updated 2017-12-26 23:11:19 +0200

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) # The conversion method is highly unorthodox !
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,

edit flag offensive delete link

add a comment

answered 2017-12-27 17:39:46 +0200

dan_fulea
5337 ●5 ●43 ●90

Let us use exact arithmetics (first), instead of numerics. (Using an inaccurate f below also works, when we only needs roots in ring=CC, which uses a given precision, usually a sufficient one.)

alpha0 =  1/4             #  0.25
alpha1 = -5/2             # -2.5
alpha2 =  69282 / 10000   #  6.9282
alpha3 = -55/10           # -5.5
alpha4 =  1/2             #  0.5
alpha5 = -1/2             # -0.5

x = var('x')
f = alpha1 + alpha2*x + alpha3*x**2 + x**4 + alpha4/(x**2 - alpha0) + alpha5*x**2/(x**2 - alpha0)
g = f.factor().numerator()
print "g = %s" % g

I have only invented g. Then one of the following is the solution:

g.roots( ring=QQbar, multiplicities=False )
g.roots( ring=CC   , multiplicities=False )

For instance:

print "g = %s" % g
print "The roots of g are as follows:"
g.roots( ring=QQbar, multiplicities=False )

Results:

g = 20000*x^6 - 115000*x^4 + 138564*x^3 - 32500*x^2 - 34641*x + 22500
The roots of g are as follows:
[-2.875215523839503?,
 -0.5433030528116862?,
 0.9541487172225938?,
 1.408897521688209?,
 0.5277361688701930? - 0.5071712949223566?*I,
 0.5277361688701930? + 0.5071712949223566?*I]

edit flag offensive delete link

Comments

Thanks for this code. This solution also works for me !

yoda789 ( 2018-01-01 10:23:23 +0200 )edit

add a comment

answered 2017-12-27 10:45:08 +0200

yoda789
13 ●1 ●1 ●5

updated 2017-12-27 10:52:22 +0200

Thank you very much ! The numerical resolution works perfectly ! But notice that you have forgotten "t^4" term in F definition. Then, the solving of the equation give 6 roots

{-2.8752155238395045 - 1.6035277768266649e-16*I,
 -0.5433030528116871 + 2.5081428804803407e-17*I,
 0.5277361688701939 - 0.5071712949223567*I,
 0.5277361688701941 + 0.5071712949223569*I,
 0.9541487172225926 + 1.4813528249977233e-16*I,
 1.408897521688212 + 3.3107296343907337e-16*I}

which are identical to the ones given by Mathematica.

{{x -> -2.87522},{x -> -0.543303},{x -> 0.527736-0.507171i},{x -> 0.527736+0.507171i},{x -> 0.95415},{x -> 1.4089}}

This suits me.

However, for information, I would to mention that I've checked the following code on https://sagecell.sagemath.org/ :

P.<t>=PolynomialRing(CDF);
FF=FractionField(P);
alpha=[CDF(complex(2*(random()-0.5),2*(random()-0.5))) for p in range(6)];

F=alpha[1]+alpha[2]*t+alpha[3]*t^2+t^4+(alpha[4]+alpha[5]*t^2)/(t^2-alpha[0]);

foo=SR(repr(F).replace("t","x")).solve(x)
[q.rhs().n() for q in foo]

and I've got the following errors :

---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
<ipython-input-1-a2132da36519> in <module>()
     19 
     20 foo=SR(repr(F).replace("t","x")).solve(x)
---> 21 [q.rhs().n() for q in foo]
     22 

/home/sc_serv/sage/src/sage/structure/element.pyx in sage.structure.element.Element.n (build/cythonized/sage/structure/element.c:8063)()
    861             0.666666666666667
    862         """
--> 863         return self.numerical_approx(prec, digits, algorithm)
    864 
    865     N = deprecated_function_alias(13055, n)

/home/sc_serv/sage/src/sage/symbolic/expression.pyx in sage.symbolic.expression.Expression.numerical_approx (build/cythonized/sage/symbolic/expression.cpp:36129)()
   5782             res = x.pyobject()
   5783         else:
-> 5784             raise TypeError("cannot evaluate symbolic expression numerically")
   5785 
   5786         # Important -- the  we get might not be a valid output for numerical_approx in

TypeError: cannot evaluate symbolic expression numerically

edit flag offensive delete link

Comments

Better open a new question for your follow-up question. But anyway, here is an answer.

Does this achieve what you wanted?

sage: [CDF['x'](q.rhs()) for q in foo]
[(5.940048125501227e+69)*x^6
+ (5.898352846771207e+68 + 1.1746606221357097e+70*I)*x^4
+ (5.292143043372301e+68 - 4.858750452229818e+69*I)*x^3
+ (-7.85794262458788e+69 + 1.084454625704207e+70*I)*x^2
+ (4.9549885732219846e+69 - 7.388262794060325e+68*I)*x
- 9.419951737570192e+69 + 1.8520244744581162e+69*I]

slelievre ( 2017-12-28 03:20:54 +0200 )edit

Thanks ! This command works but it doesn't give the list of roots, unlike the code of Emmanuel Charpentier.

yoda789 ( 2018-01-01 10:22:48 +0200 )edit

How about the following then:

sage: [CDF['x'](q.rhs()).roots() for q in foo]

sage: [q.rhs().roots(CDF) for q in foo]

or even

sage: [q.rhs().roots(QQbar) for q in foo]

slelievre ( 2018-01-04 02:14:27 +0200 )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

Getting all (complex) solutions of a non polynomial equation

3 Answers

Comments

Comments

Your Answer

Question Tools

Stats

Related questions

Getting all (complex) solutions of a non polynomial equation edit

3 Answers

Comments

Comments

Your Answer

Question Tools

Stats

Related questions

Getting all (complex) solutions of a non polynomial equation