# How to solve coupled multi-variable polynomials?

The following code generates two polynomials $q_1$ and $q_2$ in complex variables p and c:

pp[n_] := If[n > 0, pp[n - 1]^2 + c, p];
q1 = PolynomialQuotient[pp[4] - p, pp[2] - p, p];
q2 = PolynomialRemainder[D[pp[4], p] - 1, q1, p];


I would like to eliminate p and generate either a single polynomial in c or complex numerical solutions to c. Eliminate[] works for similar problems of lower order but never returns a result for this one.

Eliminate[{q1 == 0, q2 == 0}, p]


(* no result returned *)

If I get a solution to that problem, the next thing I need to do is generalize the "1" in the $q_2$ line to a variable t, so I can calculate c[t] where Abs[t]==1.

q2 = PolynomialRemainder[D[pp[4, c], p] - t, q1, p];


Thanks in advance for any help.

UPDATE: In the course of trying proposed solutions to this, I came upon what appears to be a bug, discussed HERE

• GroebnerBasis[{q1,q2},{p,c}] gives a degree 4 polynomial in p, indicating that for each value of c, viewing c as a parameter, there are 4 solutions. The basis also contains a degree 6 polynomial in c, limiting its values to the 6 roots of this polynomial. – John McGee Dec 23 '15 at 23:55
• @JohnMcGee Thanks for the help. The 4 solutions part is great news, because I suspected as much. I didn't understand your second sentence though. – Jerry Guern Dec 24 '15 at 3:50
• GroebnerBasis[{q1, q2}, {c}, {p}, MonomialOrder -> EliminationOrder, CoefficientDomain -> RationalFunctions] (* Out[315]= {4096 + 12288 c^5 + 4096 c^6 - 768 t + 48 t^2 - t^3 + c^3 (12288 + 256 t) + c^4 (12288 + 256 t) + c^2 (8192 - 256 t - 16 t^2)} *) – Daniel Lichtblau Dec 28 '15 at 22:32

## Part 1 - Eliminate p

I am not sure this is what you are looking for but both Solve and Reduce give solutions.

When I evaluate q1 and q2 I get

pp[n_Integer] := If[n > 0, pp[n - 1]^2 + c, p]

q1 = PolynomialQuotient[pp[4] - p, pp[2] - p, p]
(* 1 + 2 c^2 + 3 c^3 + 3 c^4 +
3 c^5 + c^6 + (2 c + c^2 + 2 c^3 + c^4) p + (c + 5 c^2 + 6 c^3 +
12 c^4 + 6 c^5) p^2 + (1 + 4 c^2 + 4 c^3) p^3 + (4 c + 3 c^2 +
18 c^3 + 15 c^4) p^4 + (2 c + 6 c^2) p^5 + (1 + 12 c^2 +
20 c^3) p^6 + 4 c p^7 + (3 c + 15 c^2) p^8 + p^9 + 6 c p^10 + p^12 *)

q2 = PolynomialRemainder[D[pp[4], p] - 1, q1, p]
(* 15 + 32 c^2 + 48 c^3 + 48 c^4 + 48 c^5 +
16 c^6 + (16 c + 16 c^2 + 16 c^3) p + (16 c + 48 c^2 + 80 c^3 +
160 c^4 + 80 c^5) p^2 +
32 c^2 p^3 + (16 c + 32 c^2 + 192 c^3 + 160 c^4) p^4 +
16 c p^5 + (96 c^2 + 160 c^3) p^6 + (16 c + 80 c^2) p^8 + 16 c p^10 *)


Now apply Solve

solution =
Solve[{1 + 2 c^2 + 3 c^3 + 3 c^4 + 3 c^5 +
c^6 + (2 c + c^2 + 2 c^3 + c^4) p + (c + 5 c^2 + 6 c^3 +
12 c^4 + 6 c^5) p^2 + (1 + 4 c^2 + 4 c^3) p^3 + (4 c +
3 c^2 + 18 c^3 + 15 c^4) p^4 + (2 c + 6 c^2) p^5 + (1 +
12 c^2 + 20 c^3) p^6 + 4 c p^7 + (3 c + 15 c^2) p^8 + p^9 +
6 c p^10 + p^12 == 0,
15 + 32 c^2 + 48 c^3 + 48 c^4 + 48 c^5 +
16 c^6 + (16 c + 16 c^2 + 16 c^3) p + (16 c + 48 c^2 + 80 c^3 +
160 c^4 + 80 c^5) p^2 +
32 c^2 p^3 + (16 c + 32 c^2 + 192 c^3 + 160 c^4) p^4 +
16 c p^5 + (96 c^2 + 160 c^3) p^6 + (16 c + 80 c^2) p^8 +
16 c p^10 == 0}, {c, p}, Complexes];


There is a long list of solutions (some of which have Root) but I won't use up the space showing them. Here is the first:

solution[[1]]
(* {c -> 1/4 - I/2, p -> -(I/2)} *)


## Part 2 - Introduce variable t

q2 is changed to include a variable t in place of the 1.

q2 = PolynomialRemainder[D[pp[4], p] - t, q1, p]
(* 16 + 32 c^2 + 48 c^3 + 48 c^4 + 48 c^5 +
16 c^6 + (16 c + 16 c^2 + 16 c^3) p + (16 c + 48 c^2 + 80 c^3 +
160 c^4 + 80 c^5) p^2 +
32 c^2 p^3 + (16 c + 32 c^2 + 192 c^3 + 160 c^4) p^4 +
16 c p^5 + (96 c^2 + 160 c^3) p^6 + (16 c + 80 c^2) p^8 + 16 c p^10 -
t *)


