# Mathematica: branch points for real roots of polynomial

I am doing a brute force search for "gradient extremals" on the following example function

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

This involves finding the following zeros

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

Which `Reduce` happily does for me:

`geyvals = y /. Cases[[email protected]@Reduce[gecond, {x, y}], {y -> _}]; `

`geyvals` is the three roots of a cubic polynomial, but the expression is a bit large to put here.

Now to my question: For different values of `x`, different numbers of these roots are real, and I would like to pick out the values of `x` where the solutions branch in order to piece together the gradient extremals along the valley floor (of `fv`). In the present case, since the polynomial is only cubic, I could probably do it by hand -- but I am looking for a simple way of having Mathematica do it for me?

Edit: To clarify: The gradient extremals stuff is just background -- and a simple way to set up a hard problem. I am not so interested in the specific solution to this problem as in a general hand-off way of spotting the branch points for polynomial roots. Have added an answer below with a working approach.

Edit 2: Since it seems that the actual problem is much more fun than root branching: rcollyer suggests using `ContourPlot` directly on `gecond` to get the gradient extremals. To make this complete we need to separate valleys and ridges, which is done by looking at the eigenvalue of the Hessian perpendicular to the gradient. Putting a check for "valleynes" in as a `RegionFunction` we are left with only 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}, [email protected]], PlotPoints -> 41]; `

Which gives just the valley floor line. Including some contours and the saddle point:

`fvSaddlept = {x, y} /. [email protected][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 -> [email protected] (Range[6]/3 - .01), PlotPoints -> 41, AspectRatio -> Automatic, ContourShading -> None]; gbuf["saddle"] = Graphics[{Red, Point[fvSaddlept]}]; Show[gbuf /@ {"contours", "saddle", "gevalley"}] `

We end up with a plot like this:

Not sure if this (belatedly) helps, but it seems you are interested in discriminant points, that is, where both polynomial and derivative (wrt y) vanish. You can solve this system for {x,y} and throw away complex solutions as 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].h.g] 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}} ```

If you only want to plot the result then 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 fiddle around with the plot's precision, StreamStyle, and the RegionFunction to get it perfect. Especially useful would be using the solution for the valley floor to seed `StreamPoints` programmatically.

Updated: see below.

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

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

We can then use the `EvalutionMonitor` to generate the 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 little odd, `Reap` returns a `List` of what is sown as the second item in a `List`, so this results in a nested list. Also, `Plot` doesn't sample the points in a straightforward manner, so `SortBy` is needed.) There may be a more elegant route to determine where the last two roots become complex, but since their imaginary parts are piecewise continuous, it just seemed easier to brute force it.

Edit: Since you've mentioned that you want an automatic method for generating where some of the roots become complex, I've been exploring what happens when you substitute in `y -> p + I q`. Now this assumes that `x` is real, but you've already done that in your solution. Specifically, I do the following

