My main confusion is about the difference between the two code blocks at the end of this long spiel, but the spiel contains the code to create the polynomial if its coefficients are helpful for explaining the difference. Thanks in advance!
I have this polynomial, eu1-eu2, of degree n (where n is a multiple of 3 between 3 and 93, inclusive).
solutions = Simplify[Solve[eu1 == eu2]];
solutions = p /. solutions;
Select[solutions, Element[#, Reals] &]
and
p /. Simplify[Solve[eu1 == eu2, p, Reals]]
should return the same roots, from what I can tell. And they do, for lower values of n. But as n increases, at some point they start returning slightly different roots, and later on, the first code block returns a bunch of different roots while the second one usually returns one root repeated many times. As n increases, the second one's repeated root seems to be farther away from any of the roots gotten from the first one.
The first one probably gets the true roots, since the size of "solutions" (which includes complex roots) equals n. But then what number is the second one getting?
Edit: 30 is the smallest n for which the number of roots is the same but the roots themselves are slightly different. Here the first method produces {0.333794} whereas the second produces {0.333788}.
33 is the smallest n for which the number of roots is different for each method. Here the first method produces {0.251609, 0.251683, 0.273157, 0.274114, 0.3336} whereas the second produces {0.33363}.
Edit: The code that sets up the polynomial is as follows:
n = 93;
q1 = (n - Mod[n, 3])/3;
q2 = q1;
upperbound = 0.5;
(mat1 = Array[0, {n, 3}];
prob1 = Array[0, {n, 1}];
prob2 = Array[0, {n, 1}];
prob3 = Array[0, {n, 1}];
For[m = 0, m <= n - 1, m++,
If[m + 1 - q1 != 0,
prob1[[m + 1, 1]] = q1/(m + 1);
prob2[[m + 1,
1]] = ((q2 - (n - m - 1))/(m + 1 - q1))*(1 - q1/(m + 1));
prob3[[m + 1, 1]] = 1 - prob1[[m + 1, 1]] - prob2[[m + 1, 1]];
mat1[[m + 1]] = {prob1[[m + 1, 1]], prob2[[m + 1, 1]],
prob3[[m + 1, 1]]},
mat1[[m + 1]] = {0, 0, 0}]])
(mat = Array[0, {n, 3}];
For[m = 0, m <= n - 1, m++,
If[m + 1 <= q1,
mat[[m + 1]] = {1, 0, 0},
mat[[m + 1]] = mat1[[m + 1]]]])
(mat3 = Array[0, {n, 3}];
prob4 = Array[0, {n, 1}];
prob5 = Array[0, {n, 1}];
prob6 = Array[0, {n, 1}];
For[m = 0, m <= n - 1, m++,
If[m + 1 - q2 != 0,
prob4[[m + 1,
1]] = ((q1 - (n - m - 1))/(m + 1 - q2))*(1 - q2/(m + 1));
prob5[[m + 1, 1]] = q2/(m + 1);
prob6[[m + 1, 1]] = 1 - prob4[[m + 1, 1]] - prob5[[m + 1, 1]];
mat3[[m + 1]] = {prob4[[m + 1, 1]], prob5[[m + 1, 1]],
prob6[[m + 1, 1]]},
mat3[[m + 1]] = {0, 0, 0}]])
(mat2 = Array[0, {n, 3}];
For[m = 0, m <= n - 1, m++,
If[m + 1 <= q2,
mat2[[m + 1]] = {0, 1, 0},
mat2[[m + 1]] = mat3[[m + 1]]]])
u = {{1 - p}, {p}, {0}};
(eu1 = 0;
For[m = 0, m <= n - 1, m++,
eu1 = eu1 + (mat[[m + 1]].u)*(p/upperbound)^
m*(1 - p/upperbound)^(n - m - 1)])
eu1 = Simplify[eu1];
(eu2 = 0;
For[m = 0, m <= n - 1, m++,
eu2 = eu2 + (mat2[[m + 1]].u)*(p/upperbound)^(n - m - 1)*(1 -
p/upperbound)^m])
eu2 = Simplify[eu2];
Also note that the following code results in the same solutions as the first code block:
Simplify[Solve[eu1 - eu2 == 0 && p < 15000, p, Complexes]]
Whereas this code results in the same solutions as the second code block:
Simplify[Solve[eu1 - eu2 == 0 && p < 15000, p, Reals]]
So it's something about allowing the function values to be complex vs. not allowing them to be complex. But I don't understand why my function would ever be complex when p is restricted to be real, since all of the coefficients are real.
Also, how exactly does it restrict function values to be real while it's solving for the roots? And why are the roots different from the case where complex function values are allowed?
Sorry that this is kind of a mess, but I would really appreciate any help you can provide!
upperbound=1/2will make this an exact input and the different methods will deliver the same result (around 0.333333333613, when evaluated numerically). $\endgroup$