# Solve a system of polynomial equations - getting different results using Symbolic Math Toolbox of MATLAB and Sage

I am trying to solve the following system of equations:

eq1 = 0.8 - 1.8352657742116e-5*sqrt(x1*(3071.1916*x2 + (2.19736842105263e-5 *x3 - 0.378721571752863)*(3069.38078875823*x1 - 802.539479414006*x2 - 802.539479414006*x3 + 1187113.71566341))/(2.19736842105263e-5*x3 - 0.378721571752863))

eq2 = (-0.00117574932506308*x2*sqrt(-(0.039945461792316*x2 + 0.039945461792316*x3 - 71.5306505072653)*(3071.1916*x2 + 32.9605263157895*x3 - 568082.357629294)/(2.19736842105263e-5*x3 - 0.378721571752863)) + 1.8352657742116e-5*sqrt(x1*(3071.1916*x2 + (2.19736842105263e-5*x3 - 0.378721571752863)*(3069.38078875823*x1 - 802.539479414006*x2 - 802.539479414006*x3 + 1187113.71566341))/(2.19736842105263e-5*x3 - 0.378721571752863))*(x2 + x3))/(x2 + x3)

eq3 = (0.007585843252*x2*(x2 + x3) - 0.00117574932506308*x3*sqrt(-(0.039945461792316*x2 + 0.039945461792316*x3 - 71.5306505072653)*(3071.1916*x2 + 32.9605263157895*x3 - 568082.357629294)/(2.19736842105263e-5*x3 - 0.378721571752863))*(2.19736842105263e-5*x3 - 0.378721571752863) + (x2 + x3)*(2.19736842105263e-5*x3 - 0.378721571752863)*(-0.00198227251415259*x2 - 0.00198227251415259*x3 + 42.4521708776886))/((x2 + x3)*(2.19736842105263e-5*x3 - 0.378721571752863))


When I use Symbolic Math Toolbox in Matlab through the following command:

[x1_bar,x2_bar,x3_bar] = solve(eqs,[x1,x2,x3], 'MaxDegree',4);


we get the following results:

x1_bar
ans =

1.0e+04 *

-0.0116 + 0.0000i
0.5338 + 0.0000i
-0.0232 + 0.0000i
0.2669 + 0.0000i
-0.0047 + 0.0000i
1.3191 + 0.0000i
-0.0321 - 0.5445i
0.0007 - 0.0113i
-0.0321 + 0.5445i
0.0007 + 0.0113i
-0.0012 + 0.0000i
5.1569 + 0.0000i

x2_bar =
1.0e+03 *

1.9653 + 0.0000i
1.9653 + 0.0000i
0.2480 + 0.0000i
0.2480 + 0.0000i
-0.6495 + 0.0000i
-0.6495 + 0.0000i
-0.1794 + 0.2236i
-0.1794 + 0.2236i
-0.1794 - 0.2236i
-0.1794 - 0.2236i
-0.0393 + 0.0000i
-0.0393 + 0.0000i

x3_bar =
1.0e+04 *

-0.0174 + 0.0000i
-0.0174 + 0.0000i
0.6522 + 0.0000i
0.6522 + 0.0000i
4.8817 + 0.0000i
4.8817 + 0.0000i
0.2382 + 2.1159i
0.2382 + 2.1159i
0.2382 - 2.1159i
0.2382 - 2.1159i
1.7273 + 0.0000i
1.7273 + 0.0000i


When I try to solve this problem with Sage using:

S = solve([eq1,eq2,eq3],x1,x2,x3)


I obtain the following results:

[{x1: -321.3551367159828 - 5445.311876550251*I,
x2: -179.359654094372 - 223.5730259195566*I,
x3: 2381.999033439763 - 21159.38999334523*I},
{x1: 6.685918535164614 - 113.2918302698741*I,
x2: -179.359654094372 - 223.5730259195566*I,
x3: 2381.999033439763 - 21159.38999334523*I},
{x1: 6.685918535164623 + 113.2918302698741*I,
x2: -179.359654094372 + 223.5730259195566*I,
x3: 2381.999033439763 + 21159.38999334523*I},
{x1: -321.3551367159836 + 5445.311876550251*I,
x2: -179.359654094372 + 223.5730259195566*I,
x3: 2381.999033439763 + 21159.38999334523*I}]


