Mathematics: Branch Points for Real Roots of a Polynomial

Question

I search for brute force for “gradient extremals” in the following example function

fv[{x_, y_}] = ((y - (x/4)^2)^2 + 1/(4 (1 + (x - 1)^2)))/2; 

This includes finding the next zeros

 gecond = With[{g = D[fv[{x, y}], {{x, y}}], h = D[fv[{x, y}], {{x, y}, 2}]}, g.RotationMatrix[Pi/2].hg == 0] 

Which Reduce happily does for me:

 geyvals = y /. Cases[ List@ToRules @Reduce[gecond, {x, y}], {y -> _}]; 

geyvals are the three roots of the cubic polynomial, but the expression here is a little large.

Now to my question: for different values ​​of x different numbers of these roots are valid, and I would like to highlight the values ​​of x , where are the decision branches in order to compose the gradient extremals along the valley floor ( fv ). In this case, since the polynomial is only cubic, I could do it manually - but am I looking for an easy way to get Mathematica to do this for me?

Change To clarify: the material of gradient extremals is just a background - and an easy way to fix a complex problem. I am not interested in a specific solution to this problem, as in the general method for determining branch points for polynomial roots. Added the answer below with a working approach.

Change 2 . Since it seems that the actual problem is much more interesting than root branching: rcollyer suggests using ContourPlot directly on gecond to get extreme gradients. To do this completely, we need to separate the valleys and ridges, which is done by looking at the eigenvalue of the Hessian perpendicular to the gradient. Putting a check on "valleynes" as a RegionFunction , we only stay with the valley line:

 valleycond = With[{ g = D[fv[{x, y}], {{x, y}}], h = D[fv[{x, y}], {{x, y}, 2}]}, g.RotationMatrix[Pi/2].h.RotationMatrix[-Pi/2].g >= 0]; gbuf["gevalley"]=ContourPlot[gecond // Evaluate, {x, -2, 4}, {y, -.5, 1.2}, RegionFunction -> Function[{x, y}, Evaluate@valleycond ], PlotPoints -> 41]; 

This gives only the valley floor line. Including some contours and a saddle point:

 fvSaddlept = {x, y} /. First@Solve [Thread[D[fv[{x, y}], {{x, y}}] == {0, 0}]] gbuf["contours"] = ContourPlot[fv[{x, y}], {x, -2, 4}, {y, -.7, 1.5}, PlotRange -> {0, 1/2}, Contours -> fv@fvSaddlept (Range[6]/3 - .01), PlotPoints -> 41, AspectRatio -> Automatic, ContourShading -> None]; gbuf["saddle"] = Graphics[{Red, Point[fvSaddlept]}]; Show[gbuf /@ {"contours", "saddle", "gevalley"}] 

The result is the following plot:

Answer 1

I'm not sure if this (belatedly) helps, but it seems to you that you are interested in discriminant points, that is, both the polynomial and the derivative (wrt y) vanish. You can solve this system for {x, y} and throw away complex solutions, as shown below.

 fv[{x_, y_}] = ((y - (x/4)^2)^2 + 1/(4 (1 + (x - 1)^2)))/2; gecond = With[{g = D[fv[{x, y}], {{x, y}}], h = D[fv[{x, y}], {{x, y}, 2}]}, g.RotationMatrix[Pi/2].hg] In[14]:= Cases[{x, y} /. NSolve[{gecond, D[gecond, y]} == 0, {x, y}], {_Real, _Real}] Out[14]= {{-0.0158768, -15.2464}, {1.05635, -0.963629}, {1., 0.0625}, {1., 0.0625}} 

Answer 2

If you only want to build the result, use StreamPlot[] on the gradients:

 grad = D[fv[{x, y}], {{x, y}}]; StreamPlot[grad, {x, -5, 5}, {y, -5, 5}, RegionFunction -> Function[{x, y}, fv[{x, y}] < 1], StreamScale -> 1] 

You may have to play with plot accuracy, StreamStyle and RegionFunction to achieve excellence. It would be especially useful to use the Valley Floor solution for StreamPoints software programming.

Answer 3

Updated : see below.

I would approach this first by visualizing the imaginary parts of the roots:

This immediately tells you three things: 1) the first root is always real; 2) the second - conjugated pairs; 3) there is a small region near zero in which all three are real. Also, note that the exceptions only got rid of the singular point at x=0 , and we can understand why when we zoom in:

