# How to speed or help Solve[] when it is stuck?

The code below generates a system of equations that I need to Solve[] for {c[5], c[6], c[7], c[8]} in terms of {c[1], c[2], c[3], c[4]}.

f[x_] := E^-Sum[a[i]*x^i, {i, 0, 4}];
zeroEq = {};
Do[Do[
zeroEq = Append[zeroEq, Expand[D[f[x]*x^n, {x, m}]/f[x]] /. x^r_ -> c[r] /. x -> c[1]],
{n, 0, 8 - 3*m}], {m, 1, 2}];
aElim = Eliminate[zeroEq == 0, {a[1], a[2], a[3], a[4]}];


This much works fine, but then Solve[] just runs indefinitely without returning results or error messages.

Solve[aElim, {c[5], c[6], c[7], c[8]}]


I attempted applying some known inequalities to reduce the work on Solve[], but Reduce[] bogs down in the same way, no error message or result:

Reduce[aElim && c[2] > 0 && c[4] > 0 && c[6] > 0 && c[8] > 0]


So...

1) what is going on here?
2) Can I fix it?

• What happens if you don't include the inequalities? Also, it looks to me that you are equating an entire list (eqn) to 0; without Thread[], I can imagine how it might cause some trouble. Commented May 21, 2015 at 17:18
• (1) Solve is fine with Solve[list==0] and just quietly thread over one level. (2) I suspect Eliminate is not happy with inequalities though. Commented May 21, 2015 at 19:39
• @Guesswhoitis. I changed the code to remove the inequalities, ran into the same problem. Commented May 21, 2015 at 22:06
• @DanielLichtblau I changed the code to remove the inequalities, ran into the same problem. Commented May 21, 2015 at 22:07
• Is the solution going to be used for numerical calculations in Mathematica afterwards? Or you are more interested in getting a formal solution to the problem? Commented May 22, 2015 at 2:47

One can use GroebnerBasis to eliminate variables, and it is set up in a way that tends to be more efficient than Eliminate (which may be using some dated technology).

polys = {-a[1] - 2 a[2] c[1] - 3 a[3] c[2] - 4 a[4] c[3],
1 - a[1] c[1] - 2 a[2] c[2] - 3 a[3] c[3] - 4 a[4] c[4],
2 c[1] - a[1] c[2] - 2 a[2] c[3] - 3 a[3] c[4] - 4 a[4] c[5],
3 c[2] - a[1] c[3] - 2 a[2] c[4] - 3 a[3] c[5] - 4 a[4] c[6],
4 c[3] - a[1] c[4] - 2 a[2] c[5] - 3 a[3] c[6] - 4 a[4] c[7],
5 c[4] - a[1] c[5] - 2 a[2] c[6] - 3 a[3] c[7] - 4 a[4] c[8],
a[1]^2 - 2 a[2] + 4 a[1] a[2] c[1] - 6 a[3] c[1] + 4 a[2]^2 c[2] +
6 a[1] a[3] c[2] - 12 a[4] c[2] + 12 a[2] a[3] c[3] +
8 a[1] a[4] c[3] + 9 a[3]^2 c[4] + 16 a[2] a[4] c[4] +
24 a[3] a[4] c[5] + 16 a[4]^2 c[6], -2 a[1] + a[1]^2 c[1] -
6 a[2] c[1] + 4 a[1] a[2] c[2] - 12 a[3] c[2] + 4 a[2]^2 c[3] +
6 a[1] a[3] c[3] - 20 a[4] c[3] + 12 a[2] a[3] c[4] +
8 a[1] a[4] c[4] + 9 a[3]^2 c[5] + 16 a[2] a[4] c[5] +
24 a[3] a[4] c[6] + 16 a[4]^2 c[7],
2 - 4 a[1] c[1] + a[1]^2 c[2] - 10 a[2] c[2] + 4 a[1] a[2] c[3] -
18 a[3] c[3] + 4 a[2]^2 c[4] + 6 a[1] a[3] c[4] - 28 a[4] c[4] +
12 a[2] a[3] c[5] + 8 a[1] a[4] c[5] + 9 a[3]^2 c[6] +
16 a[2] a[4] c[6] + 24 a[3] a[4] c[7] + 16 a[4]^2 c[8]};
avars = {a[1], a[2], a[3], a[4]};
cvars = {c[1], c[2], c[3], c[4], c[5], c[6], c[7], c[8]};

Timing[
gb = GroebnerBasis[polys, cvars, avars,
MonomialOrder -> EliminationOrder, Sort -> True];]

(* Out[199]= {0., Null} *)


It's not giving a small result though.

In[202]:= gb // LeafCount

(* Out[202]= 67818 *)

In[203]:= Length[gb]

(* Out[203]= 15 *)


So further processing e.g. in Reduce might still bog down.

• I had gotten into the habit of using GroebnerBasis[] to build implicit equations from the parametric representations of algebraic curves, but had always been wondering "shouldn't I be able to use Eliminate[] for this?". Now, you mention that Eliminate[] "may be using some dated technology"; I guess it was difficult to get the Gröbner stuff to play nice with everything else? Commented May 22, 2015 at 0:30
• @Guess who it is I don't actually know how much trouble it is, or even what exactly Eliminate is now doing. It might be that it was never overhauled in the way that Solve was. Or maybe it was, and is falling short for other reasons. Commented May 22, 2015 at 1:50