Finding astroids

An astroid is a particular mathematical curve: a hypocycloid with four cusps.

Given a lot of data points I'm trying to fit such a curve to the data. As an astroid is also an algebraic curve of grade 6 I used LinearModelFit:

points=Import["https://pastebin.com/raw/Bi3DPvNj","Data"];
lin ={ #1 + #2, #1^2 + #2^2, #1^3 + #2^3, #1^4 + #2^4, #1^5 + #2^5, #1^6 + #2^6, #1 #2, #1^2 #2 + #1 #2^2, #1^3 #2 + #1 #2^3, #1^4 #2 + #1 #2^4, #1^5 #2 + #1 #2^5, #1^2 #2^2, #1^3 #2^2 + #1^2 #2^3, #1^4 #2^2+ #1^2 #2^4, #1^3 #2^3} & @@@ points;
lm = LinearModelFit[lin,Prepend[Table[Subscript[a,i],{i,1,14}],1],Table[Subscript[a,i],{i,1,14}]];
Quiet@lm["ANOVATable"]


Which results in

w[x_, y_] :=lm["BestFitParameters"].{1, x+y,  x^2+y^2, x^3+ y^3, x^4+ y^4, x^5+ y^5, x^6+y^6, x y, x y^2+ x^2 y, x y^3+ x^3 y, x y^4+x^4 y, x y^5+x^5 y, x^2 y^2, x^2 y^3+x^3 y^2, x^2 y^4+ x^4 y^2}-x^3 y^3;
ContourPlot[w[x, y] == 0, {x, 0.8, 1.0}, {y, 0.8, 1.0}, Prolog -> Point[points],ContourStyle->{Red},PlotPoints->109,ImageSize->Large]


Actually, the fit is okay, but ContourPlot seems to struggle with this polynomial as can be seen near 2 of the cusps (green). Is there a work-around? (Using many PlotPoints works, but takes very long.)

• Looks like it got clipped by ContourPlot, rather than it being an issue with the fit. Could be both, though. – Michael E2 Nov 25 '20 at 19:24
• With PlotPoints -> 1000, MaxRecursion -> 4 (~7 min): i.stack.imgur.com/NIfuZ.png – Michael E2 Nov 25 '20 at 19:32
• @MichaelE2 Thanks a lot, I couldn't use so many plot points. So, we know the fit is okay. ContourPlot seems to have some issues with those curves. – fwgb Nov 25 '20 at 19:39
• ContourPlot works on the intermediate-value-theorem principle: it draws the curve between sample points when the function values lies on opposite sides of the contour value (f == 0 in this case). It takes a lot of refinement to get enough points inside and outside a cusp to draw a cusp accurately. – Michael E2 Nov 25 '20 at 20:02

An alternative method for visualization and faster than ContourPlot with PlotPoints -> 1000 for nearly the same high quality in this case:

Needs@"NDSolveFEM"
bmesh = ToBoundaryMesh[
ImplicitRegion[
w[x, y] < 0 && 0.8 < x < 1.0 && 0.8 < y < 1.0, {x, y}],
"MaxBoundaryCellMeasure" -> 0.01];
Show[
Graphics[{Point[points]}],
bmesh["Wireframe"["MeshElementStyle" -> Red]],
Frame -> True,
PlotRange -> {{0.8, 1.0}, {0.8, 1.0}}
]


Another way, using NDSolve, which can get all the way into the cusps:

icy = 0.9;
icx = x /. First@NSolve[w[x, icy] == 0, x, Reals];
param = NDSolveValue[{
D[{x[t], y[t]}, t] == Cross@D[w[x[t], y[t]], {{x[t], y[t]}}],
{x[0], y[0]} == {icx, icy},
WhenEvent[y[t] == icy && Abs[x[t] - icx] < 10^-3,
"StopIntegration"]
}, {x, y}, {t, 0, Infinity}];
Graphics[{
Black, Point@points,
Red, Line@Transpose@Through[param["ValuesOnGrid"]]},
Frame -> True]


• Very good ideas! The curve seems to fit very well. (The points are created by applying a mapping to an ellipse.) – fwgb Nov 25 '20 at 20:50