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
Check the numerical answers:
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
Improve readibility
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$$
Eliminate radicals
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()
).
Please take the time to type your questions in a machine-readable for : your equations are full of implicit products (usually a really bad idea...).
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.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!Thanks Emmanuel and SLelievre for your quick answer.
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,
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,