14
$\begingroup$

Inspired by the fantastic answers here and here, I would like to ask a questions along similar lines.

I have 6 curve types

enter image description here

that follow these rules:

  • Type A: If point 1 is positive and point 2 is negative, the curve will not pass through the origin.
  • Type B: If point 1 is negative and point 2 is positive, the curve will pass through the origin twice.
  • Types C-F: If both points are positive or both points are negative, the curve will pass through the origin once.

What is known:

  • The coordinates of points 1 and 2.
  • The gradient of tangents at points 1 and 2.
  • The approximate curve type (limaçon).

What is not known:

  • The arc length.

Data for A-F:

(* A *) {{0.000564367, 0.690525}, {-0.000689501, -0.984192}, 3.03065, -1.95699}
(* B *) {{-0.000689501, -0.984192}, {0.000664785, 1.07289}, -1.95699, 1.82419}
(* C *) {{0.000179304, 1.61576}, {0.0000936314, 0.852042}, 1.15014, 3.52804}
(* D *) {{0.000116063, 0.431337}, {0.000443491, 1.70111}, 2.88997, 1.41834}
(* E *) {{0.0000347276, 0.190688}, {0.000190634, 1.06651}, -3.77228, -2.08792}
(* F *) {{-0.000432719, -1.90935}, {-0.000142565, -0.645011}, -1.36691, -1.927}

in format: {{point 1}, {point 2}, gradient of tangent @ point 1, gradient of tangent at point 2}

Is it possible to estimate a curve fit (and hence arc length & skew) with only the data given?

Notes

I apologise for the lack of clarity on this question. The limacon is only a guess really as to the approximate shape. I have approximated the curves with limacons, but only after much fiddling. Any other approach is welcomed.

I can try to rephrase a little as follows: The problem can also be thought of like this:

A circle is cut and the points are slid up and down a vertical axis as shown:

enter image description here

Imagine that the tangents and positions of the endpoints are controllable, but there are no other parameters other than the theoretical elasticity (and tendency to move back to a circular form) are known.

I posted this image on a question on physics SE and John Rennie gave a link to an interesting PDF here on a similar subject.

It is essentially modelling a physical phenomenon. I have been playing around with Bezier curves, and have been getting fairly close, and even closer to the shapes doing it intuitively in Illustrator. I have even tried modelling it with wire - which gets the closest to the curves shown. Apart from this though, I am not having much success. As said previously, any approaches are welcome. I am looking for any insights into the problem.

Piecewise limacon example:

Show[PolarPlot[{ θ^(1/12) 1.815 Re[E^( I ( θ + 0.075))] - 
0.86}, {θ, 0.98, 2 Pi - 1.22}], 
PolarPlot[{1.875 Re[E^( I ( θ - 0.015))] - 0.95}, {θ, -1.02, 1.05}]]

Apologies for using the same image on both sites, but I thought it might be allowed since I was asking something essentially different! ;)

Data for one f the curves might help:

{{-0.000432719, -1.90935}, {-0.0572673, -1.82719}, {-0.108982,
-1.74322}, {-0.15552, -1.65773}, {-0.196844, -1.57104}, {-0.232935,
-1.48345}, {-0.263789, -1.39526}, {-0.289422, -1.30677}, {-0.309866,
-1.21829}, {-0.325172, -1.1301}, {-0.335405, -1.04249}, {-0.340648,
-0.955758}, {-0.341, -0.870171}, {-0.336576, -0.786}, {-0.327504,
-0.703507}, {-0.313929, -0.622942}, {-0.296007, -0.544546},
{-0.273908, -0.468551}, {-0.247817, -0.395174}, {-0.217925,
-0.324621}, {-0.184439, -0.257086}, {-0.147573, -0.192748},
{-0.107551, -0.131773}, {-0.0646045, -0.0743111}, {-0.0189726,
-0.0204992}, {0.0290993, 0.0295421}, {0.0793607, 
0.0757076}, {0.131557, 0.117908}, {0.185429, 0.156069}, {0.240717, 
0.190133}, {0.297158, 0.22006}, {0.35449, 0.245823}, {0.412451, 
0.267413}, {0.470781, 0.284837}, {0.529221, 0.298116}, {0.587518, 
0.307289}, {0.64542, 0.312408}, {0.702684, 0.313539}, {0.75907, 
0.310766}, {0.814345, 0.304183}, {0.868287, 0.293899}, {0.920678, 
0.280037}, {0.971313, 0.262729}, {1.02, 0.242122}, {1.06654, 
0.218372}, {1.11077, 0.191646}, {1.15252, 0.162119}, {1.19164, 
0.129975}, {1.22799, 0.0954078}, {1.26146, 0.0586159}, {1.29191, 
0.0198048}, {1.31927, -0.0208148}, {1.34344, -0.063028}, {1.36435,
-0.106616}, {1.38195, -0.151359}, {1.39619, -0.197032}, {1.40705,
-0.243414}, {1.41451, -0.29028}, {1.41858, -0.337408}, {1.41927,
-0.384578}, {1.41661, -0.431572}, {1.41063, -0.478177}, {1.4014,
-0.524182}, {1.38899, -0.569383}, {1.37347, -0.613581}, {1.35493,
-0.656586}, {1.33349, -0.698213}, {1.30925, -0.738286}, {1.28234,
-0.776638}, {1.25289, -0.81311}, {1.22105, -0.847555}, {1.18697,
-0.879835}, {1.15081, -0.909823}, {1.11274, -0.937404}, {1.07292,
-0.962474}, {1.03154, -0.984942}, {0.988782, -1.00473}, {0.944824,
-1.02176}, {0.899858, -1.03599}, {0.854076, -1.04738}, {0.80767,
-1.05589}, {0.760831, -1.06151}, {0.713753, -1.06423}, {0.666624,
-1.06407}, {0.619634, -1.06104}, {0.572969, -1.05518}, {0.526811,
-1.04654}, {0.481336, -1.03516}, {0.436719, -1.02113}, {0.393124,
-1.00451}, {0.350714, -0.985394}, {0.309639, -0.963889}, {0.270047,
-0.940099}, {0.232073, -0.914142}, {0.195846, -0.886144}, {0.161485,
-0.856239}, {0.129097, -0.824567}, {0.0987833, -0.791274},
{0.0706306, -0.756513}, {0.0447164, -0.720441}, {0.0211071,
-0.683219}, {-0.000142565, -0.645011}}