I tried also to use the PolynomialRing but it is failing to solve the problem.

EDIT : Here are the raw equations with the original quantities:

eq1 =w_G_a_in - K_a*sqrt((x1*(R*T_a + g*L_a*M_G_a)*((g*L_a*x1)/V_a -
(g*L_r*(x2 + x3 - L_bh*S_bh*rho_L))/V_r - F_riser + (R*T_a*x1)/(V_a*M_G_a)
+ (R*T_r*rho_L*x2)/(M_G_r_t*(x3 - V_r*rho_L + L_bh*S_bh*rho_L))))/(R*T_a*V_a))

eq2 =K_a*sqrt((M_G_a*((g*L_a*x1)/V_a + (R*T_a*x1)/(V_a*M_G_a))*((g*L_a*x1)/
V_a - (g*L_r*(x2 + x3 - L_bh*S_bh*rho_L))/V_r - F_riser + (R*T_a*x1)/(V_a*M_G_a)
+ (R*T_r*x2)/(M_G_r_t*(L_bh*S_bh - V_r + x3/rho_L))))/(R*T_a))-(GOR*PI*(2*F_riser - P_r + g*L_bh*rho_L + (g*L_r*(x2 + x3 - L_bh*S_bh*rho_L))/V_r
- (R*T_r*x2)/(M_G_r_t*(L_bh*S_bh - V_r + x3/rho_L))))/(GOR + 1) - (K_r*u1*x2*
sqrt(-((P0 + (R*T_r*x2)/(M_G_r_t*(L_bh*S_bh - V_r + x3/rho_L)))*(x2 + x3 - L_bh*S_bh*rho_L))/V_r))/(x2 + x3)

eq3 =PI*(GOR/(GOR + 1) - 1)*(2*F_riser - P_r + g*L_bh*rho_L + (g*L_r*(x2 + x3
- L_bh*S_bh*rho_L))/V_r - (R*T_r*x2)/(M_G_r_t*(L_bh*S_bh - V_r + x3/rho_L)))
- (K_r*u1*x3*sqrt(-((P0 + (R*T_r*x2)/(M_G_r_t*(L_bh*S_bh - V_r + x3/rho_L)))*(x2 + x3 - L_bh*S_bh*rho_L))/V_r))/(x2 + x3)

edit retag close merge delete

Your equations are not polynomial. For exemple, the second has the form :

$$\frac{\sqrt{\frac{{\left(c_{18} x_{2} + c_{8} x_{3} + c_{13}\right)} {\left(c_{5} x_{2} + c_{5} x_{3} + c_{3}\right)}}{c_{15} x_{3} + c_{0}}} c_{6} x_{2} + c_{4} {\left(x_{2} + x_{3}\right)} \sqrt{\frac{{\left(c_{16} x_{1} + c_{7} x_{2} + c_{7} x_{3} + c_{10}\right)} {\left(c_{15} x_{3} + c_{0}\right)} + c_{18} x_{2}}{c_{15} x_{3} + c_{0}}}}{x_{2} + x_{3}}$$

However,such equations can be generalized to polynomials : see slelievre's answer to this question.

( 2020-08-26 08:58:23 -0600 )edit

The question was entered with machine readable input, but not formatted as code so the * for products were interpreted as markdown formatting for "italics". That's why it looked full of implicit products and italics!

( 2020-08-26 09:14:02 -0600 )edit

Really sorry for the formatting issue. I am new to sage and this is the first time I post a question here.

So after taking a quick look on the proposed solution by SLelievre, I guess we need to transform these equations to ploynomials and then use groebner_basis() to solve it rather than roots(). Am I right?

Thanks,

( 2020-08-26 12:09:14 -0600 )edit

We want the locus where all three expressions given by eq1, eq2, eq3 are zero.

We might as well remove denominators.

For the second equation, starting from the expression given by @Emmanuel Charpentier,