Then we can use EvalutionMonitor to generate a list of roots directly:

 Map[Module[{f, fcn = #1}, f[x_] := Im[fcn]; Reap[Plot[f[x], {x, 0, 1.5}, Exclusions -> {True, f[x] == 1, f[x] == -1}, EvaluationMonitor :> Sow[{x, f[x]}][[2, 1]] // SortBy[#, First] &];] ]&, geyvals] 

(Note: the Part specification is a bit odd, Reap returns a List of what is sown as the second element in the List , so this leads to a nested list. In addition, Plot does not try points in a simple way, so SortBy is required.) There may be a more elegant way. to determine where the last two roots become complex, but since their imaginary parts are piecewise continuous, it just seemed easier for brute force.

Change Since you mentioned that you want to create an automatic generation method when some of the roots become complex, I studied what happens when you replace y -> p + I q . Now this assumes x is real, but you already did it in your solution. In particular, I do the following

 In[1] := poly = g.RotationMatrix[Pi/2].hg /. {y -> p + I q} // ComplexExpand; In[2] := {pr,pi} = poly /. Complex[a_, b_] :> a + zb & // CoefficientList[#, z] & // Simplify[#, {x, p, q} \[Element] Reals]&; 

where the second step allows me to isolate the real and imaginary parts of the equation and simplify them independently of each other. Performing the same task with the common two-dimensional polynomial f + dx + ax^2 + ey + 2 cxy + by^2 , but creating both the complex x and y ; I noted that Im[poly] = Im[x] D[poly, Im[x]] + Im[y] D[poly,[y]] , and this may be the case for your equation. Having made x real, the imaginary part of poly becomes q times some function of x , p and q . Thus, setting q=0 always gives Im[poly] == 0 . But this does not tell us anything new. However, if we

 In[3] := qvals = Cases[ List@ToRules @RReduce[ pi == 0 && q != 0, {x,p,q}], {q -> a_}:> a]; 

we get several formulas for q involving x and p . For some values ​​of x and p these formulas can be imaginary, and we can use Reduce to determine where Re[qvals] == 0 . In other words, we want the "imaginary" part of y be real, and this can be done by allowing q be zero or purely imaginary. Construction of the region where Re[q]==0 and the imposition of gradient extremal lines through

 With[{rngs = Sequence[{x,-2,2},{y,-10,10}]}, Show@ { RegionPlot[Evaluate[Thread[Re[qvals]==0]/.p-> y], rngs], ContourPlot[g.RotationMatrix[Pi/2].hg==0,rngs ContourStyle -> { Darker@Red ,Dashed}]}] 

gives

which confirms the areas in the first two graphs showing 3 real roots.

Answer 4

Ended up, tried myself, since the goal was really to do it "hands." I will leave the question open for a while to find out if anyone will find a better way.

The code below uses bisection to snap points, where CountRoots changes the value. This works for my case (defining a singularity at x = 0 is pure luck):

 In[214]:= findRootBranches[Function[x, Evaluate@geyvals [[1, 1]]], {-5, 5}] Out[214]= {{{-5., -0.0158768}, 1}, {{-0.0158768, -5.96046*10^-9}, 3}, {{0., 0.}, 2}, {{5.96046*10^-9, 1.05635}, 3}, {{1.05635, 5.}, 1}} 

Implementation:

 Options[findRootBranches] = { AccuracyGoal -> $MachinePrecision/2, "SamplePoints" -> 100}; findRootBranches::usage = "findRootBranches[f,{x0,x1}]: Find the the points in [x0,x1] \ where the number of real roots of a polynomial changes. Returns list of {<interval>,<root count>} pairs. f: Real -> Polynomial as pure function, eg f=Function[x,#^2-x&]." ; findRootBranches[f_, {xa_, xb_}, OptionsPattern[]] := Module[ {bisect, y, rootCount, acc = 10^-OptionValue[AccuracyGoal]}, rootCount[x_] := {x, CountRoots[f[x][y], y]}; (* Define a ecursive bisector w/ automatic subdivision *) bisect[{{x1_, n1_}, {x2_, n2_}} /; Abs[x1 - x2] > acc] := Module[{x3, n3}, {x3, n3} = rootCount[(x1 + x2)/2]; Which[ n1 == n3, bisect[{{x3, n3}, {x2, n2}}], n2 == n3, bisect[{{x1, n1}, {x3, n3}}], True, {bisect[{{x1, n1}, {x3, n3}}], bisect[{{x3, n3}, {x2, n2}}]}]]; (* Find initial brackets and bisect *) Module[{xn, samplepoints, brackets}, samplepoints = N@With [{sp = OptionValue["SamplePoints"]}, If[NumberQ[sp], xa + (xb - xa) Range[0, sp]/sp, Union[{xa, xb}, sp]]]; (* Start by counting roots at initial sample points *) xn = rootCount /@ samplepoints; (* Then, identify and refine the brackets *) brackets = Flatten[bisect /@ Cases[Partition[xn, 2, 1], {{_, a_}, {_, b_}} /; a != b]]; (* Reinclude the endpoints and partition into same-rootcount segments: *) With[{allpts = Join[{ First@xn }, Flatten[brackets /. bisect -> List, 2], { Last@xn }]}, {#1, Last[#2]} & @@@ Transpose /@ Partition[allpts, 2] ]]]