For a polynomial F with dependence on a parameter beta, find the value of beta for which F has 2 positive roots

Question

Good morning,

I have the following polynomial with degree 4 and its coefficients are quadratic functions of a parameter beta:

FF[z_, β_] := z^4 + z^3 (4.67687 - 0.16 β^2) + z^2 (-0.703991 - 0.652433 β^2) +
  z (0.0676982 - 0.049884 β^2) - 0.000608244 β^2 + 0.0000234502;

Consider z > 0 and β > 0. From the code

Manipulate[{Plot[FF[z, param], {z, 0, 0.2}, PlotRange -> All, 
  PlotLabel -> FF[z, param]], 
  Column[z /. NSolve[FF[z, param] == 0, z]], CountRoots[FF[z, param], {z, 0, ∞}]},
  {{param, 0.6}, 0.18, 0.8}]

I know that for some $\beta$ (close to 0.6), the quartic polynomial has exactly 2 positive roots. I'd like to find such a value of beta.

I've tried several things. For instance:

FindRoot[NMinimize[{FF[z, β], 0.05 < z < 0.15}, z], {β, 0.6}]

and

NSolve[CountRoots[FF[z, β], {z, 0, Infinity}] == 2, {β, 0.6, 0.8}]

I appreciate some suggestions.

Answer 1

You can brute force this using a method inspired by binary search.

Your function:

FF[z_, b_] := 
  z^4 + z^3 (4.67687 - 0.16 b^2) + 
   z^2 (-0.703991 - 0.652433 b^2) + 
   z (0.0676982 - 0.049884 b^2) - 0.000608244 b^2 + 
   0.0000234502;

Function which cuts the min/max range in half each time, outputting the new min/max to search:

f[{x1_, x2_}] := 
 Module[{v1, v2, v3, min = N[x1, 30], max = N[x2, 30], vals},
  v1 = NSolve[FF[z, min] == 0, {z}, Reals];
  v2 = NSolve[FF[z, max] == 0, {z}, Reals];
  v3 = NSolve[FF[z, (min + max)/2] == 0, {z}, Reals];

  If[Length@v3 <= 3,
   vals = {(min + max)/2, max},
   vals = {min, (min + max)/2}
   ];
  vals
  ]

