# Elimination of the variables

I am trying to get an equation in terms of variable z by eliminating variables g0 and g1 but

Running below code is taking time forever

In[10]:=
y = 0.01;
x = 0.057;
n0 = 0;
n1 = 1;

gE0 = 1/g0 - Sqrt[1 + z^2/(Pi*(2*n0 + 1)*y + 4.4*x*Pi*g0)^2];
gE1 = 1/g1 - Sqrt[1 + z^2/(Pi*(2*n1 + 1)*y + 4.4*x*Pi*g1)^2];
sce0 = 2*Pi*0.01*(((2*n0 + 1)*Pi*0.01*(g0 - 1))/(((2*n0 + 1)*Pi*0.01 + 4.4*Pi*x)*((2*n0 + 1)*Pi*0.01 + 4.4*Pi*x*g0))) + 0.6431; sce1 = 2*Pi*0.01*(((2*n1 + 1)*Pi*0.01*(g1 - 1))/(((2*n1 + 1)*Pi*0.01 + 4.4*Pi*x)*((2*n1 + 1)*Pi*0.01 + 4.4*Pi*x*g1)));

xed0 = Eliminate[{sce0 + sce1 == 0, gE0 == 0, gE1 == 0}, {g0, g1}]


Can someone simplify the above code to make it run quickly?

• @Bill as far as I can tell z is only part of gE0 and gE1 in the first place? I think sce0+sce1 == 0 simply serves to relate g0 and g1. – thorimur May 18 at 4:51

I think due to presence of radicals and denominators it is necessary to use exact (or at least high precision) arithmetic here.

y = 1/100;
x = 57/1000;
n0 = 0;
n1 = 1;
gE0 = 1/g0 - Sqrt[1 + z^2/(Pi*(2*n0 + 1)*y + 44/10*x*Pi*g0)^2];
gE1 = 1/g1 - Sqrt[1 + z^2/(Pi*(2*n1 + 1)*y + 44/10*x*Pi*g1)^2];
sce0 = 2*Pi*1/
100*(((2*n0 + 1)*
Pi*1/100*(g0 - 1))/(((2*n0 + 1)*Pi*1/100 +
44/10*Pi*x)*((2*n0 + 1)*Pi*1/100 + 44/10*Pi*x*g0))) +
6431/10000; sce1 =
2*Pi*1/100*(((2*n1 + 1)*
Pi*1/100*(g1 - 1))/(((2*n1 + 1)*Pi*1/100 +
44/10*Pi*x)*((2*n1 + 1)*Pi*1/100 + 44/10*Pi*x*g1)));
ee = {sce0 + sce1, gE0, gE1};


