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?

=================

1

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.

– J. M.♦

May 21 ’15 at 17:18

1

(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.

– Daniel Lichtblau

May 21 ’15 at 19:39

@Guesswhoitis. I changed the code to remove the inequalities, ran into the same problem.

– Jerry Guern

May 21 ’15 at 22:06

@DanielLichtblau I changed the code to remove the inequalities, ran into the same problem.

– Jerry Guern

May 21 ’15 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?

– Bichoy

May 22 ’15 at 2:47

=================

1 Answer

1

=================

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?

– J. M.♦

May 22 ’15 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.

– Daniel Lichtblau

May 22 ’15 at 1:50