4
$\begingroup$

Does anyone have a code for drawing tropical amoebas? For example for the polynomial $$P(z,w) = 50z^3 + 83z^2w + 24zw^2 + w^3 + 392z^2 + 414zw + 50w^2 - 28z + 59w - 100$$ it has a form

amoeba

(The picture is taken from mathoverflow.)

One of the possible approaches is implemented here, but pictures are not very sharp.

EDT Greg Hurst gave in the comments a good code which works nicely for the curve mentioned above. It works good in other examples, but something goes wrong for the simple curve $1 + x_1 + x_2 - 4 x_1 x_2 + x_1^2 x_2^2=0$:

amoeba

A picture generated by amoebas.ru

amoeba

The code given by Greg Hurst:

poly[x1_, x2_] := 1 + x1 + x2 - 4 x1 x2 + x2^2 x1^2;
Clear[try1x, try2x];
{try1x[b_, \[Theta]_], try2x[b_, \[Theta]_]} = 
  Thread[{Log[Abs[SolveValues[poly[z, w] == 0, z]]] /. 
     w -> Exp[b + I \[Theta]], b}];

Clear[try1y, try2y];
{try1y[a_, \[Theta]_], try2y[a_, \[Theta]_]} = 
  Thread[{a, 
    Log[Abs[SolveValues[poly[z, w] == 0, w]]] /. 
     z -> Exp[a + I \[Theta]]}];

showx = ParallelMap[
    ParametricPlot[#, {\[Theta], 0, 2 \[Pi]}, PlotPoints -> 500, 
      MaxRecursion -> 0] &, 
    Join @@ Table[{try1x[b, \[Theta]], try2x[b, \[Theta]]}, {b, -4, 
       6, .01}]]; // AbsoluteTiming

showy = ParallelMap[
    ParametricPlot[#, {\[Theta], 0, 2 \[Pi]}, PlotPoints -> 500, 
      MaxRecursion -> 0] &, 
    Join @@ Table[{try1y[a, \[Theta]], try2y[a, \[Theta]]}, {a, -4, 
       6, .01}]]; // AbsoluteTiming

plot = Show[showx, showy, PlotRangePadding -> None, 
   PlotRange -> {{-4, 6}, {-4, 6}}] /. {Line[{a_, ___, b_}] :> 
    Line[{a, b}]}

im = Rasterize[Show[plot, Axes -> False], RasterSize -> 720];

ImageMesh[ColorNegate[im], DataRange -> {{-4, 6}, {-4, 6}}]
Improve this question
$\endgroup$
14
  • $\begingroup$ Probably not the best way in terms of speed or quality, but maybe this will help. poly[z_, w_] := 50*z^3 + 83*z^2*w + 24*z*w^2 + w^3 + 392*z^2 + 414*z*w + 50*w^2 - 28*z + 59*w - 100; solnzlogs[w_?NumberQ] := Thread[{Log[Abs[NSolveValues[poly[z, w] == 0, z]]], Log[Abs[w]]}] bnd = 20.; skip = .01; AbsoluteTiming[ pts = Flatten[Table[solnzlogs[x + I*y], {x, -bnd, bnd, skip}, {y, -bnd, bnd, skip}], 2];] ListPlot[pts, PlotRange -> All] I expect this to not come back very soon though. $\endgroup$ Commented Aug 28 at 3:55
  • $\begingroup$ I'm a bit surprised this was closed without anyone requesting further details. Like "Is this a Mathemqtica question?" or "What have you tried thus far?" $\endgroup$ Commented Aug 28 at 3:56
  • 1
    $\begingroup$ FWIW here's a screengrab of a not-so elegant solution. I used 28 kernels, so on a more modest machine this will take awhile to run (~15 min in serial). If this question is reopened, I'll post this as an answer (with an explanation of how the code works). i.sstatic.net/OIHcDV18.png $\endgroup$ Commented Aug 28 at 14:10
  • 3
    $\begingroup$ @GregHurst Why reinventing the wheel? mathematica.stackexchange.com/questions/85081/… mathematica.stackexchange.com/questions/295955/… $\endgroup$ Commented Aug 28 at 22:22
  • 1
    $\begingroup$ @azerbajdzan Aha, I suppose this could be closed by being a duplicate question. It was closed for a different reason before I believe. Though the best plot I can make after fiddling with the solutions on the other page is i.sstatic.net/tCus1Ery.png. $\endgroup$ Commented Aug 29 at 2:16

1 Answer 1

Reset to default
7
$\begingroup$

This is not the most efficient solution, but it seems to get the job done.

Let's start with our polynomial:

poly[z_, w_] := 
  50*z^3 + 83*z^2*w + 24*z*w^2 + w^3 + 392*z^2 + 414*z*w + 50*w^2 - 28*z + 59*w - 100;

By definition of the amoeba (here), we can represent it in Mathematica by something like the union of the sets

Thread[{Log[Abs[SolveValues[poly[z, w] == 0, z]]], Log[Abs[w]]}]