Update

Here is a start, using @ybeltukov 's spline method from here.

Manipulate[ pts = {{1, -1}, {0, ya}, {1, 1}, {2, 0}, {1, -1}, {0, yb}, {1, 1}};
dist = Norm /@ Differences[pts];
k = 0.5;
coeff = Rest[dist]/(Most[dist] + Rest[dist]);
diff = Differences[pts, 1, 2];
a1 = Join[pts - Flatten[{{{0, 0}}, k (1 - coeff) diff, {{0, 0}}}, 1]];
a2 = Join[pts + Flatten[{{{0, 0}}, k coeff diff, {{0, 0}}}, 1]];
bc = Drop[Drop[Flatten[{a1, pts, a2}, {{2, 1}}][[2 ;; -2]], 3], -3];
bezier = BezierCurve@bc;
Show[Graphics[{Blue, Thick, bezier, Black, Thick, 
Line[{{0, -2}, {0, 2}}]}, PlotRange -> {{-1, 3}, {-2, 2}}], 
ListPlot[{{0, ya}, {0, yb}}, PlotStyle -> {Red}]], 
{{ya, 0}, -1, 1}, {{yb, 0}, -1, 1}]

Of course the points are only acting independently at the moment. Since a circle can be approximated with only 4 cubic splines, a solution based on that should be fine. Though the other 3 points aren't known, I think using the gradient parameter of points 1 & 2, and the added parameter of having to (generally) pass through the origin should be enough.

Further update

Below are some images overlaying the curves in OP over @Alexey Bobrick 's excellent answer:

enter image description here

They approach the solution, but are missing some of the parameters:

  • The curves must pass through the origin if possible.
  • The curves should be clockwise from point 1 to point 2.
  • If a loop occurs on $y=0$, then the intersection of the loop should pass through the origin.

I believe that if these parameters are met, this approach may well provide a solution to the problem.

$\endgroup$
  • $\begingroup$ According to the wikipedia link, the polar form is just $r=b+a\cos\theta$ - so maybe you could convert your points to polar coordinates and use Fit? $\endgroup$ – Niki Estner Sep 29 '14 at 14:14
  • $\begingroup$ @nikie, have tried :/ $\endgroup$ – martin Sep 29 '14 at 14:25
  • 1
    $\begingroup$ it is really unclear what you are asking. What is the functional form of a "basically partial and scewed limacon"? $\endgroup$ – george2079 Oct 1 '14 at 16:56
  • 1
    $\begingroup$ If you are trying to capture the physics of the problem you need to work out the associated partial differential equation and then figure how to solve it. The physics forum is probably a better place for that If you just want to draw a picture you can use Interpolation $\endgroup$ – george2079 Oct 1 '14 at 21:18
  • 1
    $\begingroup$ Before starting to look for a solution: In least squares curve fitting distance deviations are minimized as per a definition. Likewise,apart from the given boundary conditions what defined distance is to be minimized here ? $\endgroup$ – Narasimham Oct 3 '14 at 21:03
11
+200
$\begingroup$

Major update, version 2.0

What changed: Is a different, more clever way of solving the problem. Overcomes most of the issues of previous versions.

Problem summary: I will rephrase the problem in simplest terms.

There exests a named curve called limacon (french, pronounced [ˈlɪməsɒn], means snail), described by a simple equation $r=a+b\cos{\theta}$ in polar coordinates.

In the problem this curve has been shifted and rotated somehow, and $a,b$ are unknowns too.

We know that the curve passes through 3 given points and has a known derivative at 2 of them. We need to recover $a,b,\theta_0,x_0,y_0$ from this data. In particular we have six datasets from the question author.

Why hard: The equation stops being as nice to work with when shifted and rotated. In fact, the equations on parameters converge poorly and are hard to solve (see older code versions).

Main idea: Might require a bit of thinking. Is a combination of the following points:

  • If we are at the origin of limacon and transform to $r\rightarrow\dfrac{1}{r}$, the curve equation becomes $r=\dfrac{p}{1+e\cos{\theta}}$, which is a known equation of a second order curve, as seen from its focus. But we don't know, where limacon's origin is.
  • If we apply inversion tranformation with respect to a random point (x,y), we can still always fit a second order curve to our transformed data points. See e.g. here for an excellent demonstration.
  • But if and only if (x,y) are at the origin of limacon, our fitted second order curve in the inverted space will have its focus at (x,y).

