Mathematica：多项式实根的分支点

发布于 2024-10-07 14:58:29 字数 1916 浏览 1 评论 0原文

我正在对以下示例函数进行“梯度极值”的强力搜索，

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

，Reduce 对我来说很高兴：

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

geyvals 是三次多项式，但表达式放在这里有点大。

现在我的问题：对于 x 的不同值，这些根的不同数量是真实的，我想挑选出解决方案按顺序分支的 x 的值将沿谷底（fv）的梯度极值拼凑在一起。在目前的情况下，由于多项式只是三次，我可能可以手动完成它 - 但我正在寻找一种让 Mathematica 为我完成它的简单方法？

编辑：澄清一下：梯度极值只是背景——也是设置难题的简单方法。我对这个问题的具体解决方案并不感兴趣，而是对发现多项式根的分支点的一般传递方式感兴趣。在下面添加了一个带有工作方法的答案。

编辑2：因为看起来实际问题比根分支有趣得多：rcollyer建议直接在gecond上使用ContourPlot来获取梯度极值。为了完成这个任务，我们需要分离山谷和山脊，这是通过查看垂直于梯度的 Hessian 矩阵的特征值来完成的。将“valleynes”检查作为 RegionFunction 我们只剩下谷线：

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

它只给出谷底线。包括一些轮廓和鞍点：

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

我们最终得到这样的图：

Contourplot of fv with the Valley line superpped

原文

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[List@ToRules@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}, Evaluate@valleycond], 
   PlotPoints -> 41];

Which gives just the valley floor line. Including some contours and the 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"}]

We end up with a plot like this:

Contour plot of fv with the valley line superposed

分享到QQ

分享到微博

如果你对这篇内容有疑问，欢迎到本站社区发帖提问参与讨论，获取更多帮助，或者扫码二维码加入 Web 技术交流群。

发布评论

需要登录才能够评论，你可以免费注册一个本站的账号。

活泼老夫 2024-10-14 14:58:29

不确定这（迟来的）是否有帮助，但似乎您对判别点感兴趣，即多项式和导数（wrt 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}}

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

回复收藏 0 原文

影子是时光的心 2024-10-14 14:58:29

如果您只想绘制结果，请使用 StreamPlot[]< /a> 渐变：

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。

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.

回复收藏 0 原文

远昼 2024-10-14 14:58:29

更新：见下文。

我首先通过可视化根的虚部来解决这个问题：

plot of the imaginary部分的根

这告诉你立即发生三件事：1）第一个根始终是实数，2）后两个是共轭对，3）有一个接近零的小区域，其中所有三个都是实数。另外，请注意，排除仅消除了 x=0 处的奇点，当我们放大时我们可以看到原因：

zoom in of 上述照片，x 在 0 到 1.5 之间

然后我们可以使用 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< /code> 规范有点奇怪，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 + dx + ax^2 + ey + 2 cxy + by^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];

得到涉及 x 和 p 的 q 的多个公式。对于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}]}]

img src="https://i.sstatic.net/N5eo7.png" alt="xy 图显示梯度极值和有 3 个实根的区域">

这证实了前两个图中显示 3 个实根的区域。

Updated: see below.

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

plot of 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:

zoom in of above photo with x between 0 and 1.5

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[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, 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}]},
Show@{
 RegionPlot[Evaluate[Thread[Re[qvals]==0]/.p-> y], rngs],
 ContourPlot[g.RotationMatrix[Pi/2].h.g==0,rngs 
      ContourStyle -> {Darker@Red,Dashed}]}]

gives

x-y plot showing gradient extremals and region where there are 3 real roots

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

回复收藏 0 原文

幸福不弃 2024-10-14 14:58:29

最终我自己尝试了，因为目标确实是“放手”。我将把这个问题暂时搁置一段时间，看看是否有人找到更好的方法。

下面的代码使用二分法将 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]
    ]]]

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

回复收藏 0 原文

~没有更多了~