We can simplify things for the plotter by setting Log[Abs[w]] == b, which implies w == Exp[b + I θ] for some θ between 0 and 2π. This will let us represent our amoeba as

aparam1 = Thread[{Log[Abs[SolveValues[poly[z, w] == 0, z]]] /. w -> Exp[b + I θ], b}];

So in theory we then can plot over b and θ, just as this and this answer. However the best plot I can get is something like

ParametricPlot[aparam1, {b, -4, 6}, {θ, 0, 2π}, 
  PlotStyle -> Directive[ColorData[110, 1], Opacity[1]], 
  BoundaryStyle -> None, 
  PlotRange -> {{-4, 4}, {-4, 6}}, 
  PlotPoints -> 200
]

It's pretty good, but the exclusions handling (which we rely on) introduce some defects and we also don't get the full boundary in case we'd like to style it different etc.

Instead we can fix an explicit value of b and get a 1D plot over θ. We do this for many values of b and combine the plots.

There are of course features we will miss by plotting over a fixed increment of b. We can catch those by also plotting in the orthogonal direction, in the form {a, Log[Abs[root[a, θ]]]}.

aparam2 = Thread[{a, Log[Abs[SolveValues[poly[z, w] == 0, w]]] /. z -> Exp[a + I θ]}];

It's not the most efficient way, especially because we need a lot of plot points. Also note I played with settings to get the exclusions engine to not make any mistakes.

showx = ParallelMap[
  ParametricPlot[#, {θ, 0, 2π}, PlotPoints -> 500, MaxRecursion -> 0, 
    PlotHighlighting -> None, PlotStyle -> ColorData[110, 1]]&,
  Join @@ Table[aparam1, {b, -4, 6, .01}]
];

showy = ParallelMap[
  ParametricPlot[#, {θ, 0, 2π}, PlotPoints -> 200, MaxRecursion -> 1, 
    PlotHighlighting -> None, PlotStyle -> ColorData[110, 1]]&,
  Join @@ Table[aparam2, {a, -4, 4, .01}]
];

We can then show things together:

plot = Show[showx, showy, PlotRangePadding -> None, PlotRange -> {{-4, 4}, {-4, 6}}]

We still don't have the boundary though. We can infer one by discretizing an image of this plot:

im = Rasterize[Show[plot, Axes -> False], RasterSize -> 1440];

Show[ImageMesh[ColorNegate[im], MeshCellStyle -> {1 -> Black}, DataRange -> {{-4, 4}, {-4, 6}}], Axes -> True]

Edit

It turns out to be about 12 times faster to combine plots over argument instead of magnitude:

Show[
 ParallelTable[
  Quiet @ ParametricPlot[
    Evaluate[Thread[Log[Abs[{z, SolveValues[poly[z, w] == 0, w]}]] /. z -> Exp[θ I] z]], 
    {z, Exp[-4], Exp[4]}, 
    PlotPoints -> 50, PlotRange -> {{-4, 4}, {-4, 6}}, 
    AspectRatio -> Automatic, PlotStyle -> ColorData[110, 1], 
    PlotHighlighting -> None
  ],
  {θ, 0., 2π, π/55}
 ],
 ParallelTable[
  Quiet @ ParametricPlot[
    Evaluate[Thread[Log[Abs[{SolveValues[poly[z, w] == 0, z], w}]] /. w -> Exp[θ I] w]], 
    {w, Exp[-4], Exp[6]}, 
    PlotPoints -> 100, PlotRange -> {{-4, 4}, {-4, 6}}, 
    AspectRatio -> Automatic, PlotStyle -> ColorData[110, 1], 
    PlotHighlighting -> None
  ],
  {θ, 0., 2π, π/55}
 ]
]

This can then also be passed to ImageMesh to extract a boundary too.

Improve this answer
$\endgroup$
5
  • $\begingroup$ Can you plot it with PlotRange -> {{-4, 4}, {-5, 6}}]? I am not trying to use your code as I suspect it can freeze my PC. $\endgroup$ Commented Aug 30 at 8:59
  • $\begingroup$ Sure, here you go: i.sstatic.net/YjCFL9Dx.png $\endgroup$ Commented Aug 31 at 3:00
  • $\begingroup$ Because the left downward tentacle is in fact two tentacles close to each other, so I was curious whether in your plots they are distinguishable. From the plots it is not evident but maybe just because they are not zoomed enough. $\endgroup$ Commented Aug 31 at 9:20
  • $\begingroup$ It could be that the lines are too thick to show fine features like that. $\endgroup$ Commented Aug 31 at 20:55
  • 1
    $\begingroup$ @azerbajdzan If I take that same plot I commented above, and plot it with thickness AbsoluteThickness[.1] I get this plot, which does show that there are indeed two tentacles: i.sstatic.net/fzHU1b76.png. If I sample over θ at 1/10th spacing and keep AbsoluteThickness[.1] we get the plot you're after: i.sstatic.net/bZr2GYfU.png $\endgroup$ Commented Aug 31 at 22:30

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.