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.

  • With brute force search: $\beta =0.60705233167945$ – Mariusz Iwaniuk Oct 20 at 22:11
  • @user6014 Well, I'm interested about the positive roots. $\beta = 0.6$ and $\beta=0.601$ provide only one positive solution ;) – Alex Pereira Oct 20 at 23:46
  • @mariusz-iwaniuk, can you show us how to do it? I'm asking this because I'll probably have to consider this problem again but with different numbers on the coefficients of F – Alex Pereira Oct 20 at 23:51
  • 2
    This article might be useful. – J. M. is computer-less Oct 21 at 9:33

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.

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}} *)

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

number of real roots

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

discriminant plot

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

  • Again, thanks for following through! :) – J. M. is computer-less Oct 21 at 14:59
  • @J.M. You always seem to have the right intuition for a solution ;) First I wanted to write a solution based on the Total[Boole[Im[#]==0]&/@zeros]==4 idea together with existential quantifiers and Resolve but the Im/@zeros expressions are way too messy, so the discriminant was the right way to go! :) – Thies Heidecke Oct 21 at 15:05
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} *)
  • Thank you all guys. I've learned a lot with your suggestions and comments. I'm pretty satisfied with the possibilities... (Y) – Alex Pereira Oct 21 at 20:17

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

plot

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

plot

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

That's it. Hope it helps!

Your Answer

 
discard

By clicking "Post Your Answer", you acknowledge that you have read our updated terms of service, privacy policy and cookie policy, and that your continued use of the website is subject to these policies.

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