Using roots to get original equation

I have a list of four rules, that I would like to use to find the original polynomial that produced them.

The rules are of the form {{y -> ...}, {y -> ...}, {y -> ...}, {y -> ...}} and the equation is 0 == A y^4 +B y^3 + C y^2 + D y + E.

I am interested in finding the coefficients of the polynomial. The rules are complex, but the coefficients should not be.

I have tried using SolveAlways, but I don't believe I am using it correctly.

• You are probably aware of this, but it bears repeating anyway. If you only know the roots, at best you can recover a family of polynomials, all multiples of each other: if $a\ y^2+b\ y=0$, then of course $c\ (a\ y^2+b\ y)=0$ as well (for $c\neq0$). Mar 22 '16 at 21:28
• @user38729, you may use Rule -> Subtract to massage the rules. Mar 22 '16 at 21:39
• Interesting question. Actually, made me think. Does such a polynomial with real coefficients always exist? I mean, given complex roots can we always find a polynomial with real coefficients? Mar 22 '16 at 23:52
• @BehzadNazari so long as the complex numbers are in conjugate pairs there is a real-coefficient polynomial with those roots. Mar 23 '16 at 1:43
• @george2079 Yeah, OK now it makes sense and based on that condition I can see my answer always giving real coefficients as well. Mar 23 '16 at 1:49

poly = Total[MapIndexed[#1  x^(First@#2 - 1 ) &, RandomReal[ 1, {5}] ]]

0.310085 + 0.86543 x + 0.199379 x^2 + 0.630797 x^3 + 0.998793 x^4

roots = Solve[poly == 0, x]

{{x -> -0.992501}, {x -> -0.375095}, {x -> 0.368018 - 0.835763 I}, {x -> 0.368018 + 0.835763 I}}

result = Times @@ ((x - y) /. roots) // FullSimplify // Expand

0.310459 + 0.866476 y + 0.19962 y^2 + 0.63156 y^3 + y^4

Of course given only the roots, the polynomial can only be recovered within a multiplicative factor, but we can see this is the same:

Expand[result (poly /. x -> 0)/(result /. y -> 0)]

0.310085 + 0.86543 y + 0.199379 y^2 + 0.630797 y^3 + 0.998793 y^4

• Alternatively: Expand[First[Apply[Times, roots /. Rule -> Subtract]]] Mar 22 '16 at 22:24

As an alternative, check the method used in the help page of InterpolatingPolynomial (Link1). Assuming your rules are defined:

rules = {{y -> -1}, {y -> 2}, {y -> -3}, {y -> 4}}

Riffle with zeros (the function's value at the roots), add 0,0,1 to the end of the list and partition 2 by 2 to get the coordinates:

data = Partition[Join[Riffle[y /. rules, {0}], {0, 0, 1}], 2]

Then use InterpolatingPolynomial:

f[x_] := InterpolatingPolynomial[data, x]

which, for this example, results in:

1/24 (-4 + x) (-2 + x) (1 + x) (3 + x)

after expanding with f[x] // Expand :

1 + (7 x)/12 - (13 x^2)/24 - x^3/12 + x^4/24

Plotting this function shows that it passes through all the roots: Note that you need 5 points to get the exact interpolating polynomial of order 4. Your system has "an order of freedom" (you have four rules/equations but 5 unknowns/coefficients). I added point (0,1) to make it a closed system.

An example for complex roots in conjugate pairs:

rules = {{y -> -1 + 2 I}, {y -> -1 - 2 I}, {y -> -3 I}, {y -> +3 I}}
data = Partition[Join[Riffle[y /. rules, {0}], {0, 0, 1}], 2]
f[x_] := InterpolatingPolynomial[data, x]

results in

1 + (2 x)/5 + (14 x^2)/45 + (2 x^3)/45 + x^4/45

with the coefficients being:

CoefficientList[g = f[x], x]

resulting in

{1, 2/5, 14/45, 2/45, 1/45}

Make up some solutions:

SeedRandom["MMA.SE"]
sols = {y -> #} & /@ RandomInteger[{-10, 10}, 4]

(*Out: {{y -> 5}, {y -> 6}, {y -> -3}, {y -> -7}} *)

Then use them to recover the coefficients of the relative polynomial in $y$ (minus a factor, of course):

CoefficientList[Times @@ ((y - #) & /@ (y /. sols)), y]

(* Out: {630, 69, -59, -1, 1} *)

Notice that CoefficientList returns the coefficients in order of increasing power of $y$, from $y^0$ on up, so the polynomial represented above would be: $y^4-y^3-59\ y^2+69\ y+630$:

Solve[y^4 - y^3 - 59 y^2 + 69 y + 630 == 0, y]
(* Out: {{y -> -7}, {y -> -3}, {y -> 5}, {y -> 6}} *)

Let's do this for $n=3$, you can always assume that $A=1$:

f[y_] := y^2 + B y + C;
roots = {{y -> 1 - I}, {y -> 1 + I}};
S = (f[y] == 0) /. roots
Solve[S, {B, C}] // Simplify

Alternatively you can multiply $(y-y_1)..(y-y_n)$ as george2079 suggests.