1
$\begingroup$

I wish to be able to scatter slightly (by at most some fixed $\varepsilon>0$) the zeros of a given complex polynomial $p(z)$ in a random way, and then see the effect on the critical points.

For example, I wish to pick $100$ monic complex polynomials of degree $5$, each of which has one zero within distance $1/100$ of the fifth roots of unity, and then plot the critical points of each of these polynomials. (The case from the first sentence with $p(z)=z^5-1$, and $\varepsilon=1/100$.)

I am relatively new to using Mathematica in this way, so any suggestions or pointers (or of course full solutions!) would be very helpful.

$\endgroup$
2
  • $\begingroup$ Before I write an answer that is not correct let me ask something: The five roots of z^5-1 can be constructed with e.g. Table[Exp[-2/5 Pi I i], {i, 0, 4}]. If you draw a random complex number that lies within the disk with radius epsilon, you can add them to each root. Then you can reconstruct the polynomial by multiplication (z-root1-rand1)(z-root2-rand2)... You can do this without explicitly specifying the random shifts. The critical points are the zeros of D[poly,z], right? These too can be calculated without values for the random shifts and you get an analytical solution for them. $\endgroup$
    – halirutan
    Mar 18, 2018 at 5:09
  • $\begingroup$ If above is what you want, then you have a solution for your critical points, where you only need to insert the random values for the shifts. Plotting them can be done by using Graphics, Point and ReIm and maybe with a nice ContourPlot beneath it. Is that what you seek? $\endgroup$
    – halirutan
    Mar 18, 2018 at 5:11

1 Answer 1

4
$\begingroup$

Borrowing DumpsterDoofus's routine from this answer, and safening it slightly:

With[{γ = 0.12, β = 1.},
     fLor = Compile[{{x, _Integer}, {y, _Integer}},
                    (γ/(γ + x^2 + y^2))^β, RuntimeAttributes -> {Listable}]];

PlotComplexPoints[list_, magnification_, paddingX_, paddingY_, brightness_, vec_] :=
    Module[{dimX, dimY, RePos, ImPos, lor, posf, sparse},
           posf = 1 + Round[magnification (# - Min[#])] &;
           RePos = paddingX + posf[Re[list]]; ImPos = paddingY + posf[Im[list]];
           dimX = paddingX + Max[RePos]; dimY = paddingY + Max[ImPos];
           With[{spopt = SystemOptions["SparseArrayOptions"]}, 
                Internal`WithLocalSettings[
                SetSystemOptions["SparseArrayOptions" -> {"TreatRepeatedEntries" -> 1}],
                Image[Outer[Times,
                            brightness Abs[InverseFourier[Fourier[
                            SparseArray[Thread[Transpose[{ImPos, RePos}] -> 
                                        ConstantArray[1, Length[list]]], {dimY, dimX}]]
                            Fourier[RotateRight[fLor[#[[All, All, 1]],
                                                     #[[All, All, 2]]] & @
                                    Outer[List, 
                                          Range[-Quotient[dimY, 2],
                                                Quotient[dimY - 1, 2]], 
                                          Range[-Quotient[dimX, 2], 
                                                Quotient[dimX - 1, 2]]],
                                          {Quotient[dimY, 2], Quotient[dimX, 2]}]],
                                    FourierParameters -> {-1, 1}]], 
                            Developer`ToPackedArray[N[vec]]], Magnification -> 1],
                SetSystemOptions[spopt]]]]

Generate randomly perturbed critical points given the roots:

makeRandomCriticalPoints[roots_?VectorQ, h_] := Block[{n = Length[roots], z}, 
    z /. NSolve[D[Product[z - zk,
                          {zk, roots + RandomReal[h, n] Exp[I RandomReal[2 π, n]]}], z],
                z]]

Finally:

penta = Exp[2 π I Range[0, 4]/5];

BlockRandom[SeedRandom[42]; (* for reproducibility *)
            With[{nPolys = 200, ε = 1/100},
                 PlotComplexPoints[Flatten[Table[makeRandomCriticalPoints[penta, ε],
                            {nPolys}]], 600, 20, 20, 10, {1., 0.3, 0.1}]]]

critical points of 200 polynomials with perturbed roots

A bonus picture:

BlockRandom[SeedRandom[42];
            With[{nPolys = 200, ε = 1/10},
                 PlotComplexPoints[Flatten[Table[makeRandomCriticalPoints[
                            Riffle[-penta/(1 + GoldenRatio), penta], ε], {nPolys}]],
                                   200, 20, 20, 10, {1., 0.3, 0.1}]]]

starburst

$\endgroup$
9
  • $\begingroup$ Omg. Bonus picture doesn't evaluate though: Thread::tdlen: Objects of unequal length in (...) cannot be combined. $\endgroup$ Mar 18, 2018 at 15:45
  • 1
    $\begingroup$ @Vsevolod, I had inadvertently copied an old version of makeRandomCriticalPoints[]. Please try it now. $\endgroup$ Mar 18, 2018 at 16:52
  • $\begingroup$ @J.M. Thanks very much. However I have a question. When I changed the definition of "penta" to do the same thing for $p(z)=z^{50}-1$, by replacing "Range[0, 4]/5" in the definition of "penta" with "Range[0, 49]/50", the behavior is very strange. Could you try this and see if you can explain the behavior? Note that when I do the same for $p(z)=z^{30}-1$, the picture is exactly what we expect. The strange behavior only begins when the number of roots approach $50$. $\endgroup$ Mar 20, 2018 at 1:48
  • $\begingroup$ I get a picture like this for PlotComplexPoints[Flatten[Table[makeRandomCriticalPoints[Exp[2 π I Range[0, 49]/50], 1/100], {100}]], 300, 20, 20, 10, {1., 0.3, 0.1}], @Trevor. I see nothing amiss; what were you expecting, exactly? $\endgroup$ Mar 20, 2018 at 3:31
  • $\begingroup$ @J.M. , Sorry I should have been more specific. That is exactly the picture I got and the one I was expecting. It is the second picture, which shows both the zeros and the critical points, where the behavior gets strange when the number of zeros approach 50. Sorry I don't know how to post pictures as you did (very helpfully), but it looks like some of the critical points scatter towards $1$, and some of the zeros scatter towards $0$. Please let me know if you see the same. $\endgroup$ Mar 21, 2018 at 20:13

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.