Hence, the algorithm:

  1. Make a function which inverts the coordinates and gradients around a defined point (x,y)
  2. Fit a conic section to inverted data
  3. Plot the distance to the nearest focus as a function of (x,y)
  4. For the points with zero distance return limacon parameters

Further, implementation (may be skipped):

This gets the data set:

(*Input data*)

Inputs = {(*A*){{0.000564367, 0.690525}, {-0.000689501, -0.984192}, 
    3.03065, -1.95699},
   (*B*){{-0.000689501, -0.984192}, {0.000664785, 1.07289}, -1.95699, 
    1.82419},
   (*C*){{0.000179304, 1.61576}, {0.0000936314, 0.852042}, 1.15014, 
    3.52804},
   (*D*){{0.000116063, 0.431337}, {0.000443491, 1.70111}, 2.88997, 
    1.41834},
   (*E*){{0.0000347276, 0.190688}, {0.000190634, 
     1.06651}, -3.77228, -2.08792},
   (*F*){{-0.000432719, -1.90935}, {-0.000142565, -0.645011}, \
-1.36691, -1.927}};
(*Auxiliary*)
Params[InputNum_] := Block[{Inp = Inputs[[InputNum]]},
   {X1 -> Inp[[1, 1]], X2 -> Inp[[2, 1]], Y1 -> Inp[[1, 2]], 
    Y2 -> Inp[[2, 2]], F1 -> Inp[[3]], F2 -> Inp[[4]], X3 -> 0, 
    Y3 -> 0}];

This inverts a general point, a general gradient and a particular selected dataset case around (x,y):

(*Inversion operations for points and gradients*)