• this means $\sqrt{A} \cdot B + C \cdot \sqrt{D} = 0$
• in other words $\sqrt{A} \cdot B = -C \cdot \sqrt{D}$
• which is a subset of where $A \ge 0$ and $D \ge 0$ and $A \cdot B^2 = C^2 \cdot D$
( 2020-08-27 01:22:28 -0600 )edit

Sort by » oldest newest most voted

NOTE : This is a totally re-written answer : the previous version were garbled by too much edits and possibly a typo-induced error.

Original data taken from the question

## Sage

# Compute
%time SolsSM=solve([eq1, eq2, eq3], [x1, x2, x3], solution_dict=True)
print(SolsSM)
# Check
CheckSM=[[abs(v.subs(u)) for v in [eq1, eq2, eq3]] for u in SolsSM]
print(CheckSM)

CPU times: user 2min 41s, sys: 1.26 s, total: 2min 42s # Not *that* fast...
Wall time: 2min 4s
[{x1: -321.3551367159828 - 5445.311876550251*I, x2: -179.359654094372 - 223.5730259195566*I,
x3: 2381.999033439763 - 21159.38999334523*I},
{x1: 6.685918535164614 - 113.2918302698741*I, x2: -179.359654094372 - 223.5730259195566*I,
x3: 2381.999033439763 - 21159.38999334523*I},
{x1: 6.685918535164623 + 113.2918302698741*I, x2: -179.359654094372 + 223.5730259195566*I,
x3: 2381.999033439763 + 21159.38999334523*I},
{x1: -321.3551367159836 + 5445.311876550251*I, x2: -179.359654094372 + 223.5730259195566*I,
x3: 2381.999033439763 + 21159.38999334523*I}]

[[2.9139133938393855e-15, 2.853183107612774e-15, 9.826835244968086e-15],
[2.220446049250313e-16, 2.1154867529750935e-17, 9.826835244968086e-15],
[2.344603589928194e-16, 5.288716882437734e-17, 9.826835244968086e-15],
[3.3513599103077758e-15, 3.422400574978585e-15, 9.826835244968086e-15]]


Close enough for government's work...

## Matlab

Not having Matlab, I copied the question's data and checked them :

# Solution Matlab, telle qu'imprimée dans l'énoncé
x1bar=[ -0.0116 + 0.0000*I,  0.5338 + 0.0000*I, -0.0232 + 0.0000*I,
0.2669 + 0.0000*I, -0.0047 + 0.0000*I,  1.3191 + 0.0000*I,
-0.0321 - 0.5445*I,  0.0007 - 0.0113*I, -0.0321 + 0.5445*I,
0.0007 + 0.0113*I, -0.0012 + 0.0000*I,  5.1569 + 0.0000*I]
x2bar=[  1.9653 + 0.0000*I,  1.9653 + 0.0000*I,  0.2480 + 0.0000*I,
0.2480 + 0.0000*I, -0.6495 + 0.0000*I, -0.6495 + 0.0000*I,
-0.1794 + 0.2236*I, -0.1794 + 0.2236*I, -0.1794 - 0.2236*I,
-0.1794 - 0.2236*I, -0.0393 + 0.0000*I, -0.0393 + 0.0000*I]
x3bar=[ -0.0174 + 0.0000*I, -0.0174 + 0.0000*I,  0.6522 + 0.0000*I,
0.6522 + 0.0000*I,  4.8817 + 0.0000*I,  4.8817 + 0.0000*I,
0.2382 + 2.1159*I,  0.2382 + 2.1159*I,  0.2382 - 2.1159*I,
0.2382 - 2.1159*I,  1.7273 + 0.0000*I,  1.7273 + 0.0000*I]
SolsML=map(lambda u,v,w:{x1:1.0e4*u, x2:1.0e3*v, x3:1.0e4*w},
x1bar, x2bar, x3bar)
# Numerical check
CheckML=[[abs(v.subs(u)) for v in [eq1, eq2, eq3]] for u in SolsML]
print(CheckML)