The best way to eliminate variables is (alas) not with Eliminate (old old technology) but rather with GroebnerBasis. This even has a term order specialized for the task (though it's not really needed for this example).

GroebnerBasis[ee, {z}, {g0, g1},
MonomialOrder -> EliminationOrder]

(* Out[61]= {\
1087191781255118852432205908237263695371306687898414793553937916276108\
8539220193254804961973869866459309790986240000000000000000000000000000\
0000 \[Pi]^12 -
16268517781507860032737188758040268612029306311347272704315305344023\
7238897918646667014719758990108604114670256499168765721664552960000000\
0000000000000 \[Pi]^10 z^2 +
56738570910688790592087477644695709810029509155667904340817842292123\
9347364218848001458326552347662710033232181322227278181276799464374272\
0000000000000000 \[Pi]^8 z^4 -
76141485891245220851706051052591264935429766097405661030697946539103\
8953552797608365697130284448476658681402247606314777047661399316777410\
9663232000000000000 \[Pi]^6 z^6 +
48091488360591912991048204087577606236749539610147285181149875603683\
1168149454360470756370124698044548746378946137263823892891719840553131\
1438872323891200000000 \[Pi]^4 z^8 -
13638620110911998467924632042373320416852081399506396254722968526570\
9692816709167930328435898547536877212793691508552018024025486557122137\
3396013486887615269120000 \[Pi]^2 z^10 +
12072464157868041559939145105552971974917159481422285017012503613471\
4944532397589113777181341371047147797682391535913140689377031298240042\
493086282606446471987952001 z^12} *)


In finite precision we see there is a scale difference in the coefficients of the result that might also play a role in messing up any machine arithmetic attempt.

GroebnerBasis[N[ee, 25], {z}, {g0, g1},
MonomialOrder -> EliminationOrder]

(* Out[63]= {8.3235522*10^-16 - 1.26197511*10^-9 z^2 +
4.4594515*10^-7 z^4 - 0.000060635178 z^6 + 0.0038803579 z^8 -
0.111499842 z^10 + 1.00000000 z^12} *)


At the slightly lower precision of 20 digits this fails, with a message that there was not sufficient precision for the computation. Usijng machine arithmetic it silently gives a bad result of {1.} (meaning the expressions are not consistent). They are of course, but there is some amount of unwinding to do in order to make polynomials out of the radicals and denominators, and that can make for trouble at finite precision.

• Actually it worked for n2 = 2 as well. Actually I am trying to solve for several n's equations self consistently and get the values for z where I need to have convergence after certain value of n. Is there any other way above problem can be performed so as to solve for value of z numerically? Already for nmax = 3 its running forever. – Tiku May 18 at 23:20
• I do not know what n2 or nmax are, or what would be the equations that are problematic. – Daniel Lichtblau May 19 at 14:12

Ok, I'm not sure why, but after rationalizing your equations, performing the elimination in sequence seemed to work:

y = Rationalize@0.01;
x = Rationalize@0.057;
n0 = 0;
n1 = 1;

gE0 = 1/g0 -
Sqrt[1 + z^2/(Pi*(2*n0 + 1)*y + Rationalize[4.4]*x*Pi*g0)^2];
gE1 = 1/g1 -
Sqrt[1 + z^2/(Pi*(2*n1 + 1)*y + Rationalize[4.4]*x*Pi*g1)^2];
sce0 = 2*Pi*(1/
100)*(((2*n0 + 1)*
Pi*(1/100)*(g0 - 1))/(((2*n0 + 1)*Pi*(1/100) +
Rationalize[4.4]*Pi*x)*((2*n0 + 1)*Pi/100 +
Rationalize[4.4]*Pi*x*g0))) + Rationalize[0.6431]; sce1 =
2*Pi*(1/100) (((2*n1 + 1)*
Pi/100*(g1 - 1))/(((2*n1 + 1)*Pi/100 +
Rationalize[4.4]*Pi*x)*((2*n1 + 1)*Pi/100 +
Rationalize[4.4]*Pi*x*g1)));

xed0 = Eliminate[Eliminate[{sce0 + sce1 == 0, gE0 == 0, gE1 == 0}, {g1}], {g0}]


I'm not totally sure if this is "right", though. I'd take it with a grain of salt, as some one-way implications might have been used; this might not be an equivalent condition.

• Interesting: changing the elimination order gives the same result, but doing both simultaneously does not terminate within say 1minute. May be worth reporting to Wolfram. – A.G. May 18 at 4:07

Using the equations with rationalized constants.

The elimination can be simplified by equating the numerators of the LHS of each equation to zero.

(xed0 = Eliminate[{
Numerator@Together[sce0 + sce1] == 0,
Numerator@Together@gE0 == 0,
Numerator@Together@gE1 == 0}, {g0, g1}]) // N

(* -1.52351*10^155 z^2 + 5.38366*10^157 z^4 - 7.32016*10^159 z^6 +
4.68455*10^161 z^8 - 1.34608*10^163 z^10 +
1.20725*10^164 z^12 == -1.00486*10^149 *)


This is the same as the result provided by thorimur

All of the roots are real

(solz = Solve[xed0, z]) // N

(* {{z -> -0.000812231}, {z -> 0.000812231}, {z -> -0.107171}, {z ->
0.107171}, {z -> -0.108051}, {z -> 0.108051}, {z -> -0.109268}, {z ->
0.109268}, {z -> -0.110874}, {z -> 0.110874}, {z -> -0.253192}, {z ->
0.253192}} *)


EDIT: To add additional equations use Daniel Lichtblaus answer. For example,

Clear["Global*"]

y = 1/100;
x = 57/1000;


Generalize the definitions

gE[n_Integer?NonNegative] :=
1/g[n] - Sqrt[1 + z^2/(Pi*(2*n + 1)*y + 44/10*x*Pi*g[n])^2];

sce[0] = 2*
Pi*1/100*((Pi*1/
100*(g[0] - 1))/((Pi*1/100 + 44/10*Pi*x)*(Pi*1/100 +
44/10*Pi*x*g[0]))) + 6431/10000;

sce[n_Integer?Positive] :=
2*Pi*1/100*(((2*n + 1)*
Pi*1/100*(g[n] - 1))/(((2*n + 1)*Pi*1/100 +
44/10*Pi*x)*((2*n + 1)*Pi*1/100 + 44/10*Pi*x*g[n])));

nmax = 2;

nValues = Range[0, nmax];

ee = {Total[sce /@ nValues], gE /@ nValues} // Flatten;

xed0 = GroebnerBasis[ee, {z}, g /@ nValues,
MonomialOrder -> EliminationOrder];


Solving for z

(solz = Solve[xed0 == 0, z]) // N

(* {{z -> -0.000761198}, {z -> 0.000761198}, {z -> -0.000812231}, {z ->
0.000812231}, {z -> -0.00159705}, {z ->
0.00159705}, {z -> -0.00161735}, {z ->
0.00161735}, {z -> -0.00234828}, {z ->
0.00234828}, {z -> -0.00243983}, {z ->
0.00243983}, {z -> -0.101386}, {z ->
0.101386}, {z -> -0.102136}, {z ->
0.102136}, {z -> -0.102655}, {z ->
0.102655}, {z -> -0.10322}, {z -> 0.10322}, {z -> -0.103766}, {z ->
0.103766}, {z -> -0.104686}, {z ->
0.104686}, {z -> -0.106111}, {z ->
0.106111}, {z -> -0.107184}, {z ->
0.107184}, {z -> -0.107299}, {z ->
0.107299}, {z -> -0.108065}, {z ->
0.108065}, {z -> -0.109282}, {z ->
0.109282}, {z -> -0.110889}, {z -> 0.110889}, {z -> -0.1439}, {z ->
0.1439}, {z -> -0.145999}, {z -> 0.145999}, {z -> -0.15033}, {z ->
0.15033}, {z -> -0.236372}, {z ->
0.236372}, {z -> -0.240141}, {z ->
0.240141}, {z -> -0.243525}, {z ->
0.243525}, {z -> -0.24677}, {z -> 0.24677}, {z -> -0.250029}, {z ->
0.250029}, {z -> -0.253403}, {z ->
0.253403}, {z -> -0.438011}, {z -> 0.438011}} *)


As expected, the complexity grows with nmax

• Thank you for the reply and it worked for n0=0 , n1=1 but if I add one more equation for n2=2 and repeat the same process its taking forever. In fact I have to vary n(let's say from n=0 -10) and solve for z. – Tiku May 18 at 13:51
• What equation are you adding? We can only work with the information that you provide. – Bob Hanlon May 18 at 14:30
• I am adding gE2 = 1/g2 - Sqrt[1 + z^2/(Pi*(2*n2 + 1)*y + Rationalize[4.4]*xPig2)^2]; sce2 = 2*Pi*(1/100)*(((2*n2 + 1)*(Pi/100)*(g2 - 1))/(((2*n2 + 1)*(Pi/100) + Rationalize[4.4]*Pix)*((2*n2 + 1)*(Pi/100) + Rationalize[4.4]*Pix*g2))); xed0 = N[Eliminate[{Numerator[Together[sce0 + sce1 + sce2]] == 0, Numerator[Together[gE0]] == 0, Numerator[Together[gE1]] == 0, Numerator[Together[gE2]] == 0}, {g0, g1, g2}]]; and other equations as above and this time I have to eliminate g2 as well – Tiku May 18 at 15:42