PointInvert[X_, Y_, X0_, 
   Y0_] := {X0, Y0} + (#/(Norm[#])^2) &[{X - X0, Y - Y0}];
GradInvert[X_, Y_, F_, X0_, Y0_] := 
  Block[{DX, DY}, DX = X - X0; DY = Y - Y0;
   (2 DX DY + DY^2 F - DX^2 F)/(DX^2 - DY^2 + 2 DX DY F)];
(*For inverting all the data for given input number with respect to \
point (X0,Y0)*)

DataSetInvertedParams[InputNum_, X0_, Y0_] := 
 Block[{X1New, X2New, X3New, Y1New, Y2New, Y3New, F1New, F2New, 
   ParamsSet},
  ParamsSet = Params[InputNum];
  {X1New, Y1New} = PointInvert[X1, Y1, X0, Y0] /. ParamsSet;
  {X2New, Y2New} = PointInvert[X2, Y2, X0, Y0] /. ParamsSet;
  {X3New, Y3New} = PointInvert[X3, Y3, X0, Y0] /. ParamsSet;
  F1New = GradInvert[X1, Y1, F1, X0, Y0] /. ParamsSet;
  F2New = GradInvert[X2, Y2, F2, X0, Y0] /. ParamsSet;
  {X1 -> X1New, X2 -> X2New, Y1 -> Y1New, Y2 -> Y2New, F1 -> F1New, 
   F2 -> F2New, X3 -> X3New, Y3 -> Y3New}
  ]

This constructs the equations used for fitting conic sections for inverted data.

(*Fitting a conic section to inverted points. Type 1: Conic section \
goes through the point, Type 2: Conic section has a given derivative \
at a given point*)

InvertedEqType1[X1_, Y1_] := 
  1 + Binv X1 + Cinv Y1 + Dinv X1 Y1 + Einv X1^2 + Finv Y1^2 == 0;
InvertedEqType2[X1_, Y1_, F1_] := 
  Binv + Cinv F1 + Dinv (Y1 + X1 F1) + 2 Einv X1 + 2 Finv Y1 F1 == 0;

(*Assemble the explicit equations for a given input number*)

InvertedEqs[InputNumber_, X0_, Y0_] := 
  Block[{ParamSet, InvertedEqs1, InvertedEqs2, InvertedEqsAll},
   ParamSet = DataSetInvertedParams[InputNumber, X0, Y0];
   InvertedEqs1 = 
    And[InvertedEqType1[X1, Y1], InvertedEqType1[X2, Y2], 
      InvertedEqType1[X3, Y3]] /. ParamSet;
   InvertedEqs2 = 
    And[InvertedEqType2[X1, Y1, F1], InvertedEqType2[X2, Y2, F2]] /. 
     ParamSet;
   InvertedEqsAll = And[InvertedEqs1, InvertedEqs2]];

(*Auxiliary. Constructs the equations for conic section. May be used \
for visualising results*)

InvertedCurveEq[InputNumber_, X0_, Y0_] := 
  Block[{Eqs, Sol, X0Y0Eqs, X0New, Y0New},
   Eqs = InvertedEqs[InputNumber, X0, Y0];
   Sol = NSolve[Eqs, {Binv, Cinv, Dinv, Einv, Finv}];
   1 + Binv X + Cinv Y + Dinv X Y + Einv X^2 + Finv Y^2 == 0 /. Sol];

This is a big routine, which, for a given input and given (x,y) gets out all the canonical parameters of the second order curve. It is just simple geometry inside, which makes it big:

(*Fits the conic section, brings it to canonic form, extracts the \
distance to the nearest focus, and reconstructs the parameters of the \
limacon, constructed around the nearest focus. If the nearest focus \
coinsides with (X0,Y0), the limacon passes through the original \
dataset.*)
InvertedParams[InputNumber_?NumberQ, X0_?NumberQ, 
  Y0_?NumberQ] := 
 Block[{Eqs, QFormNew, Sol, X0Y0Eqs, X0New, Y0New, RotMatrix, 
   EigenVals, EigenVecs, RotationRules, Tmp, ANew, BNew, CNew, ENew, 
   FNew, CenterCurrPos, OriginalCenter, OffsetRules, FocalOffset, 
   FocalPointsOriginal, FocalPointLocal, DistToNearestFocus, 
   DistancesToFoci, FocusSign, Ecc, LCoef, aParamFin, bParamFin, 
   ThetaParamFin},
  (*Get the solution*)
  Eqs = InvertedEqs[InputNumber, X0, Y0];
  Sol = First[NSolve[Eqs, {Binv, Cinv, Dinv, Einv, Finv}]];
  (*Rotation transformation for the section to align it with the \
coordinate axes.*)
  {EigenVals, EigenVecs} = 
   Eigensystem[{{Einv, Dinv/2}, {Dinv/2, Finv}} /. Sol];
  (*Ensure that the rotation is such that the "biggest" axis is \
aligned with X axis.*)

  If[Not[Sort[#] == # &@
     Abs[EigenVals]], {EigenVals, 
      EigenVecs} = (Reverse[#] & /@ {EigenVals, EigenVecs}) Sign[
       Times @@ EigenVals];];
  (*Apply rotation to the quadratic form*)

  RotationRules = Thread[Rule[{X, Y}, Inverse[EigenVecs].{X, Y}]];
  QFormNew = (1 + Binv X + Cinv Y + Dinv X Y + Einv X^2 + Finv Y^2 /. 
       Sol) /. RotationRules // Simplify;
  (*Extract coefficients explicitly for detailed work*)
  {{ANew, 
     CNew, FNew}, {BNew, Tmp, Tmp}, {ENew, Tmp, Tmp}} = 
   CoefficientList[QFormNew, {X, Y}];
  (*Now get the position of the center of the quadratic form: 
  before and after the rotation*)

  CenterCurrPos = {-BNew/(2 ENew), -CNew/(2 FNew)};
  OriginalCenter = Inverse[EigenVecs].{-BNew/(2 ENew), -CNew/(2 FNew)};
  (*Reset the conic section so that it has its center at the origin*)

    OffsetRules = Thread[Rule[{X, Y}, {X, Y} + CenterCurrPos]];
  QFormNew = QFormNew /. OffsetRules // Simplify;
  (*Extract the remaining coefficients after the last transform*)
  \
{{ANew, Tmp, FNew}, {Tmp, Tmp, Tmp}, {ENew, Tmp, Tmp}} = 
   CoefficientList[QFormNew, {X, Y}];
  {ENew, FNew} = {ENew, FNew}/Abs[ANew];
  (*Determine geometric properties of the section*)
  (*How far are \
the foci from the center of the figure*)

  FocalOffset = Sqrt[Abs[1/ENew - 1/FNew]];
  (*Where are the foci located. In all cases, on X-axis, 
  except for if we are dealing with an inverted hyperbola.*)

  FocalPointLocal = 
   If[And[ENew > 0, ENew FNew < 0], {0, FocalOffset}, {FocalOffset, 
     0}];
  (*Where were the original foci were located*)

  FocalPointsOriginal = ({OriginalCenter + 
      Inverse[EigenVecs].FocalPointLocal, 
     OriginalCenter - Inverse[EigenVecs]. FocalPointLocal});
  (*Distances to each of them*)

  DistancesToFoci = Norm[# - {X0, Y0}] & /@ FocalPointsOriginal;
  (*This keep track, which of the two foci was the closest one*)

  FocusSign = If[DistancesToFoci == Sort[DistancesToFoci], 1, -1];
  (*Eccentricity of the curve, its polar parameter (auxiliary)*)

  Ecc = FocalOffset Sqrt[Abs[ENew]];
  LCoef = Sqrt[Abs[1./ENew]] (1. - Ecc^2);
  (*And finally, the original parameters of the limacon*)

  aParamFin = 1./LCoef;
  bParamFin = -(Ecc/LCoef) FocusSign;
  (*And the rotation angle*)

  ThetaParamFin = -ArcSin[EigenVecs[[1, 2]]];
  (*Assemble everything together*)
  {Min[DistancesToFoci], {X0, Y0, 
    ThetaParamFin, aParamFin, bParamFin}}]

Now the tasty part. For example, we can make the map of roots for a given data point (see below):

RootMap[InputNumber_, Range_] := 
 DensityPlot[
  InvertedParams[InputNumber, X0, Y0][[1]], {X0, -Range, 
   Range}, {Y0, -Range, Range}, PlotRange -> {0, 1}, 
  ColorFunction -> ColorData[{"DeepSeaColors", "Reverse"}], 
  PlotLegends -> Automatic, ImageSize -> 600, MaxRecursion -> 3, 
  ClippingStyle -> Automatic]

When we have found the interesting area, we can find a root with:

FindMin[InputNum_, {MinX0_, MaxX0_}, {MinY0_, MaxY0_}] := 
 Block[{DistanceFunction, Res, ResFin},
  DistanceFunction[X0_?NumberQ, Y0_?NumberQ] := 
   InvertedParams[InputNum, X0, Y0][[1]];
  Res = NMinimize[{(DistanceFunction[X0, Y0]), MinX0 < X0 < MaxX0, 
     MinY0 < Y0 < MaxY0}, {X0, Y0}, PrecisionGoal -> 5, 
    AccuracyGoal -> 5];
  ResFin = (InvertedParams[InputNum, X0, Y0] /. Res[[2]])[[2]];
  Print["Min value: ", Res[[1]]];
  Print["X0, Y0, \[Theta]0: ", ResFin[[1 ;; 3]]];
  Print["a, b: ", ResFin[[4 ;; 5]]];
  ResFin
  ]

If the root was successfully found, the function will print out the distance to the nearest root being zero (up to the data point precision).

And finally, we can see, how good was the fit with:

PointsShow[InputNumber_, Res_, PRange_] :=

 Block[{Curve, Curve2, Curve3, X0, Y0, a, b, t0},
  (*Get the parameters in*)
  {X0, Y0, t0, a, b} = Res;
  (*The conic section in inverted space. Might be used for debugging.*)

    Curve = 
   Block[{Eq, F}, Eq = InvertedCurveEq[InputNumber, X0, Y0][[1]];
    F[XC_, YC_] := Eq[[1]] /. {X -> XC, Y -> YC};
    ContourPlot[
     F[X1, Y1] == 0, {X1, -PRange, PRange}, {Y1, -PRange, PRange}, 
     PlotRange -> {{-PRange, PRange}, {-PRange, PRange}}, 
     Frame -> False, ImagePadding -> 0]];
  (*May be used for self-consistency check. 
  Backwards inversion of the fitted section. 
  Should be limacon if the center of inversion is in one of the foci*)
\
  (*Curve2=Block[{Eq,F,F2},Eq=InvertedCurveEq[InputNumber,X0,Y0][[
  1]];
  F[XC_,YC_]:=Eq[[1]]/.{X\[Rule]XC,Y\[Rule]YC};
  F2[XC_,YC_]:=Block[{XCNew,YCNew},{XCNew,YCNew}=PointInvert[XC,YC,X0,
  Y0];Eq[[1]]/.{X\[Rule]XCNew,Y\[Rule]YCNew}];
  ContourPlot[F2[X1,Y1]\[Equal]0,{X1,-PRange,PRange},{Y1,-PRange,
  PRange},PlotRange\[Rule]{{-PRange,PRange},{-PRange,PRange}},
  ContourStyle\[Rule]Red,Frame\[Rule]False,
  ImagePadding\[Rule]0]];*)
  (*Limacon, 
  constructed using the recovered parameters*)

  Curve3 = ParametricPlot[{X0 + (a + b Cos[t - t0]) Cos[t], 
     Y0 + (a + b Cos[t - t0]) Sin[t]}, {t, 0, 2 Pi}, 
    PlotRange -> {{-PRange, PRange}, {-PRange, PRange}}, 
    PlotStyle -> Red, Frame -> False, ImagePadding -> 0, 
    Axes -> False, PlotStyle -> Thick];
  Rasterize[
   Graphics[{{Inset[
        Show[Curve,(*Curve2,*)
         Curve3], {-PRange, -PRange}, {-PRange, -PRange}, 2 PRange], 
       PointSize -> Large, Point[{X1, Y1}], Point[{X2, Y2}], 
       Point[{X3, Y3}], Arrow[{{X1, Y1}, {X1 + 0.1, Y1 + 0.1 F1}}], 
       Arrow[{{X2, Y2}, {X2 + 0.1, Y2 + 0.1 F2}}]} /. 
      DataSetInvertedParams[InputNumber, X0, Y0], {Red, 
      PointSize -> 
       Large, {Point[{X1, Y1}], Point[{X2, Y2}], Point[{0, 0}], 
        Arrow[{{X1, Y1}, {X1 + 0.1, Y1 + 0.1 F1}}], 
        Arrow[{{X2, Y2}, {X2 + 0.1, Y2 + 0.1 F2}}]} /. 
       Params[InputNumber]}, Disk[{X0, Y0}, 0.05], Red, 
     Circle[{X0, Y0}, 1]}, 
    PlotRange -> {{-PRange, PRange}, {-PRange, PRange}}, 
    Frame -> True], ImageSize -> 800]]

Demonstration:

Now we can do many interesting things. Say, take system 1. Now we can easily visualise all the existing roots for it:

RootMap[1, 1]

Root map

We can see a lot from here. Brightest regions, where the distance to the focus is close to zero, indicate the possible limacon solutions. We can already expect that there is more than one solution.

For example:

Res11 = FindMin[1, {0.2, 0.5}, {-0.2, 0.2}]

reports that it has found a solution, which I can visualize:

PointsShow[1, Res11, 3]

Fit 11

We see the point of inversion, the fitted curve in the inverted space (here its ellipse) and the original limacon (red), neatly passing through the dataset.

But note: for the same dataset we can find another solution, as prompted by the root map:

Res12 = FindMin[1, {-0.2, -0.1}, {-0.2, 0.2}];
PointsShow[1, Res12, 3]

Fit12

This gives us another limacon passing through the dataset as neatly.

The problem now is easy to deal with.

Coda: To finish up I will comment on a few data points from the question and give a grid of the fits I got.

  • Data set 1 and a few more, hence, clearly have more than one fit. If the fit is to have its node at the origin, one needs to find a limacon with (x0,y0) close to (0,0). Whether such a root exists, one can tell from a root map.

  • Data set 4 has point positions, which do not coincide with the points shown on the plots in the original question.

  • Data set 5 seems to have no roots. Rootmaps with good recursion in a wide region indicate minimal distance to the focus of order 0.1. In the most promising region the solver tends to find close-to-degenerate solutions.

Finally, using:

Res11 = FindMin[1, {0.2, 0.5}, {-0.2, 0.2}];
Res2 = FindMin[2, {0.2, 0.5}, {-0.2, 0.2}];
Res3 = FindMin[3, {0., 0.3}, {0.8, 0.9}];
Res4 = FindMin[4, {0.1, 0.4}, {0.5, 0.7}];
Res6 = FindMin[6, {0.2, 0.3}, {-0.8, -0.75}];

One gets the following fits. To remind, set 5 does not have roots:

Final

What can be further improved: Among obvious points: 1) I didn't check what happens to degenerate conic sections, e.g. lines, and in other kinds of degenerate cases. 2) The code may be optimised to work faster. 3) And can be easily generalised to other types of input data.


Update, version 1.1:

Main point: The update is to see, what changes if we allow the curve to be rotated arbitrarily and impose an additional condition that it goes through the origin.

Main result: The modifications required are rather minimal. The method converges in some cases, but it gets much harder to specify good initial values.

Mathematical formulae are now, for type 1 equations: $$r^\prime-(x^\prime \cos{\theta_0} + y^\prime \sin{\theta_0})\dfrac{b}{r^\prime}=a$$ And for type 2 equations: $$b(\dfrac{y^\prime}{r^\prime}\cos{\theta_0}-\dfrac{x^\prime}{r^\prime}\sin{\theta_0})(\dfrac{y^\prime}{r^\prime}-f^\prime\dfrac{x^\prime}{r^\prime})=f^\prime y^\prime+x^\prime$$

Setting $\theta_0=0$ returns us back to the original equations.

Updated blocks:

Initialisation, as before:

(*Input data*)

Inputs = {(*A*){{0.000564367, 0.690525}, {-0.000689501, -0.984192}, 
    3.03065, -1.95699},
   (*B*){{-0.000689501, -0.984192}, {0.000664785, 1.07289}, -1.95699, 
    1.82419},
   (*C*){{0.000179304, 1.61576}, {0.0000936314, 0.852042}, 1.15014, 
    3.52804},
   (*D*){{0.000116063, 0.431337}, {0.000443491, 1.70111}, 2.88997, 
    1.41834},
   (*E*){{0.0000347276, 0.190688}, {0.000190634, 
     1.06651}, -3.77228, -2.08792},
   (*F*){{-0.000432719, -1.90935}, {-0.000142565, -0.645011}, \
-1.36691, -1.927}};
(*Auxiliary*)

Params[Inp_] := {X1 -> Inp[[1, 1]], X2 -> Inp[[2, 1]], 
   Y1 -> Inp[[1, 2]], Y2 -> Inp[[2, 2]], F1 -> Inp[[3]], 
   F2 -> Inp[[4]]};

New equations:

(*Two types of the equations and a function to simplify them.*)

EqSimplify[
   Eq_] := (Eq /. {R -> Sqrt[X^2 + Y^2]}) /. {X -> X - X0, 
     Y -> Y - Y0} // Simplify;
EqType1[XInp_, YInp_] := 
  Evaluate[EqSimplify[(R - (b/R) (X Cos[t0] + Y Sin[t0])) == 
      a]] /. {X -> XInp, Y -> YInp};
EqType2[XInp_, YInp_, FInp_] := 
  Evaluate[EqSimplify[
     b (Y Cos[t0]/R - X Sin[t0]/R) (Y/R - F X/R) == (F Y + 
        X)]] /. {X -> XInp, Y -> YInp, F -> FInp};

These functions apply a given equation type to a given data set.

Eqs1[InputNumber_] := 
  And[EqType1[X1, Y1], EqType1[X2, Y2], EqType1[0., 0.]] /. 
   Params[Inputs[[InputNumber]]];
Eqs2[InputNumber_] := 
  And[EqType2[X1, Y1, F1], EqType2[X2, Y2, F2]] /. 
   Params[Inputs[[InputNumber]]];
EqsAll[InputNumber_] := 
  And[Eqs1[InputNumber], Eqs2[InputNumber]] // Simplify;

And these two functions form the trial-and-error block:

FindSol[InputNumber_, GuessA_, GuessB_, GuessX0_, GuessY0_, 
  GuessT0_] := Block[{EqsFin},
  EqsFin = And[Eqs1[InputNumber], Eqs2[InputNumber]] // Simplify;
  FindRoot[
   EqsFin, {{a, GuessA}, {b, GuessB}, {X0, GuessX0}, {Y0, 
     GuessY0}, {t0, GuessT0}}, AccuracyGoal -> 10, 
   PrecisionGoal -> 10, MaxIterations -> 100]]

Test[InputNumber_, GuessA_, GuessB_, GuessX0_, GuessY0_, GuessT0_] := 
 Block[{Sol, Curve, CurveWas, PosNew},
  Sol = FindSol[InputNumber, GuessA, GuessB, GuessX0, GuessY0, 
    GuessT0];
  Print[Sol];
  PosNew = {X0, Y0} /. Sol;
  CurveWas = 
   Rasterize[
    ParametricPlot[{X0 + (a + b Cos[t - t0]) Cos[t], 
       Y0 + (a + b Cos[t - t0]) Sin[t]} /. {X0 -> GuessX0, 
       Y0 -> GuessY0, a -> GuessA, b -> GuessB, t0 -> GuessT0}, {t, 0,
       2 Pi}, PlotRange -> {{-2, 2}, {-2, 2}}, PlotStyle -> Red, 
     Frame -> False, ImagePadding -> 1, Axes -> False, 
     PlotStyle -> Thick], ImageSize -> 800];
  Curve = 
   Rasterize[
    ParametricPlot[{{X0 + (a + b Cos[t - t0]) Cos[t], 
        Y0 + (a + b Cos[t - t0]) Sin[t]} /. 
       Sol, {X0 + (a + b Cos[t - t0]) Cos[t], 
        Y0 + (a + b Cos[t - t0]) Sin[t]} /. {X0 -> GuessX0, 
        Y0 -> GuessY0, a -> GuessA, b -> GuessB, t0 -> GuessT0}}, {t, 
      0, 2 Pi}, PlotRange -> {{-2, 2}, {-2, 2}}, Frame -> False, 
     ImagePadding -> 1, Axes -> False, 
     PlotStyle -> {{Thick, Black}, {Red}}], ImageSize -> 800];
  Rasterize[
   Graphics[{Inset[Curve, {-2, -2}, {-2, 2}, 
      4], {PointSize -> Large, Point[{X1, Y1}], Point[{X2, Y2}], 
       Point[{0, 0}], Arrow[{{X1, Y1}, {X1 + 0.1, Y1 + 0.1 F1}}], 
       Arrow[{{X2, Y2}, {X2 + 0.1, Y2 + 0.1 F2}}]} /. 
      Params[Inputs[[InputNumber]]], Red, 
     Arrow[{{GuessX0, GuessY0}, PosNew}]}, 
    PlotRange -> {{-2, 2}, {-2, 2}}, Frame -> True], 
   ImageSize -> 400]]

The initial value for one more variable (curve rotation angle) needs to be specified. Additionally, it plots the originally guessed curve, which might help with chosing better starting values.

For the curves of type A, B it works very well.

For example,

GraphicsRow[{Test[1, 1.3, -1, 0.4, 0., Pi],Test[2, 1, 2, 0.01, 0.5, 0.01]}]

gives this:

Good fit example

Here red color denotes the original guess, and black is the fit. Arrow indicates the offset of the curve center.

But for some other data sets the fit is worse:

Test[3, -.5, 2.1, 1.7, -0.05, 2.8]

gives this, for example:

Example of a bad fit

Maybe it would converge eventually, or maybe the root doesn't exist for these points.


Version 1

Here is one possible solution.

Mathematical preamble:

The curve is known to be described by $r=a+b\cos{\theta}$ in polars. Two parameters $a,b$ are to be determined. We have four constraints, however. Therefore, the curve must have been shifted: $x^\prime=r\cos{\theta}, y^\prime=r\sin{\theta}$, $x=x_0+x^\prime$, $y=y_0+y^\prime$.

We have then two types of equations:

Type 1 states that the data point belongs to the curve. Rewriting the curve equation in cartesians, we get $(r^2 - b x)^2 = a^2 r^2$, with $r^2 = x^{\prime2} + y^{\prime2}$.

Type 2 equation states that at a given point $(x,y)$ one has $dy/dx=f^\prime$. In cartesians this statement is formulated as: $b y^\prime (y^\prime - f^\prime x^\prime) = r^2 (f^\prime y^\prime + x^\prime)$.

In Mathematica these equations are:

(*Two types of the equations and a function to simplify them.*)
EqSimplify[
   Eq_] := (Eq /. {R^2 -> X^2 + Y^2}) /. {X -> X - X0, Y -> Y - Y0} //
    Simplify;
EqType1[XInp_, YInp_] := 
  Evaluate[EqSimplify[(R^2 - b X)^2 == a^2 R^2]] /. {X -> XInp, 
    Y -> YInp};
EqType2[XInp_, YInp_, FInp_] := 
  Evaluate[EqSimplify[b Y (Y - F X) == R^2 (F Y + X)]] /. {X -> XInp, 
    Y -> YInp, F -> FInp};

We also want to initialise the input data:

(*Input data*)
Inputs = {(*A*){{0.000564367, 0.690525}, {-0.000689501, -0.984192}, 
    3.03065, -1.95699},
   (*B*){{-0.000689501, -0.984192}, {0.000664785, 1.07289}, -1.95699, 
    1.82419},
   (*C*){{0.000179304, 1.61576}, {0.0000936314, 0.852042}, 1.15014, 
    3.52804},
   (*D*){{0.000116063, 0.431337}, {0.000443491, 1.70111}, 2.88997, 
    1.41834},
   (*E*){{0.0000347276, 0.190688}, {0.000190634, 
     1.06651}, -3.77228, -2.08792},
   (*F*){{-0.000432719, -1.90935}, {-0.000142565, -0.645011}, \
-1.36691, -1.927}};
(*Auxiliary*) 
Params[Inp_] := {X1 -> Inp[[1, 1]], X2 -> Inp[[2, 1]], 
   Y1 -> Inp[[1, 2]], Y2 -> Inp[[2, 2]], F1 -> Inp[[3]], 
   F2 -> Inp[[4]]};

Then, we can formulate explicitly both types of equations for any given point:

Eqs1[InputNumber_] := 
  And[EqType1[X1, Y1], EqType1[X2, Y2]] /. 
   Params[Inputs[[InputNumber]]];
Eqs2[InputNumber_] := 
  And[EqType2[X1, Y1, F1], EqType2[X2, Y2, F2]] /. 
   Params[Inputs[[InputNumber]]];
EqsAll[InputNumber_] := 
  And[Eqs1[InputNumber], Eqs2[InputNumber]] // Simplify;

For example, equations on $a,b,x_0,y_0$ for the first point are:

   EqsAll[1]//TableForm

Output in Table form

One can see, these are four polynomial equations, two of order 3 and two of order 4. One may, therefore expect several solutions to them to exist.

The best way to deal with this that I could find is to use FindRoot. This requires, however, some guesswork on starting points.

FindSol[InputNumber_, GuessA_, GuessB_, GuessX0_, GuessY0_] := 
 Block[{EqsFin},
  EqsFin = And[Eqs1[InputNumber], Eqs2[InputNumber]] // Simplify;
  FindRoot[
   EqsFin, {{a, GuessA}, {b, GuessB}, {X0, GuessX0}, {Y0, GuessY0}}, 
   PrecisionGoal -> 5, AccuracyGoal -> 5]]@ 

This means, that each point has to be dealt with one-by-one, trying several initial guesses sometimes. In other words, it is hard to automatise this procedure.

This can be made interactive, however. This function will output the found solutions and plot them. If the plot is not good enough, the function can be called again:

Test[InputNumber_, GuessA_, GuessB_, GuessX0_, GuessY0_] := 
 Block[{Sol, Curve},
  Sol = FindSol[InputNumber, GuessA, GuessB, GuessX0, GuessY0];
  Print[Sol]; 
  Curve = Rasterize[
    ParametricPlot[{X0 + (a + b Cos[t]) Cos[t], 
       Y0 + (a + b Cos[t]) Sin[t]} /. Sol, {t, 0, 2 Pi}, 
     PlotRange -> {{-2, 2}, {-2, 2}}, Frame -> False, Axes -> False, 
     PlotStyle -> Thick], ImageSize -> 800];
  Rasterize[
   Graphics[{Inset[Curve, {-2, -2}, {-2, 2}, 
      4], {PointSize -> Large, Point[{X1, Y1}], Point[{X2, Y2}], 
       Arrow[{{X1, Y1}, {X1 + 0.1, Y1 + 0.1 F1}}], 
       Arrow[{{X2, Y2}, {X2 + 0.1, Y2 + 0.1 F2}}]} /. 
      Params[Inputs[[InputNumber]]]}, PlotRange -> {{-2, 2}, {-2, 2}},
     Frame -> True], ImageSize -> 400]]

For example, I have found these solutions:

GraphicsGrid[{{Test[1, 2, -2, 3, -0.3], 
   Test[2, 3, 3, 1.5, 1]}, {Test[3, -2, -1, -2, 0], 
   Test[4, 1, 3, -1, 2]},
  {Test[5, -4, -5, -2, 1], Test[6, 1, -2, 0, -2]}}]

Which produces the following:

Resulting images

The method thus allows one to manually find the solutions by trying initial conditions. The convergence gets a bit harder to achieve if one of the points ends up on bifurcation part of the curve. Also, more then one solution may exist for a given set of points.

One other comment, of course, is that the curve may be rotated by an angle $\theta_0$. To solve for this would require one more constraint, though.

$\endgroup$
  • $\begingroup$ Thanks for your answer, these are great! The curves have not only been shifted however, they have been skewed. Is there any way of introducing this into the solution? I can include data if that would help. $\endgroup$ – martin Oct 2 '14 at 18:01
  • $\begingroup$ A third point and its gradient could be given. Would this be enough? $\endgroup$ – martin Oct 2 '14 at 18:13
  • $\begingroup$ Aha! The other parameter is that the curve must pass through the origin if possible (with the exception of the first one - I don't know if there is a solution for this, but I doubt it). Please see update. $\endgroup$ – martin Oct 2 '14 at 20:01
  • $\begingroup$ @martin: Let's see! So, (a,b,x0,y0,\theta0) would be sought for, right? The conditions would be that the curve passes through the origin, points 1 and 2, and has known gradients at points 1, 2. $\endgroup$ – Alexey Bobrick Oct 3 '14 at 9:26
  • $\begingroup$ And then, two other constraints would allow one to select a solution, I suppose. $\endgroup$ – Alexey Bobrick Oct 3 '14 at 9:27

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.