Mathematica: точки ветвления действительных корней многочлена

Не уверен, что это (с опозданием) помогает, но, похоже, вас интересуют дискриминантные точки, то есть, где и полиномиальные, и производные (относительно y) исчезают. Вы можете решить эту систему для {x,y} и отбросить сложные решения, как показано ниже.

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

Если вы хотите только отобразить результат, используйте StreamPlot[] для градиентов:

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]

Возможно, вам придется повозиться с точностью графика, StreamStyle и RegionFunction, чтобы добиться идеального результата. Особенно полезным было бы использование решения для дна долины для программного заполнения StreamPoints.

Обновлено: см. ниже.

Я бы подошел к этому сначала, визуализируя воображаемые части корней:

Это сразу говорит вам о трех вещах: 1) первый корень всегда действителен, 2) вторые два являются сопряженными парами и 3) существует небольшая область около нуля, в которой все три действительны. Кроме того, обратите внимание, что исключения избавились только от особой точки в x=0, и мы можем понять, почему, когда увеличим масштаб:

Затем мы можем использовать EvalutionMonitor для генерирования списка корней напрямую:

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]

(Обратите внимание, что спецификация Part немного странная, Reap возвращает List того, что посеяно как второй элемент в List, так что получается вложенный список. Кроме того, Plot не выбирает точки прямым способом, поэтому SortBy.) Может быть более элегантный способ определить, где два последних корня становятся комплексными, но, поскольку их мнимые части кусочно-непрерывны, кажется, что проще выполнить его методом грубой силы.

Правка: поскольку вы упомянули, что вам нужен автоматический метод для генерации, когда некоторые корни становятся сложными, я исследовал, что происходит, когда вы подставляете y -> p + I q. Теперь это предполагает, что x реально, но вы уже сделали это в своем решении. В частности, я делаю следующее

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]&;

где второй шаг позволяет мне изолировать действительную и мнимую части уравнения и упростить их независимо друг от друга. Делаем то же самое с общим двумерным полиномом f + d x + a x^2 + e y + 2 c x y + b y^2, но делаем как x, так и y комплексными; Я отметил, что Im[poly] = Im[x] D[poly, Im[x]] + Im[y] D[poly,[y]], и это может относиться и к вашему уравнению. Делая x реальным, мнимая часть poly становится q раз некоторой функцией x, p и q. Итак, установка q=0 всегда дает Im[poly] == 0. Но это не говорит нам ничего нового. Однако, если мы

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

мы получаем несколько формул для q с участием x и p. Для некоторых значений x и p эти формулы могут быть воображаемыми, и мы можем использовать Reduce, чтобы определить, где Re[qvals] == 0. Другими словами, мы хотим, чтобы «мнимая» часть y была реальной, а этого можно добиться, разрешив q быть равным нулю или чисто мнимым. Построение области, где Re[q]==0, и наложение экстремальных линий градиента через

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

дает

что подтверждает области на первых двух графиках, показывающие 3 действительных корня.

Закончилось тем, что я попробовал себя, так как целью действительно было сделать это «без рук». Я оставлю вопрос открытым на некоторое время, чтобы посмотреть, найдет ли кто-нибудь лучший способ.

В приведенном ниже коде используется деление пополам для заключения в скобки точек, в которых CountRoots изменяет значение. Это работает для моего случая (обнаружение сингулярности при x = 0 - чистая удача):

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

Выполнение:

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