[[0.000102245046124216, 0.0504273057558141, 0.00407651153817091],
[0.000145426550936412, 0.0501796341587534, 0.00407651153817091],
[0.000167832577447902, 0.0000380920578297709, 0.00350920387509823]
[0.000718604724398997, 0.000924529359676793, 0.00350920387509823],
[0.000590312714312535, 1.60057187570945, 120.256029930173],
[0.00152913372345476, 1.60151069671859, 120.256029930173],
[7.89390598709807, 7.89386741627798, 0.00223166033868648],
[1.10797685571357, 1.10813035022203, 0.00223166033868648],
[7.89390598709807, 7.89386741627798, 0.00223166033868648],
[1.10797685571357, 1.10813035022203, 0.00223166033868648],
[0.000561991548113761, 1.59838923451504, 702.100941033462],
[1.06578983072011, 2.66474105678326, 702.100941033462]]


The first four answers might be correct (the original solutions were printed with a limited precision) ; the eight last seem (widely) off the mark.

# Look for an "exact" answer

"Exact" has a special meaning here : the original question involves (a lot of) numerical floating-point constants, presumably approximations of some (unspecified) quantities. An "exact" solution implies treating them as "exact" decimal numbers, i. e.converting them to rationals.

## Transform the equations in rational fractions

In order to use Sage's polynomial handling, we have to convert those expressions to rational fractions, by isolating the terms involving square roots and squaring the equations. This is not an automated process ; in order to understand the structure of the expression, it is convenient to replace (temporarily) the numerical constants by symbols :

# Utility : get the set of numerical constants occurring in an expression
def GetNums(expr):
if expr.operator() is None:
if expr.is_numeric(): return set([expr])
return set()
return reduce(union, list(map(GetNums, expr.operands())))
# Set of the numerivcal constants in the three original equations
Consts=reduce(union, map(lambda u:set(GetNums(u)), [eq1, eq2, eq3]))
# Generate a set of symbolic constants...
SymConsts=[var("c_{}".format(u)) for u in range(len(Consts))]
# ... and build a dictionary
Dic=dict(zip(Consts,SymConsts))
# We do not want to loose the signs (represented by multiplication by -1)...
Dic.pop(-1)
# ...nor the square roots (represented by exponentiation to 1/2)
Dic.pop(1/2)


Now we can represent our original equations (without right-and side) as symbolic equations, much more readable :

Sys=[u.subs(Dic)==0 for u in [eq1, eq2, eq3]]


which gives us the following system :

$$c_{14} \sqrt{\frac{{\left({\left(c_{16} x_{1} + c_{7} x_{2} + c_{7} x_{3} + c_{10}\right)} {\left(c_{15} x_{3} + c_{0}\right)} + c_{18} x_{2}\right)} x_{1}}{c_{15} x_{3} + c_{0}}} + c_{1} = 0$$

$$\frac{\sqrt{\frac{{\left(c_{18} x_{2} + c_{8} x_{3} + c_{13}\right)} {\left(c_{5} x_{2} + c_{5} x_{3} + c_{3}\right)}}{c_{15} x_{3} + c_{0}}} c_{6} x_{2} + c_{4} {\left(x_{2} + x_{3}\right)} \sqrt{\frac{{\left(c_{16} x_{1} + c_{7} x_{2} + c_{7} x_{3} + c_{10}\right)} {\left(c_{15} x_{3} + c_{0}\right)} + c_{18} x_{2}}{c_{15} x_{3} + c_{0}}}}{x_{2} + x_{3}} = 0$$

$$\frac{{\left(c_{15} x_{3} + c_{0}\right)} \sqrt{\frac{{\left(c_{18} x_{2} + c_{8} x_{3} + c_{13}\right)} {\left(c_{5} x_{2} + c_{5} x_{3} + c_{3}\right)}}{c_{15} x_{3} + c_{0}}} c_{6} x_{3} + {\left(c_{12} x_{2} + c_{12} x_{3} + c_{11}\right)} {\left(c_{15} x_{3} + c_{0}\right)} {\left(x_{2} + x_{3}\right)} + c_{9} x_{2} + c_{9} x_{3}}{{\left(c_{15} x_{3} + c_{0}\right)} {\left(x_{2} + x_{3}\right)}} = 0$$

For this task, two "external" helpers are really useful :

• Maxima's part function is a great shortcut for expression surgery ;
• Sympy's simplify is often the fastest way to a concise rewriting of a complicated expression.