```In[1] := poly = g.RotationMatrix[Pi/2].h.g /. {y -> p + I q} // ComplexExpand; In[2] := {pr,pi} = poly /. Complex[a_, b_] :> a + z b & // 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 independent of each other. Doing this same thing with the generic 2D polynomial, `f + d x + a x^2 + e y + 2 c x y + b y^2`, but making both `x` and `y` complex; I noted that `Im[poly] = Im[x] D[poly, Im[x]] + Im[y] D[poly,[y]]`, and this may hold for your equation, also. By making `x` real, the imaginary part of `poly` becomes `q` times some function of `x`, `p`, and `q`. So, setting `q=0` always gives `Im[poly] == 0`. But, that does not tell us anything new. However, if we

```In[3] := qvals = Cases[[email protected]@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`, those formulas may be imaginary, and we can use `Reduce` to determine where `Re[qvals] == 0`. In other words, we want the "imaginary" part of `y` to be real and this can be accomplished by allowing `q` to be zero or purely imaginary. Plotting the region where `Re[q]==0` and overlaying the gradient extremal lines via

```With[{rngs = Sequence[{x,-2,2},{y,-10,10}]}, [email protected]{ RegionPlot[Evaluate[Thread[Re[qvals]==0]/.p-> y], rngs], ContourPlot[g.RotationMatrix[Pi/2].h.g==0,rngs ContourStyle -> {[email protected],Dashed}]}] ```

gives

which confirms the regions in the first two plots showing the 3 real roots.

Ended up trying myself since the goal really was to do it 'hands off'. I'll leave the question open for a good while to see if anybody finds a better way.

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

```In[214]:= findRootBranches[Function[x, [email protected][[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, e.g 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 = [email protected][{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[{[email protected]}, Flatten[brackets /. bisect -> List, 2], {[email protected]}]}, {#1, Last[#2]} & @@@ Transpose /@ Partition[allpts, 2] ]]] ```

Category:wolfram mathematica Views:0 Time:2010-12-14

## Related post

• SVN - Branched tag back to root 2010-01-04

I have a svn repository. I have now created a branch for this repository to do some major dev work. Now some of the features that I can created in my project can be committed directly to the main repository instead of the branch. Is there a way to do

• How to find the roots of polynomial equation of 'n'th degree in Java? 2011-05-07

Help me solving the above problem. --------------Solutions------------- as I said in the comment it is proven that there is no algorithm for n>=5. however, you can find a numerical approximation (i.e. Newton`s Method) If n >= 5, then according

• Is there a special name (like HEAD, FETCH_HEAD) for the root commit on the current branch in a git repository? 2010-10-01

I find myself referring to the root commit for various reasons. I usually make a tag called first or root to make it easier, but I have been dealing with quite a few repositories lately and even this manual tagging is getting quite tedious. Is there

• TFS 2008 Branching always adds root folder 2010-12-15

I'm trying to create a branch for QA using TFS 2008 I have my main code in the trunk folder. When I perform the branch, I select a qa folder which I created to be mirror the trunk folder. However, when I perform the branch, it puts the trunk folder i

• What is a simple way to find real roots of a (cubic) polynomial? 2011-02-05

this seems like an obvious question to me, but I couldn't find it anywhere on SO. I have a cubic polynomial and I need to find real roots of the function. What is THE way of doing this? I have found several closed form formulas for roots of a cubic f

• Traverse a branch top to bottom. [root and destination node given] 2011-06-09

How to traverse a branch of n-ary tree if you just have a source node and a destination node. Rest everything can be assumed. Given that source node is ancestor of that destination node. --------------Solutions------------- yes dfs is one way we can

• File permissions resetting when switching branches in Git 2009-12-09

I created a Git repository on a folder that had a different Linux owner than my user. It wasn't until much later that I set the group permissions to write so that my user could make changes and commits to Git. However, whenever I switch from a branch

• SVN naming convention: repository, branches, tags 2010-05-26

Just curious what your naming conventions are for the following: Repository name Branches Tags Right now, we're employing the following standards with SVN, but would like to improve on it: Each project has its own repository Each repository has a set

• Labelling mercurial revision from teamcity fails with "push creates new remote heads on branch 'default'" 2010-09-17

I have a build set up in Teamcity that builds and tests a mercurial branch, and should then tag that branch. The building works correctly but when it comes to labelling it fails with the error "push creates new remote heads on branch 'default'". I fi

• How do I use git-tfs and idiomatic git branching against a TFS repository? 2011-02-26

How Do I Use git-tfs Idiomatically? The git idiom is to check out branches to the root directory of the repository. Checking out a branch will replace the contents of the directory with the contents of that branch. The TFS idiom is to check out each

• Howto create a branch from tag in CVS with intelliJ 2011-02-28

I checked a specific tag representing some release from CVS using intelliJ Idea 9.0. Now I made few fixes to the code and want to create a branch for the same release starting from this tag. So what I want to do is exactly this: Start branch (before

• Git svn: how do I work with "partial" branches? 2011-05-23

By "partial", I mean the branches in my SVN repository were not copied from trunk/, but from a subdirectory of trunk/. Here's the svn repo layout: trunk/ core projectA projectB ... branches/ core_new_feature (branched from "core" on trunk) another_br

• git: How to split off library from project? filter-branch, subtree? 2011-06-19

So, I've a bigger (closed source) project, and in the context of this project created a library which could also be useful elsewhere, I think. I now want to split off the library in its own project, which could go as open source on github or similar.

• What is the branch in the destructor reported by gcov? 2011-08-26

When I use gcov to measure test coverage of C++ code it reports branches in destructors. struct Foo { virtual ~Foo() { } }; int main (int argc, char* argv[]) { Foo f; } When I run gcov with branch probabilities enabled (-b) I get the following output

• How can I commit specific directory to a different than a working branch on git? 2011-08-27

I have a simple static website, based on HTML 5 boilerplate and the whole thing is in a github repository. The project structure is something like this build css img js publish # this directory isn't in source control! ... index.html & other file

• How do I find the minimum/maximum values to use with Root Solvers? 2011-10-03

I want to use the root solvers (ex: BrentSolver) in Commons Math to find roots for polynomial functions, but they all seem to require using an initial estimate for minimum/maximum, where the function has different signals. So how do I go about doing

• Root finding with quadrature 2011-10-16

I have a problem with root finding and am having difficulty in getting it to work in this instance. Some complicated function I need. f[x_, lambda_, alpha_, beta_, mu_] = Module[{gamma}, gamma = Sqrt[alpha^2 - beta^2]; (gamma^(2*lambda)/((2*alpha)^(l

• Mathematica; NDSolve; Can I use InterpolatingFunction as initial conditions? 2011-10-28

The questions arises because I want to use the solution of one PDE as initial condition to solve another. Since in NDSolve, the solution is given by InterpolatingFunction, I have to use InterpolatingFunction in the second PDE. Is this possible? Why i

• In TFS is the direction of a branch's parent-child relationship important 2012-04-25

Is the direction of the parent-child relationship between branches important? Is parent/child just an abstract concept that identifies which was the source and which was the destination, or are there specific operations that can only be performed in