Mathematica: điểm chi nhánh cho rễ thực sự của đa thức

Question

Tôi đang làm một tìm kiếm brute force cho "cực trị gradient" trên ví dụ hàm sau

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

này liên quan đến việc tìm kiếm các số không sau

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

nào Reduce vui vẻ cho tôi:

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

geyvals là ba gốc của một đa thức khối, nhưng biểu thức là một chút lớn để đặt ở đây.

Bây giờ câu hỏi của tôi: Đối với các giá trị khác nhau của x, các số khác nhau của những gốc này là có thật và tôi muốn chọn các giá trị của x nơi chi nhánh giải pháp để ghép lại với nhau. (trong số fv). Trong trường hợp hiện tại, vì đa thức chỉ là khối, tôi có thể làm điều đó bằng tay - nhưng tôi đang tìm kiếm một cách đơn giản để Mathematica làm điều đó cho tôi?

Chỉnh sửa: Để làm rõ: Công cụ Gradient extremals chỉ là nền - và cách đơn giản để thiết lập sự cố khó. Tôi không quan tâm đến giải pháp cụ thể cho vấn đề này như trong một cách thức tổng quát để phát hiện các điểm chi nhánh cho các rễ đa thức. Đã thêm câu trả lời dưới đây với phương pháp làm việc.

Chỉnh sửa 2: Kể từ khi có vẻ như vấn đề thực sự là niềm vui nhiều hơn phân nhánh rễ: rcollyer gợi ý sử dụng ContourPlot trực tiếp trên gecond để có được cực trị gradient. Để thực hiện điều này, chúng ta cần tách các thung lũng và rặng núi, được thực hiện bằng cách nhìn vào giá trị riêng của Hessian vuông góc với gradient. Đưa một tấm séc cho "valleynes" trong như một RegionFunction chúng tôi là trái với chỉ dòng thung lũng:

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

Mà cho chỉ dòng đáy thung lũng. Trong đó có một số đường nét và điểm yên ngựa:

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

Chúng tôi kết thúc với một cốt truyện như thế này:

Answer 1

Không chắc nếu điều này (muộn) giúp, nhưng có vẻ như bạn đang quan tâm các điểm phân biệt, có nghĩa là, cả hai đa thức và phái sinh (wrt y) biến mất. Bạn có thể giải quyết hệ thống này cho {x, y} và vứt bỏ các giải pháp phức tạp như dưới đâ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}} 

Answer 2

Nếu bạn chỉ muốn âm mưu kết quả sau đó sử dụng StreamPlot[] trên gradient:

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] 

Bạn có thể phải loay hoay với độ chính xác của lô, StreamStyle và RegionFunction để hoàn thiện nó. Đặc biệt hữu ích sẽ được sử dụng các giải pháp cho sàn thung lũng để hạt giống StreamPoints lập trình.

Answer 3

Đã cập nhật: xem bên dưới.

tôi muốn tiếp cận này đầu tiên bằng cách hình dung những phần ảo của rễ:

này sẽ cho bạn ba điều ngay lập tức: 1) vào thư mục gốc đầu tiên luôn luôn là có thật, 2) thứ hai hai là các cặp liên hợp, và 3) có một khu vực nhỏ gần bằng không trong đó cả ba là có thật. Ngoài ra, lưu ý rằng các loại trừ chỉ đã thoát khỏi những điểm kỳ dị tại x=0, và chúng ta có thể thấy tại sao khi chúng ta phóng to:

Sau đó chúng tôi có thể sử dụng EvalutionMonitor để tạo ra danh sách các rễ trực tiếp:

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] 

(Lưu ý, các đặc điểm kỹ thuật Part là một chút kỳ quặc, Reap trả về một List về những gì đang gieo như mục thứ hai trong một List, vì vậy kết quả này trong một danh sách lồng nhau. Ngoài ra, 01.237.không lấy mẫu điểm một cách đơn giản, vì vậy cần có SortBy.) Có thể có một tuyến đường thanh lịch hơn để xác định nơi hai rễ cuối cùng trở nên phức tạp, nhưng vì phần ảo của chúng liên tục lặp đi lặp lại, nó có vẻ dễ dàng hơn buộc nó.

Chỉnh sửa: Vì bạn đã đề cập rằng bạn muốn một phương pháp tự động để tạo nơi một số rễ trở nên phức tạp, tôi đã khám phá những gì sẽ xảy ra khi bạn thay thế trong y -> p + I q. Bây giờ điều này giả định rằng x là có thật, nhưng bạn đã làm điều đó trong giải pháp của bạn. Cụ thể, tôi làm như sau:

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

nơi bước thứ hai cho phép tôi cách ly các phần thực và tưởng tượng của phương trình và đơn giản hóa chúng độc lập với nhau. Làm điều tương tự với đa thức 2D chung, f + d x + a x^2 + e y + 2 c x y + b y^2, nhưng làm cho cả phức hợp x và y; Tôi lưu ý rằng Im[poly] = Im[x] D[poly, Im[x]] + Im[y] D[poly,[y]] và điều này cũng có thể giữ cho phương trình của bạn. Bằng cách thực hiện x thực, phần ảo của poly trở thành q lần một số chức năng của x, p và q. Vì vậy, hãy đặt q=0 luôn cung cấp Im[poly] == 0. Nhưng, điều đó không cho chúng ta biết điều gì mới mẻ. Tuy nhiên, nếu chúng ta

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

chúng tôi nhận được một số công thức cho q liên quan đến x và p. Đối với một số giá trị của x và p, các công thức đó có thể là tưởng tượng và chúng tôi có thể sử dụng Reduce để xác định vị trí Re[qvals] == 0. Nói cách khác, chúng tôi muốn phần "tưởng tượng" của y trở thành hiện thực và điều này có thể được thực hiện bằng cách cho phép q bằng không hoặc hoàn toàn là tưởng tượng. Âm mưu khu vực nơi Re[q]==0 và bao phủ các đường extremal dốc qua

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

cho

trong đó khẳng định các khu vực trong hai lô đầu tiên cho thấy 3 rễ thật.

Answer 4

Đã kết thúc cố gắng bản thân mình vì mục tiêu thực sự là thực hiện 'tắt tay'. Tôi sẽ để câu hỏi mở trong một thời gian dài để xem liệu có ai tìm được cách tốt hơn không.

Mã bên dưới sử dụng bisection để gắn các điểm trong đó CountRoots giá trị thay đổi. Này làm việc cho trường hợp của tôi (đốm dị tại x = 0 là tinh khiết may mắn):

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

Thực hiện:

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