Plotting the Star of Bethlehem

Question

Christmas is around the corner! In preparation, consider the polynomial

f[x_,y_]:= 1/10 (4 - x^2 - 2 y^2)^2 + 1/2 (x^2 + y^2) 

and define $r=(f_x,f_y)$ using the first-order partial derivatives such that

r[x_,y_]:={x (-3 + 2 x^2 + 4 y^2)/5, y (-11 + 4 x^2 + 8 y^2)/5} 

and $F=f_{xx}\cdot f_{yy}-f_{xy}^2$ using the second-order partial derivatives such that

F[x_,y_]:= (33 + 24 x^4 - 116 y^2 + 96 y^4 - 78 x^2 + 96 x^2 y^2)/25 

The question is how to plot the following 'implicit' function efficiently and sharp in Mathematica?

$$Z(s,t)= \sum_{(x,y)\in \mathbb{R}^2:\, r(x,y)=(s,t)} \frac{1}{\left| F(x,y) \right|}$$

It should look like this in $[-1,1]^2$: (when black indicates zero and white large values)

Actually, the fastes contributions are from

• @xzczd (backward way, tracing $r^{-1}$)
• @George Varnavides (forward way, tracing $r$)

Most elegant solution from

• @Roman

Addendum: If you successfully plotted the star you may try the following function as input

