10
$\begingroup$

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]

enter image description here

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.)

$\endgroup$
  • $\begingroup$ Looks like it got clipped by ContourPlot, rather than it being an issue with the fit. Could be both, though. $\endgroup$ – Michael E2 Nov 25 '20 at 19:24
  • 3
    $\begingroup$ With PlotPoints -> 1000, MaxRecursion -> 4 (~7 min): i.stack.imgur.com/NIfuZ.png $\endgroup$ – Michael E2 Nov 25 '20 at 19:32
  • $\begingroup$ @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. $\endgroup$ – fwgb Nov 25 '20 at 19:39
  • 2
    $\begingroup$ 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. $\endgroup$ – Michael E2 Nov 25 '20 at 20:02
11
$\begingroup$

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

Needs@"NDSolve`FEM`"
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}}
 ]

enter image description here

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]

enter image description here

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

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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