So :

import sympy
E1=(Sys[0]-Sys[0].maxima_methods().part(1,2))^2


$$\frac{{\left({\left(c_{16} x_{1} + c_{7} x_{2} + c_{7} x_{3} + c_{10}\right)} {\left(c_{15} x_{3} + c_{0}\right)} + c_{18} x_{2}\right)} c_{14}^{2} x_{1}}{c_{15} x_{3} + c_{0}} = c_{1}^{2}$$

E2=(sympy.sympify(Sys[1]-Sys[1].maxima_methods().part(1,1,1)/\
Sys[1].maxima_methods().part(1,2)).simplify()._sage_())^2


$$\frac{{\left({\left(c_{16} x_{1} + c_{7} x_{2} + c_{7} x_{3} + c_{10}\right)} {\left(c_{15} x_{3} + c_{0}\right)} + c_{18} x_{2}\right)} c_{4}^{2}}{c_{15} x_{3} + c_{0}} = \frac{{\left(c_{18} x_{2} + c_{8} x_{3} + c_{13}\right)} {\left(c_{5} x_{2} + c_{5} x_{3} + c_{3}\right)} c_{6}^{2} x_{2}^{2}}{{\left(c_{15} x_{3} + c_{0}\right)} {\left(x_{2} + x_{3}\right)}^{2}}$$

E3=(sympy.sympify(Sys[2]-Sys[2].maxima_methods().part(1,1,1)/\
Sys[2].maxima_methods().part(1,2)).simplify()._sage_())^2


$$\frac{{\left(c_{18} x_{2} + c_{8} x_{3} + c_{13}\right)} {\left(c_{5} x_{2} + c_{5} x_{3} + c_{3}\right)} c_{6}^{2} x_{3}^{2}}{{\left(c_{15} x_{3} + c_{0}\right)} {\left(x_{2} + x_{3}\right)}^{2}} = \frac{{\left({\left(c_{12} x_{2} + c_{12} x_{3} + c_{11}\right)} {\left(c_{15} x_{3} + c_{0}\right)} {\left(x_{2} + x_{3}\right)} + c_{9} x_{2} + c_{9} x_{3}\right)}^{2}}{{\left(c_{15} x_{3} + c_{0}\right)}^{2} {\left(x_{2} + x_{3}\right)}^{2}}$$

which can be rewritten as a system of expressions (equations with implicit null right-hand-side)

Sp=[sympy.sympify((u.lhs()-u.rhs()).factor()).simplify()._sage_().numerator() for u in [E1, E2, E3]]


(The LaTeX representation is too large for the available engine...)

## Solve the system of fractions

As far as I know, Sage does not have (built-in) facilities for solving rational fractions. So we will solve the fraction denominators, while keeping in mind to eliminate "solutions" where any of the denominators is null.

# Solution ring
R1=PolynomialRing(QQbar,"x1, x2, x3")
# Translation of our symbolic constants to *rationals*
Qid={Dic[u]:QQbar(QQ(u)) for u in Dic.keys()}
SP=[R1(u.subs(Qid)) for u in Sp]
# Ideal generated by our system.
J=R1.ideal(*SP)

J.dimension()
0


So far, so good... But there is a fly in the ointment :

J.variety(QQbar)


hasn't returned for > 8 hours... CPU usage shows that one process (Sage...) eats 100% of its CPU, RAM usage is absolutely stable. I begin to suspect a bug. Advice welcome...

## What should follow

• Check for null denominators for all candidate solutions (see above...).
• Check the proposed solutions against the original equations :

Our transformations may have introduced extraneous solutions ($x^2=a^2$ has more solutions than $x=a$). The simplest way is to check the solutions against the original equations (or, better, against a version of these equations where the numerical constants are replaced by their rational transcriptions, thus (theoretically) allowing to use is_null()).

more

I just edited my initial question with the raw eqautions.

( 2020-09-08 08:44:44 -0600 )edit

Waw...Thanks a lot Emmanual for all the efforts you are spending to help on this issue.

I have been trying your proposed solution but I am getting :

J.dimension()
1


I am probably doing something wrong but I will be posting the equations with the original quantities if this can be helpful to debug.

more