f=ListInterpolation[RandomReal[2,{9,9}]+Outer[#1^2+#2^2&, Range[9], Range[9]]/2];

leading to $Z$ looking like this:

Answer 1

Hint: $Z$ is the Jacobian for the transformation between $(x,y)$ and $(s,t)$ coordinates. No need to even define $Z$, no need to invert polynomials, no branch cuts.

f[x_, y_] = 1/10 (4 - x^2 - 2 y^2)^2 + 1/2 (x^2 + y^2);
r[x_, y_] = D[f[x, y], {{x, y}}];
With[{span = 2, step = 0.001, binsize = 1/100},
  DensityHistogram[
    Join @@ Table[r[x, y], {x, -span, span, step}, {y, -span, span, step}],
    {-1, 1, binsize},
    ColorFunction -> GrayLevel]]

Same idea but much faster (takes 1.9 seconds on my laptop):

With[{span = 1.72, step = 0.001, binsize = 1/100},
  ArrayPlot[
    Transpose@BinCounts[Transpose[
      r @@ Transpose[Tuples[Range[-span, span, step], 2]]], {-1,1,binsize}, {-1,1,binsize}],
    ColorFunction -> GrayLevel]]

where the number 1.72 is Root[-5-3#+2#^3&, 1], as gleaned from Reduce[Thread[-1 <= r[x,y] <= 1], {x,y}].

Addendum

This method works with the addendum question as well:

f = ListInterpolation[
      RandomReal[2, {9, 9}] + Outer[#1^2 + #2^2 &, Range[9], Range[9]]/2];
r[x_, y_] = D[f[x, y], {{x, y}}];
With[{step = 0.005, binsize = 1/10},
  ArrayPlot[Transpose@BinCounts[Transpose[
    r @@ Transpose[Tuples[Range[1, 9, step], 2]]],
    {0, 10, binsize}, {0, 10, binsize}],
  ColorFunction -> GrayLevel]]

I'll leave it to skillful experts to make prettier plots. Here I'm only addressing the question of plotting efficiency.

Answer 2

The noise in the picture given by OP seems to suggest random number has been used, so here's my trial:

data = RandomReal[{-2, 2}, {2, 4 10^5}];
s = 1000;
func = Function[{x, y}, Round@Rescale[r[x, y], {-1, 1}, {1, s}] -> 1/Abs@F[x, y], 
   Listable];
stFlst = func @@ data; // AbsoluteTiming
(* {5.87326, Null} *)
SetSystemOptions["SparseArrayOptions" -> {"TreatRepeatedEntries" -> 1}];
array = SparseArray[
    DeleteCases[stFlst, ({x_, y_} -> _) /; ! (1 < x < s && 1 < y < s)], {s, 
     s}]; // AbsoluteTiming
(* {0.655598, Null} *)
ArrayPlot[array // Transpose, PlotRange -> {0, 10}, 
 ColorFunction -> "AvocadoColors"]

---

The following is a fast smooth plot, based on the discussion here:

Clear[s, x, y]
sol = {x, y} /. Solve[r[x, y] == {s, t}, {x, y}];
symbolicroots = Cases[sol, a__Root, Infinity] // Union
coef = CoefficientList[symbolicroots[[1, 1, 1]], #]
{x, y} = sol[[1]] /. __Root -> ro
nroots[c_List] := Block[{a}, a = DiagonalMatrix[ConstantArray[1., Length@c - 2], 1];
   a[[-1]] = -Most@c/Last@c;
   Eigenvalues[a]];
check = Function @@ {If[Abs@Im@# < 10^-6, 1, 0] // PiecewiseExpand // 
    Simplify`PWToUnitStep}
eps = 10^-6; dt = 2 10^-3; 
 data = ParallelTable[
   Block[{ro = Pick[#, check@#, 1] &@nroots[coef]}, 
    1/F[x, y] // Abs // Total], {s, -1 + eps, 1 + eps, dt}, {t, -1 + eps, 1 + eps, 
    dt}]; // AbsoluteTiming
(* {34.7344, Null}, dual core laptop *)
ArrayPlot[data // Transpose, ColorFunction -> "AvocadoColors", PlotRange -> {0, 30}]

---

Solution to Addendum

The method above for smooth plot no longer works well on the f in addendum, but the first method using random number still works. I've used InterpolationToPiecewise to transform the InterpolatingFunction to something compilable.

SeedRandom["star"]; 
dataf = RandomReal[2, {9, 9}] + Outer[#1^2 + #2^2 &, Range[9], Range[9]]/2;
InterpolationToPiecewise[if_, x_] := Module[{main, default, grid}, grid = if["Grid"];
   Piecewise[{if@"GetPolynomial"[#, x - #], x < First@#} & /@ grid[[2 ;; -2]], 
   if@"GetPolynomial"[#, x - #] &@grid[[-1]]]] /; if["InterpolationMethod"] == "Hermite"
poly[x_, y_] = 
  InterpolationToPiecewise[
   ListInterpolation[InterpolationToPiecewise[ListInterpolation[#], y] & /@ dataf], x];
compile = Compile[{x, y}, #, RuntimeAttributes -> {Listable}] &;
r = compile@D[poly[x, y], {{x, y}}];
absdivideF = 
  Abs@Divide[1, (D[#, x, x] D[#, y, y] - D[#, x, y]^2)] &@poly[x, y] // compile;
domain = {0, 10};
{datax, datay} = RandomReal[domain, {2, 10^6}];
s = 1000;
func = Function[{x, y}, 
  Clip[Round@Rescale[r[x, y], domain, {1, s}], {1, s}] -> absdivideF[x, y]];
stFlst = func[datax, datay]; // AbsoluteTiming
(* {2.93232, Null} *)
SetSystemOptions["SparseArrayOptions" -> {"TreatRepeatedEntries" -> 1}];
array = SparseArray[stFlst, {s, s}];
ArrayPlot[array // Transpose, PlotRange -> {0, 2}, ColorFunction -> "AvocadoColors"]

Remark

If you have a C compiler installed, consider using the following definition of compile:

compile = Last@
    Compile[{{xlst, _Real, 1}, {ylst, _Real, 1}}, 
     Block[{x, y}, Table[x = xlst[[i]]; y = ylst[[i]]; #, {i, Length@xlst}]], 
     CompilationTarget -> "C", RuntimeOptions -> "Speed"] &;

It will vastly speed up the generation of stFlst, at the expense of long compilation time of r and absdivideF.

Answer 3

A compiled version of @Roman's idea (binning code adapted from this answer):

starCompiled = 
 Compile[{{xyspan, _Real}, {delta, _Real}, {binspan, _Real}, \
{binsize, _Real}}, 
  Block[{bins, dimx, dimy, x, y, tx, ty}, 
   bins = Table[
     0, {Floor[binspan 2/binsize] + 1}, {Floor[binspan 2/binsize] + 
       1}];
   {dimx, dimy} = Dimensions[bins];
   Do[{x, y} = {xx - 2/5 xx (4 - xx^2 - 2 yy^2), 
      yy - 4/5 yy (4 - xx^2 - 2 yy^2)};
    tx = Floor[(x + binspan)/binsize] + 1;
    ty = Floor[(y + binspan)/binsize] + 1;
    If[tx >= 1 && tx <= dimx && ty >= 1 && ty <= dimy, 
     bins[[tx, ty]] += 1], {xx, -xyspan, xyspan, delta}, {yy, -xyspan,
      xyspan, delta}];
   bins], CompilationTarget -> "C", RuntimeOptions -> "Speed"]

I seem to get a ~3.5x speedup on my machine, with this relatively sharp image taking ~1.5 seconds:

ArrayPlot[starCompiled[1.72,0.0005,1.,0.005]^\[Transpose],Frame->False,ColorFunction->"AvocadoColors"]

Addendum

For completeness, here's a compiled version of the binning along idea on the addendum. Credit for the compiled interpolating function goes to @xzczd and reference therein.

dataf = RandomReal[2, {9, 9}] + 
   Outer[#1^2 + #2^2 &, Range[9], Range[9]]/2;
InterpolationToPiecewise[if_, x_] := 
 Module[{main, default, grid}, grid = if["Grid"];
   Piecewise[{if@"GetPolynomial"[#, x - #], x < First@#} & /@ 
     grid[[2 ;; -2]], if@"GetPolynomial"[#, x - #] &@grid[[-1]]]] /; 
  if["InterpolationMethod"] == "Hermite"
poly[x_, y_] = 
  InterpolationToPiecewise[
   ListInterpolation[
    InterpolationToPiecewise[ListInterpolation[#], y] & /@ dataf], 
   x];
compile = 
  Compile[{x, y}, #, RuntimeAttributes -> {Listable}, 
    CompilationTarget -> "C", RuntimeOptions -> "Speed"] &;
r = compile@D[poly[x, y], {{x, y}}];
addendumCompiled = 
 Compile[{{xmin, _Real}, {xmax, _Real}, {ymin, _Real}, {ymax, _Real}, \
{delta, _Real}, {binsize, _Real}}, 
  Block[{bins, dimx, dimy, x, y, tx, ty}, 
   bins = Table[
     0, {Floor[(xmax - xmin)/binsize] + 
       1}, {Floor[(ymax - ymin)/binsize] + 1}];
   {dimx, dimy} = Dimensions[bins];
   Do[{x, y} = r[xx, yy];
    tx = Floor[(x - xmin)/binsize] + 1;
    ty = Floor[(y - ymin)/binsize] + 1;
    If[tx >= 1 && tx <= dimx && ty >= 1 && ty <= dimy, 
     bins[[tx, ty]] += 1], {xx, xmin, xmax, delta}, {yy, ymin, ymax, 
     delta}];
   bins], CompilationTarget -> "C", RuntimeOptions -> "Speed", 
  CompilationOptions -> {"InlineExternalDefinitions" -> True}]

Seems to be relatively fast (and CompilePrint confirms no calls to MainEvaluator):

AbsoluteTiming[bins = addendumCompiled[1, 9, 1, 9, 0.0005, 0.005];]
Image@Rescale@Log[bins + 1]

PS

Potentially (un)related (any politics aside), here's the same idea on a complex strange attractor given by the equation (from the book Symmetry in Chaos: A Search for Pattern in Mathematics, Art and Nature): $$ F(z) = \left(\lambda + \alpha z z^* + \beta \Re(z^n) + \delta \Re\left[\left(\frac{z}{|z|}\right)^{np}\right]|z|\right)z + \gamma z^{*n-1} $$

attractorNonLinear = 
 Compile[{{xmin, _Real}, {xmax, _Real}, {ymin, _Real}, {ymax, _Real}, \
{delta, _Real}, {itmax, _Integer}, {n, _Integer}, {lambda, _Real}, \
{a, _Real}, {b, _Real}, {c, _Real}, {d, _Real}, {p, _Integer}}, 
  Block[{bins, dim, x, y, tx, ty, z, b1, radii, normed, coordinates},
   bins = 
    Table[0, {Floor[(xmax - xmin)/delta] + 
       1}, {Floor[(ymax - ymin)/delta] + 1}];
   dim = Dimensions[bins];
   z = -0.3 + 0.2 I;
   Do[z = (lambda + a z Conjugate[z] + b Re[z^n] + 
         d Re[(z/Abs[z])^n p] Abs[z]) z + c Conjugate[z]^(n - 1);
    {x, y} = {Re[z], Im[z]};
    tx = Floor[(x - xmin)/delta] + 1;
    ty = Floor[(y - ymin)/delta] + 1;
    If[tx >= 1 && tx <= dim[[1]] && ty >= 1 && ty <= dim[[2]], 
     bins[[tx, ty]] += 1], {i, 1, itmax}];
   bins], CompilationTarget :> "C", RuntimeOptions -> "Speed"]

which gives the Star of David with the following parameter values:

AbsoluteTiming[bins=N[attractorNonLinear[-1.6,1.6,-1.6,1.6,0.001,5 10^7,6,-2.42,1.0,-0.04,0.14,0.088,0]];]
ArrayPlot[Log[bins+1],ColorFunction->"AvocadoColors",Frame->False]

Answer 4

Some code, built for speed, building on @Roman's idea of projecting r[x, y]:

rr[x_, y_] := {x (-3 + 2 x^2 + 4 y^2)/5, y (-11 + 4 x^2 + 8 y^2)/5};
FF[x_, y_] := (33 + 24 x^4 - 116 y^2 + 96 y^4 - 78 x^2 + 96 x^2 y^2)/
   25;
PrintTemporary@Dynamic@{Clock@Infinity};
Block[{n, x, y, r, F, classes, clusters},
   {r, F} = {rr[x, y], FF[x, y]} /. {Power[v_, n_] :> v[n], 
      v : x | y :> v[1]};
   n = 1000;
   (* x,y coordinates )
   {x[1], y[1]} = Transpose@Tuples[Subdivide[-1.8, 1.8, n], 2];
   ( x,y dithering )
   {x[1], y[1]} += RandomReal[2/n {-1, 1}, Dimensions@{x[1], y[1]}];
   ( powers of x,y for r,F )
   {x[2], y[2]} = {x[1], y[1]}^2;
   {x[4], y[4]} = {x[2], y[2]}^2;
   ( opacity, point size depend on class = magnitude(F) *)
   classes = Round[Abs[F], 5];
   clusters = GatherBy[Range[(n + 1)^2], classes[[#]] &];
   star = Graphics[{White,
      GraphicsComplex[
       r // Transpose,
       With[{class = classes[[First@#]]},
          With[{s = (4. + 4. class^1.5)/n},
           {PointSize[s], Opacity[0.0004/n/s^2], Point@#}
           ]] & /@ clusters]
      },
     Background -> Black, PlotRange -> 1.1]
   ]; // AbsoluteTiming
(*  {0.225483, Null}  *)

There's a Moiré effect if you leave out the dithering (second {x[1], y[1]} line), which may be seen as pretty by some. It can be changed somewhat by manipulating the size parameter s. The size and opacity do not quite scale with n, and I gave trying to make it work exactly. Calculating n = 2000 took only 1 sec. Rendering n = 2000 led to a crash. I tried it twice, because my laptop often crashes for other reasons; but two crashes suggest n = 2000 exceeds the limits of my machine. Here is n = 1000:

img = Image@star; // AbsoluteTiming
img
(*  {2.4411, Null}  *)
star (* renders nicer than img, but takes a little longer *)

---

Interpolating functions are slower, but the operate the same way as other functions:

ff = BlockRandom[
   ListInterpolation[
     RandomReal[2, {9, 9}] + 
      Outer[#1^2 + #2^2 &, Range[9], Range[9]]/2][x, y],
   RandomSeeding -> 0];
rr[x_, y_] = D[ff, {{x, y}}];
FF[x_, y_] = Det@D[ff, {{x, y}, 2}];
PrintTemporary@Dynamic@{Clock@Infinity};
Block[{n, x, y, r, F, classes, clusters, pr},
   {r, F} = {rr[x, y], FF[x, y]};
   n = 600;
   {x, y} = Transpose@Tuples[Subdivide[1., 9., n], 2];
   {x, y} += RandomReal[2/n {-1, 1}, Dimensions@{x[1], y[1]}];
   {x, y} = Clip[{x, y}, {1., 9.}];
classes = Round[Abs[F], 5];
   clusters = GatherBy[Range[(n + 1)^2], classes[[#]] &];
   r = r;
   pr = MinMax /@ r;
   star = Graphics[{White,
      GraphicsComplex[
       r // Transpose,
       With[{class = classes[[First@#]]},
          With[{s = (2. + 2. class^1.8)/n},
           {PointSize[s], Opacity[0.0003/n/s^2], Point@#}
           ]] & /@ clusters]
      },
     Background -> Black, PlotRange -> pr]
   ]; // AbsoluteTiming
(*  {14.0745, Null}  *)
PrintTemporary@Dynamic@{Clock@Infinity};
(img0 = Image[star];
  img = ImageApply[#^{3, 2, 1/3} &, img0];) // AbsoluteTiming
img
(*  {20.255, Null}  *)

Answer 5

I'll show a picture for now, and add code after others have had time to play around with this.

--- edit ---

[Adding code as promised]

First the input polynomials.

poly[x_, y_] := (33 + 24 x^4 - 116 y^2 + 96 y^4 - 78 x^2 + 
    96 x^2 y^2)/25
r[x_, y_] := {x (-3 + 2 x^2 + 4 y^2)/5, y (-11 + 4 x^2 + 8 y^2)/5}

Now create a function z that is the reciprocal of poly, but expressed in terms of {s,t} coordinates. We get this via a GroebnerBasis computation that eliminates {x,y} and gives the reciprocal. (Remark: This method will not work on the addendum problem.) An advantage here is that we get a function directly in terms of the new coordinates and thus do not need to consider the solutions in terms of {x,y}.

gb = GroebnerBasis[
   Flatten[{{s, t} - r[x, y], z*poly[x, y] - 1}], {z, s, t}, {x, y}];

Construct a function that sums the absolute value over the roots of our reciprocal.

new[s_, t_] = RootSum[(Evaluate[gb[[1]]] /. z -> #) &, Abs[#] &];

And a similar one that only sums the real roots.

newReal[s_, t_] = 
  RootSum[(Evaluate[gb[[1]]] /. z -> #) &,
  (1 - Abs[Sign[Im[#]]])*Abs[#] &];

--- end edit ---

ContourPlot[new[s, t], {s, -1.2, 1.2}, {t, -1.6, 1.6}, 
 ColorFunction -> GrayLevel, Contours -> {Automatic, 120}]

One alternative is to use a list plot.

n = 100000;
randata = With [{pts = RandomReal[{-1, 1}, {n, 2}]},
   Join[pts, Transpose[{Apply[new, pts, 1]}], 2]];
ListDensityPlot[randata, ColorFunction -> GrayLevel]

No particular improvement in artistry, I'm afraid.

I add a version that enforces that the roots in question be real-valued. Plot settings dutifully cribbed from response by @BobHanlon.

DensityPlot[newReal[s, t], {s, -1, 1}, {t, -1, 1}, PlotPoints -> 75, 
 MaxRecursion -> 3, ColorFunction -> GrayLevel]

Here is an alternative rendition as proposed by @MichaelE2.

DensityPlot[newReal[s, t], {s, -1., 1.}, {t, -1., 1.}, 
 PlotPoints -> 75, MaxRecursion -> 4, ColorFunctionScaling -> False, 
 ColorFunction -> (GrayLevel[ArcTan[0.25 #]/(Pi/2)] &), 
 PlotRange -> All]

Here is one these contours.

ContourPlot[newReal[s, t] == 3, {s, -1., 1.}, {t, -1., 1.}, MaxRecursion -> 4]

And here is one rotation, to get the 3D requested by @dominic (I realize the quality is poor).

star3D[s_, t_, u_] := newReal[s, t^2 + u^2]

Another rotation looks somewhat like Saturn. In deference to the upcoming solstice planetary event I will not show it; just pretend it is being blocked by Jupiter.

After clarification that requests a surface plot, here is one variant. I sum over all (real and complex) roots in this variant, and also take a square root so as to fit most of the plot in the default range by lowering the peaks. I like the effect from the avocado color scheme others have used,

Plot3D[new[s, t]^(1/2), {s, -1, 1}, {t, -1, 1}, PlotPoints -> 150, 
 MaxRecursion -> 4, ColorFunction -> "AvocadoColors"]

Note that it is somewhat "cuspy". Also there are some spiky parts that get clipped. To see them (and lose much of the rest) one can use PlotRange->All. If one uses only real roots the plot is not so nice. I suspect this is due to discontinuities wherein complex-valued roots merge into real values on some curve in the s-t plane (technically, it's a part of the discriminant variety).

Answer 6

Clear["Global`*"]

f[x_, y_] := 1/10 (4 - x^2 - 2 y^2)^2 + 1/2 (x^2 + y^2)

r[x_, y_] = D[f[x, y], {{x, y}}] // Simplify;

F[x_, y_] = 
  D[f[x, y], x, x]*D[f[x, y], y, y] - D[f[x, y], x, y]^2 // Simplify;

Z[s_?NumericQ, t_?NumericQ, max_ : Infinity] := 
  Total[1/Abs[
     F @@@ ({x, y} /.
        NSolve[Thread[r[x, y] == {s, t}], {x, y}])]];

(dp1 = DensityPlot[
     Z[s, t], {s, -1.2, 1.2}, {t, -1.2, 1.2},
     PlotPoints -> 75,
     MaxRecursion -> 3,
     ColorFunction -> GrayLevel,
     Frame -> False,
     ImageSize -> Medium]) // AbsoluteTiming // Column

You can sharpen the edges by rescaling the color function

(dp2 = DensityPlot[
     Z[s, t], {s, -1.2, 1.2}, {t, -1.2, 1.2},
     PlotPoints -> 75,
     MaxRecursion -> 3,
     ColorFunction -> (GrayLevel[Max[0, Log[#]]] &),
     Frame -> False,
     ImageSize -> Medium]) // AbsoluteTiming // Column

EDIT: Revised since the OP was edited to restrict {x, y} to reals.

Z[s_?NumericQ, t_?NumericQ, max_ : Infinity] := 
  Total[1/Abs[
     F @@@ ({x, y} /. NSolve[Thread[r[x, y] == {s, t}], {x, y}, Reals])]];
(dp1 = DensityPlot[Z[s, t], {s, -1.2, 1.2}, {t, -1.2, 1.2}, PlotPoints -> 75, 
     MaxRecursion -> 3, ColorFunction -> GrayLevel, Frame -> False, 
     ImageSize -> Medium]) // AbsoluteTiming // Column

With sharp edges

(dp2 = DensityPlot[Z[s, t], {s, -1.2, 1.2}, {t, -1.2, 1.2}, PlotPoints -> 75, 
     MaxRecursion -> 3, ColorFunction -> (GrayLevel[Max[0, Log[#]]] &), 
     Frame -> False, ImageSize -> Medium]) // AbsoluteTiming // Column

Answer 7

Using Roman's method to get an array of bin counts, using Image or Raster + Graphics gives better pictures and faster:

bincounts = With[{span = 1.72, step = 10^-3}, 
   Transpose @ BinCounts[Transpose[r @@ Transpose[Tuples[Range[-span, span, step], 2]]], 
  {-1, 1, 5 step}, {-1, 1, 5 step}]];

Image

Image  @ Rescale @ bincounts

To change the default color space ("Grayscale") colors, we can map desired colors on matrix entries:

Style[#, ImageSizeMultipliers -> {1, 1}] & @ Grid @ 
  Partition[Image[Map[#, Rescale @ bincounts, {-1}]] & /@ 
    {ColorData["AvocadoColors"], ColorData["DeepSeaColors"], 
     ColorData["SolarColors"], Blend[{Black, Red, White}, #] &}, 2]

We can also use ImageMultiply:

Style[#, ImageSizeMultipliers -> {1, 1}] & @ Grid @
 Partition[ImageMultiply[Image[Rescale @ bincounts], #] & /@ 
  {Red, Yellow, Green, Magenta}, 2]

Raster + Graphics

Graphics @ Raster @ Rescale @ bincounts

With Raster we can use the option ColorFunction with built-in and custom color functions:

Style[#, ImageSizeMultipliers -> {1, 1}] & @ Grid @ 
 Partition[Graphics[Raster[Rescale@bincounts, ColorFunction -> #]] & /@ 
   {ColorData["AvocadoColors"], ColorData["DeepSeaColors"], 
    ColorData["SolarColors"], Blend[{Black, Red, White}, #] &}, 2]

Answer 8

The question is how to plot the following 'implicit' function efficiently and sharp in Mathematica?

I am entirely sure that the code below is not what OP wants, but this question did make me to write a document describing an elaborated approach of using several types of Machine Learning (ML) workflows. See this post or this flow chart .

This code and results are from the "Simplistic approach" section of the post (it is better to try it in a new notebook):

ResourceFunction["DarkMode"][True]
lsStarReferenceSeeds = 
  DeleteDuplicates@{2021, 697734, 227488491, 296515155601, 328716690761, 
    25979673846, 48784395076, 61082107304, 63772596796, 128581744446, 
    254647184786, 271909611066, 296515155601, 575775702222, 595562118302, 
    663386458123, 664847685618, 680328164429, 859482663706};
Multicolumn[
 Table[BlockRandom[
    ResourceFunction["RandomMandala"]["RotationalSymmetryOrder" -> 2, 
     "NumberOfSeedElements" -> Automatic, 
     "ConnectingFunction" -> FilledCurve@*BezierCurve], RandomSeeding -> rs], {rs, lsStarReferenceSeeds}] /. GrayLevel[0.25`] -> White, 6, 
 Appearance -> "Horizontal"]

Here is a link to the corresponding notebook.

Answer 9

Here's my surface plot of the star showing analytically the magnitude of F(s,t). Note in particular I am summing both real and complex roots. Maybe could use the height to then color-code the 2D plots above.

r[x_, y_] = D[f[x, y], {{x, y}}] // Simplify
theF[{x_, 
   y_}] := (33 + 24 x^4 - 116 y^2 + 96 y^4 - 78 x^2 + 96 x^2 y^2)/25
mySumFun[s_?NumericQ, t_?NumericQ] := Module[{mySol},
   mySol = {x, y} /. 
     NSolve[{1/5 x (-3 + 2 x^2 + 4 y^2) == s, 
       1/5 y (-11 + 4 x^2 + 8 y^2) == t}, {x, y}];
   Plus @@ ((1/Abs[theF[#]]) & /@ mySol)
   ];
p1 = Plot3D[mySumFun[s, t], {s, -1, 1}, {t, -1, 1}, PlotRange -> 20, 
  ClippingStyle -> None, BoxRatios -> {1, 1, 1},PlotPoints->100]
GraphicsRow[{p1, p1}]

And here is the under-side which looks more like a star. Might be a challenge to remove the four corners and reveal only the imbedded star-shape.

Update: My method to extract the star-shape is to define a BoundaryMeshRegion that I can pass to Plot3D. In order to do this, I first create a color-coded root-map of the region: LightGreen: One real root, DarkGreen: Three real roots, Brown: Five real roots. As you can see in the plot, this gives an easy way to delimit the region: Simply go through the array and pick out where the lightgreen changes to darkgreen. This will then give the boundary of the star-shape.

mySumFun[s_?NumericQ, t_?NumericQ] := Module[{mySol},
   mySol = {x, y} /. 
     NSolve[{1/5 x (-3 + 2 x^2 + 4 y^2) == s, 
       1/5 y (-11 + 4 x^2 + 8 y^2) == t}, {x, y}];
   Plus @@ ((1/Abs[theF[#]]) & /@ mySol)
   ];
myCountFun[{s_?NumericQ, t_?NumericQ}] := Module[{mySol},
   mySol = {x, y} /. 
     NSolve[{1/5 x (-3 + 2 x^2 + 4 y^2) == s, 
       1/5 y (-11 + 4 x^2 + 8 y^2) == t}, {x, y}];
   {{s, t}, Count[mySol, {x_, y_} /; Element[x, Reals]]}
   ];
matrixSize = 100;
xMin = -1.;
xMax = 1.;
yMin = -1.;
yMax = 1.;
plotRange6 = {{xMin, xMax}, {yMin, yMax}};
pointData = 
  Array[{#2, #1} &, {matrixSize, 
    matrixSize}, {{yMin, yMax}, {xMin, xMax}}];
myTable = Table[
   myCountFun[pointData[[i, j]]],
   {i, 1, matrixSize}, {j, 1, matrixSize}];
points1 = {};
i = 0;
While[i < matrixSize,
  i++;
  loc = FirstPosition[myTable[[i]], {{_, _}, 3}] // First;
  If[loc != "NotFound",
   AppendTo[points1, myTable[[i, loc, 1]]];
   ];
  ];
points2 = {};
i = 0;
While[i < matrixSize,
  i++;
  loc = FirstPosition[myTable[[i]], {{_, _}, 3}] // First;
  If[loc != "NotFound",
   thePoint = {-myTable[[i, loc, 1, 1]], myTable[[i, loc, 1, 2]]};
   AppendTo[points2, thePoint];
   ];
  ];
comboPoints = Join[points1, Reverse@points2];
    regionRange = Range[Length@comboPoints];
    myBoundary = 
      BoundaryMeshRegion[comboPoints, Line@(AppendTo[regionRange, 1])];
colorDataRules = {1 -> RGBColor[
   0.7725490196078432, 0.8509803921568627, 0.7333333333333333], 
  2 -> RGBColor[
   0.8117647058823529, 0.22745098039215686`, 0.8549019607843137], 
  3 -> RGBColor[
   0.3764705882352941, 0.6313725490196078, 0.5843137254901961], 
  4 -> RGBColor[
   0.7294117647058823, 0.7137254901960784, 0.7058823529411764], 
  5 -> RGBColor[
   0.8431372549019608, 0.5882352941176471, 0.3411764705882353]}
starRootMap = 
 ArrayPlot[myTable[[All, All, 2]], ColorRules -> colorDataRules, 
  FrameTicks -> Automatic, DataRange -> plotRange6, ImageSize -> 800, 
  Axes -> True, AxesStyle -> White, PlotLabel -> Style["Root Map", 16]]
comboPoints = Join[points1, Reverse@points2];
regionRange = Range[Length@comboPoints];
myBoundary = 
  BoundaryMeshRegion[comboPoints, Line@(AppendTo[regionRange, 1])];
boundaryMesh = 
 Show[myBoundary, Axes -> True, 
  PlotLabel -> Style["BoundaryMeshRegion", 16]]

I next shade the star with a chrome color with Specularity of 1 and place it in a black background with a few extra far-distant stars:

myColor = RGBColor[0.6667, 0.6627, 0.6784]
p1C = Plot3D[mySumFun[s, t], {s, -1, 1}, {t, -1, 1}, 
  PlotRange -> {{-1, 1}, {-1, 1}, {0, 9}}, ClippingStyle -> None, 
  BoxRatios -> {1, 1, 1}, 
  RegionFunction -> 
   Function[{s, t}, RegionMember[myBoundary, {s, t}]], 
  PlotStyle -> {myColor, Specularity[White, 1]}, Mesh -> None, 
  Background -> Black, Boxed -> False, Axes -> None
  ]
backGroundStars = 
  Table[Graphics3D@{Specularity[White, 5], myColor, 
     Sphere[{RandomReal[{-5, 5}], RandomReal[{-5, 5}], 
       RandomReal[{10, 15}]}, RandomReal[{0.001, 0.02}]]}
   , {250}];
Show[{p1C, backGroundStars}, PlotRange -> {{-5, 5}, {-5, 5}, {0, 20}}]

Update 2:

Note: I too realize this is not what the OP desired but still I learned a lot working on this. As my final version, I learned I could compose my star with another picture, preferably the dunes of Bethlehem. So that I next found a stock picture and composed my star of what I thought would be an idealized version for the three wise men to follow using the following code:

ImageCompose[image1, star, {920, 300}]

Note I can offset the star anywhere on the background with the third parameter so did so in the direction of the sunrise:

Answer 10

Not a complete solution but something to share nonetheless. Enjoy.

Rather than trying to invert the nasty polynomials, using ContourPlot functionality to generate points.

More work required but feel free to use. I am not in search of the bounty.

F[x_, y_] := (33 + 24 x^4 - 116 y^2 + 96 y^4 - 78 x^2 + 96 x^2 y^2)/25
Z[s_, t_] := Module[{},
  ptsmap = 
   ContourPlot[{x (-3 + 2 x^2 + 4 y^2)/5 == s, 
     y (-11 + 4 x^2 + 8 y^2)/5 == t}, {x, -50, 50}, {y, -50, 50}];
  If[Length[ptsmap[[1, 1]]] != 0, 
   Total[(1/Abs[F[#[[1]], #[[2]]]]) & /@ ptsmap[[1, 1, 1]]], 0]
  ]
dd = Table[Z[ind1, ind2], {ind1, -1, 1, .1}, {ind2, -1, 1, .1}];
ListDensityPlot[dd]

Answer 11

Some combination of code by Roman and xzczd gives plot same quality as second plot from xzczd answer but 50-100 times faster (on my ASUS laptop without compilation)

r[x_, y_] := {x (-3 + 2 x^2 + 4 y^2)/5, y (-11 + 4 x^2 + 8 y^2)/5};
fr[n_] := 
  Module[{nn = n}, 
   a = Transpose@
     BinCounts[
      Transpose[r @@ Transpose[RandomReal[{-2, 2}, {nn, 2}]]], {-1.01,
        1.01, .005}, {-1.01, 1.01, .005}]; a];
With[{array = fr[10^7]}, max = Max[array];
  ArrayPlot[array, ColorFunction -> "AvocadoColors", Frame -> None, 
   PlotRange -> {0, max/2}]] 

The nice property of function fr is that computation time linear depends on n.