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:
- Make a function which inverts the coordinates and gradients around a defined point (x,y)
- Fit a conic section to inverted data
- Plot the distance to the nearest focus as a function of (x,y)
- 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]
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]
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]
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:
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:
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:
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
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:
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.
Fit? $\endgroup$Interpolation$\endgroup$