The expression handed to Solve is modified accordingly

solT = Solve[{1 + 2 c^2 + 3 c^3 + 3 c^4 + 3 c^5 +
c^6 + (2 c + c^2 + 2 c^3 + c^4) p + (c + 5 c^2 + 6 c^3 +
12 c^4 + 6 c^5) p^2 + (1 + 4 c^2 + 4 c^3) p^3 + (4 c +
3 c^2 + 18 c^3 + 15 c^4) p^4 + (2 c + 6 c^2) p^5 + (1 +
12 c^2 + 20 c^3) p^6 + 4 c p^7 + (3 c + 15 c^2) p^8 + p^9 +
6 c p^10 + p^12 == 0,
16 + 32 c^2 + 48 c^3 + 48 c^4 + 48 c^5 +
16 c^6 + (16 c + 16 c^2 + 16 c^3) p + (16 c + 48 c^2 + 80 c^3 +
160 c^4 + 80 c^5) p^2 +
32 c^2 p^3 + (16 c + 32 c^2 + 192 c^3 + 160 c^4) p^4 +
16 c p^5 + (96 c^2 + 160 c^3) p^6 + (16 c + 80 c^2) p^8 +
16 c p^10 == t}, {c, p}, Complexes];


Again there are many solutions.

I don't believe there is any hope to create a functionc[t] for that can be used as an input to Plot.

The first solution was used to form a table of values that could be input to ListPlot.

First a table of complex numbers whose absolute value is unity is formed.

tList = Table[k + I Sqrt[1 - k^2], {k, -1, 1, 0.1}]
(* {-1. + 0. I, -0.9 + 0.43589 I, -0.8 + 0.6 I, -0.7 +
0.714143 I, -0.6 + 0.8 I, -0.5 + 0.866025 I, -0.4 +
0.916515 I, -0.3 + 0.953939 I, -0.2 + 0.979796 I, -0.1 + 0.994987 I,
0. + 1. I, 0.1 + 0.994987 I, 0.2 + 0.979796 I, 0.3 + 0.953939 I,
0.4 + 0.916515 I, 0.5 + 0.866025 I, 0.6 + 0.8 I, 0.7 + 0.714143 I,
0.8 + 0.6 I, 0.9 + 0.43589 I, 1. + 0. I} *)


This is used as the input t to get a table of c[t].

solT1 = Table[{t, N[solT[[1, 1, 2]], 30]}, {t, tList}]
(* {{-1. + 0. I, -1.93971 + 0. I}, {-0.9 +
0.43589 I, -1.94168 + 0.000929621 I}, {-0.8 + 0.6 I, -1.94071 -
0.0000163376 I}, {-0.7 + 0.714143 I, -1.94177 +
0.000473936 I}, {-0.6 + 0.8 I, -1.94109 + 0.000180776 I}, {-0.5 +
0.866025 I, -1.94149 + 0.00108421 I}, {-0.4 +
0.916515 I, -1.93945 - 0.000428769 I}, {-0.3 +
0.953939 I, -1.93898 - 0.000332329 I}, {-0.2 +
0.979796 I, -1.93938 + 0.000494942 I}, {-0.1 +
0.994987 I, -1.94207 - 0.000668456 I}, {0. + 1. I, -1.93976 +
0.000452186 I}, {0.1 + 0.994987 I, -1.94271 -
0.00109377 I}, {0.2 + 0.979796 I, -1.9416 + 0.00229018 I}, {0.3 +
0.953939 I, -1.93757 + 0.00229704 I}, {0.4 +
0.916515 I, -1.94232 + 0.00142531 I}, {0.5 +
0.866025 I, -1.93943 + 0.00165118 I}, {0.6 + 0.8 I, -1.94212 -
0.00232099 I}, {0.7 + 0.714143 I, -1.9319 + 0.000202267 I}, {0.8 +
0.6 I, -1.92091 + 0.00870952 I}, {0.9 + 0.43589 I, -1.93432 +
0.00732884 I}, {1. + 0. I, -1.92309 + 0. I}} *)


I don't know how your want to plot the two complex numbers (t vs c[t]).

I chose to plot the real part of t on the x - axis and make two plots with the imaginary part of c[t] on the y -axis.

ListLinePlot[Transpose[{Re[solT1[[All, 1]]], Im[cList[[All, 2]]]}]]


ListLinePlot[Transpose[{Re[solT1[[All, 1]]], Re[solT1[[All, 2]]]}]]


• Oh, I was trying to eliminate p and didn't think of just solving for {c,p} as a pair. Any thoughts on the generalization part, so I can plot c[t] where Abs[t] in some reasonable time? – Jerry Guern Dec 24 '15 at 0:26
• It also really surprised me that this returns complex roots that aren't paired with their complex conjugate. I know that if a single-variable polynomial has real coefficients, it's complex roots will be paired with their conjugate. Was I confused when I expected that here? – Jerry Guern Dec 24 '15 at 0:27
• Actually, this is returning incorrect results. I posted a Question about it. – Jerry Guern Dec 24 '15 at 3:49
• I think you meant replace 1 in q1 with t rather than q2. I don't see a 1 in q2. I tried replacing 1 in q1 with t and Solve spit out a result but I don't think there is any hope of being able to plot it. – Jack LaVigne Dec 24 '15 at 15:42
• Your posted question addresses NSolve rather than Solve. Solve is giving the correct answer. If you numerically evaluate them some answers appear wrong. In particular solution[[5]]. {q1, q2} /. N[solution[[5]]] is bad. It turns out that this is a matter of precision. {q1, q2} /. N[solution[[5]], 30] gives the correct solution. – Jack LaVigne Dec 24 '15 at 15:45