And using our custom function above with FixedPoint we get a range that is hopefully specific enough for your application (with the floating point numbers in your original equations we're restricted to $MachinePrecision):

In[9]:= FixedPoint[f, {.6, .61}]
Out[9]= {0.6070523316794493, 0.6070523316794494}

So as @Mariusz Iwaniuk said in the comments, you're looking for a number in the ballpark of $0.607052331679449$.

Note that my implementation isn't very general and relies on your specific FF being the format it is. No doubt it could be generalized to fit a broader scope of function definitions, though.

There may be a more elegant math-y way to solve it, but I just resort to computational/numerical methods above.

Answer 2

There is no region in the positive parameter space for which there are exactly two positive roots. We can show this as follows. First note that the number of positive roots can only change by a root crossing the axis or by two roots joining (a multiple roots) and becoming complex valued as the parameter changes. The first case is accounted for by solving for the parameter when z is zero. The second arises at roots of the discriminant polynomial. All we need to do is find all these relative change points, and count positive roots at the centers of the intervals between them (since the number of positive roots is constant in these intervals).

fF[z_, \[Beta]_] := 
  z^4 + z^3 (4.67687 - 0.16 \[Beta]^2) + 
   z^2 (-0.703991 - 0.652433 \[Beta]^2) + 
   z (0.0676982 - 0.049884 \[Beta]^2) - 0.000608244 \[Beta]^2 + 
   0.0000234502;

discrim = Discriminant[fF[z, b], z]

(* Out[27]= -0.0914115 + 0.386364 b^2 - 0.353399 b^4 - 0.062053 b^6 + 
 0.00686865 b^8 + 0.000346896 b^10 + 9.23814*10^-6 b^12 *)

changePts1 = 
 Select[b /. NSolve[discrim], FreeQ[#, Complex] && # > 0 &]
changePts2 = 
 Select[b /. NSolve[fF[z, b] /. z -> 0], FreeQ[#, Complex] && # > 0 &]
changePts = Sort[Join[changePts1, changePts2]]

(* Out[42]= {0.607052, 0.779012, 3.00021}

Out[43]= {0.196352}

Out[44]= {0.196352, 0.607052, 0.779012, 3.00021} *)

extendedPts = Join[{0}, changePts, {Last[changePts] + 1}];
testPts = Map[Mean, Partition[extendedPts, 2, 1]]

(* Out[46]= {0.0981758, 0.401702, 0.693032, 1.88961, 3.50021} *)

solns = Map[z /. NSolve[fF[z, #]] &, testPts]

(* Out[47]= {{-4.82541, -0.000260933, 0.0751716 - 0.0912008 I, 
  0.0751716 + 0.0912008 I}, {-4.82147, 0.00127418, 
  0.0845699 - 0.070761 I, 0.0845699 + 0.070761 I}, {-4.81328, 
  0.00736127, 0.0480394, 
  0.157853}, {-4.74059, -0.0177988 - 0.0189466 I, -0.0177988 + 
   0.0189466 I, 0.670623}, {-4.58696, -0.0412625, -0.0203156, 
  1.93191}} *)

We see that no interval has exactly two positive roots (between 0.607052 and 0.779012 there are three positive roots). If you ignore multiplicity, there are exactly two at the endpoints of that particular interval.

Map[z /. NSolve[fF[z, #]] &, changePts]

(* Out[48]= {{-4.82466, 0., 0.0769773 - 0.0877929 I, 
  0.0769773 + 0.0877929 I}, {-4.81613, 0.00443875, 0.0968922, 
  0.0968922}, {-4.81006, 0.0193691, 0.0193691, 
  0.191554}, {-4.63723, -0.0284017, -0.0284017, 1.45737}} *)

Answer 3

Let's try a more algebraic solution! First we Rationalize the polynomial so we can work symbolically until the end and be sure we get maximum precision. In case you have symbolic expressions for your polynomial coefficients you can insert those instead and solve the problem fully analytically right until you want the decimal representation of the solution.

poly = Numerator[Rationalize[
z^4 + z^3 (4.67687 - 0.16 \[Beta]^2) + 
 z^2 (-0.703991 - 0.652433 \[Beta]^2) + 
 z (0.0676982 - 0.049884 \[Beta]^2) - 0.000608244 \[Beta]^2 + 
 0.0000234502, 0] // Simplify]

117251 + 5000000000 z^4 - 3041220 \[Beta]^2 - 50000 z^3 (-467687 + 16000 \[Beta]^2) - 1000 z (-338491 + 249420 \[Beta]^2) - 5000 z^2 (703991 + 652433 \[Beta]^2)

Here we cleaned up a bit by throwing away the denominator because it doesn't affect the roots.

To get a crude idea how the number of zeros changes depending on $\beta$ let's plot an approximation:

zeros = z /. Solve[poly == 0, z];
Plot[
  Total[Boole[Im[#] == 0] & /@ zeros]
  ,{\[Beta], 0, 1}, WorkingPrecision -> 30, PlotRange -> {0, 4}
]

Here we can see we are interested in the change around 0.6 when the number of real roots changes from two to four. We shouldn't use this form directly though since the Im[#]==0 part is numerically unstable and therefore not reliable for a solution.

Now, as J.M. hinted to in the comments we can use the Discriminant of the polynomial to get that exact $\beta$ value where the number of real roots changes. Let's look at a plot first:

Plot[Discriminant[poly, z], {\[Beta], 0, 1}]

which looks exactly like what we want. So now the only thing left to do is to get the root of the discriminant close to 0.6 and we are done!

\[Beta] /. First@Solve[
  Discriminant[expr, z] == 0 \[And] 1/2 < \[Beta] < 3/4, \[Beta]
]
N[%, 20]

Root[-228528823144524583756372813176499152547 + 965910347668582228294215219952243710680 #1^2 - 883498027761558572591308192199443792200 #1^4 - 155132407672471157729945150782311491400 #1^6 + 17171627184838703747107136205337068000 #1^8 + 867239413153808414882278532192000000 #1^10 + 23095356254911115872300032000000000 #1^12 &, 4]

0.60705233167945003248

Answer 4

FF[z_, β_] := Rationalize[z^4 + z^3 (4.67687 - 0.16 β^2) + z^2 (-0.703991 - 0.652433 β^2) + 
z (0.0676982 - 0.049884 β^2) - 0.000608244 β^2 + 0.0000234502, 0];

f[β_?NumericQ] := (NMinimize[FF[z, β], 5/100 < z < 15/100, z, WorkingPrecision -> 20, Method -> "RandomSearch"][[1]]); 
FindRoot[f[β] == 0, {β, 6/10}, WorkingPrecision -> 20]

(* {β -> 0.60705233167944997062} *)

Answer 5

First an extended comment:.

You could use ContourPlot to get an overview.

  eq=z^4 + z^3 (4.67687 - 0.16 β^2) +z^2 (-0.703991 - 0.652433 β^2) +z (0.0676982 - 0.049884 β^2) -0.000608244 β^2 +0.0000234502
  ContourPlot[0 ==eq, {β, 0, 1}, {z, -6, 1}, MaxRecursion -> 5 ,FrameLabel -> {β, z}]

The solutionplot shows two ranges β < 0.6 and β > 0.78 with only two solutions!

solution

The intersting branch of solution is in the range z0 < z < z1

{β0, z0} = {0, z /. Solve[eq == 0 /. β -> 0, z][[2]]} 
(*{0, -0.000345152}*)
{β1, z1} = {1, z /. Solve[eq == 0 /. β -> 1, z][[4]]}
(*{1, 0.271234}*)

Now it's easy to calculate

beta = β /. Solve[eq == 0, β][[2]] (* beta[z]*)

ParametricPlot[{z, beta}, {z, z0, z1}, AxesLabel -> {"z", "β"}]

The two extrema define the β-values we are looking for

beta /. NSolve[{D[beta, z] == 0, 0 < z < z1}, z]
(*{0.779012, 0.607052}*)

That's it. Hope it helps!