3
$\begingroup$

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

$\endgroup$
  • $\begingroup$ 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. $\endgroup$ – John McGee Dec 23 '15 at 23:55
  • $\begingroup$ @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. $\endgroup$ – Jerry Guern Dec 24 '15 at 3:50
  • $\begingroup$ 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)} *) $\endgroup$ – Daniel Lichtblau Dec 28 '15 at 22:32
4
$\begingroup$

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]]]}]]

Mathematica graphics

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

Mathematica graphics

|improve this answer|||||
$\endgroup$
  • $\begingroup$ 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? $\endgroup$ – Jerry Guern Dec 24 '15 at 0:26
  • $\begingroup$ 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? $\endgroup$ – Jerry Guern Dec 24 '15 at 0:27
  • $\begingroup$ Actually, this is returning incorrect results. I posted a Question about it. $\endgroup$ – Jerry Guern Dec 24 '15 at 3:49
  • $\begingroup$ 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. $\endgroup$ – Jack LaVigne Dec 24 '15 at 15:42
  • $\begingroup$ 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. $\endgroup$ – Jack LaVigne Dec 24 '15 at 15:45

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.