Este grupo redactor ha estado mucho tiempo buscando la respuesta a tu interrogante, te ofrecemos la respuestas así que esperamos que resulte de gran apoyo.
Solución:
Aquí hay otro enfoque. Se podría mejorar (estoy seguro) para determinar correctamente los ejes principales y la traducción (si tengo tiempo, intentaré actualizar):
lin = #1^2, #1, #2, 2 #1 #2, #2^2 & @@@ points;
lm = LinearModelFit[lin, 1, a, b, c, d, a, b, c, d]
Explorando el modelo:
lm["ParameterTable"]
Determinación de la fórmula cuadrática:
pa = lm["BestFitParameters"];
w[x_, y_] := pa.1, x^2, x, y, 2 x y - y^2;
ContourPlot[w[x, y] == 0, x, -400, 400, y, -200, 500,
Epilog -> Point[points]]
Fórmula:
TraditionalForm[-w[x, y] == 0]
ACTUALIZAR
Publico esta actualización para determinar la traducción desde el origen, los ejes principales y $ a $ y $ b $ de la forma deseada $ frac x ^ 2 a ^ 2 + frac y ^ 2 b ^ 2 = 1 $ y, por tanto, evaluación del área de la elipse: a saber, $ pi ab $, etc. Los ejes se pueden intercambiar como se desee.
La traducción se puede derivar de la siguiente manera:
dx, dy =
0.5 Inverse[-pa[[2]], -pa[[5]], -pa[[5]], 1].pa[[3]], pa[[4]]
flexible:
0.0000180983, 221.
Los ejes se pueden derivar de la matriz de forma cuadrática:
mat = -pa[[2]], -pa[[5]], -pa[[5]], 1;
const = (-pa[[2]] dx^2 + dy^2 - 2 pa[[5]] dx dy) - pa[[1]]
trf[x_, y_] := -pa[[2]] (x - dx)^2 -
2 pa[[5]] (x - dx) (y - dy) + (y - dy)^2 - const
Luego, traduciendo la elipse al origen:
tran = [email protected][x + dx, y + dy]
rinde:
-49550.5 + 1.02504 x ^ 2 + 0.498521 xy + y ^ 2 = 0
Luego, mirar el sistema propio de la matriz permitirá la visualización de los ejes principales:
fun[poly_, a_] := Module[mat, val, vc, vcn, gr,
mat = #1, #2/2, #2/2, #3 & @@ (Coefficient[
poly, x^2, x y, y^2]);
val, vc = Eigensystem[mat];
vcn = Normalize /@ vc;
gr = Graphics[Line[0, 0, Sqrt[a] vcn[[1]]/Sqrt[[email protected][[1]]]],
Line[0, 0, Sqrt[a] vcn[[2]]/Sqrt[[email protected][[2]]]]];
Show[ContourPlot[poly == a, x, -500, 500, y, -500, 500], gr]
];
Los valores de $ a $ y $ b $ pueden derivarse:
ev, vec = Eigensystem[mat];
st = x^2, y^2.ev;
a, b = Sqrt[1/Coefficient[Expand[st/const], x^2, y^2]]
rinde:
198.143, 254.846
Poniendo todas las visualizaciones juntas:
cntr = ContourPlot[trf[x, y] == 0, x, -500, 500, y, -200, 500,
Epilog -> Point[points], Red, PointSize[0.03], Point[dx, dy]];
axes = fun[1.0250354570455185` x^2 + 0.49852055996040495` x y + y^2,
49550.46446156194`];
norm = ContourPlot[st == const, x, -500, 500, y, -500, 500];
Estos gráficos y las ecuaciones correspondientes se enlazaron y luego se exportaron como un gif animado:
grap = Show[#, PlotRange -> -500, 500, -500, 500] & /@ cntr,
axes, norm
eqns = Style[#, 20] & /@ TraditionalForm[trf[x, y] == 0],
TraditionalForm[
1.0250354570455185` x^2 + 0.49852055996040495` x y + y^2 ==
49550.46446156194`], TraditionalForm[ell[x, y] == 0]
El gif recorre (i) los datos con ajuste y el punto rojo es dx, dy (ii) el segundo es la elipse trasladada al origen (iii) los ejes de la elipse están “alineados” con los ejes cartesianos.
A continuación, se muestra una forma directa estándar de obtener los principales exes y otros datos de transformación. Encuentre la media de los puntos, reste para centrarlos y tome la descomposición del valor singular. El tercer y segundo componentes del mismo proporcionan los datos de rotación y escala necesarios para formar un círculo en el que se encuentra aproximadamente el primer componente, visto como un conjunto de puntos. La distancia media desde el origen de este nuevo conjunto de puntos da un radio. Ahora invierta todo eso para obtener la ecuación de elipse.
Si bien estoy seguro de que eso no tiene sentido (y menos para mí), la codificación es sencilla.
mean = Mean[points];
newpts = Map[# - mean &, points];
uu, ww, vv = SingularValueDecomposition[newpts, 2]
(* Out[195]= 0.154787303139,
0.24681671674, 0.264618646868, -0.0802782275335,
-0.244940409458,
0.132178292936, -0.2488804106, -0.141810231081,
-0.21297321559, -0.198934383898, -0.286761033433,
-0.0175688333828, -0.138307266939,
0.244940128897, 0.171124630754, -0.250504589808,
0.0288583749054, -0.287472943299, -0.0558601605337,
-0.280401389223, 0.186667941,
0.221277578296, 0.263211882778, -0.148978381541,
-0.154787157436, -0.246816544024, -0.264618872812,
0.0802779596975, 0.244940436909, -0.132178260395,
0.248880505624, 0.141810343723, 0.21297315013,
0.1989343063, 0.286760891953,
0.0175686656724, 0.13830729439, -0.244940096356,
-0.171124603303, 0.250504622349, -0.0288583474542,
0.28747297584, 0.0558601879848,
0.280401421764, -0.186667913548, -0.221277545756,
-0.263211855327, 0.148978414082, 876.375887309, 0., 0.,
671.508170856, -0.672351542658,
0.740231992746, 0.740231992746, 0.672351542658 *)
Aquí está el círculo, centrado en el origen, que los nuevos puntos (del factor de la izquierda, uu), ocupar.
ListPlot[uu, AspectRatio -> Automatic]
Ahora obtenga el radio al cuadrado de este nuevo círculo.
rsqr = Mean[Map[#.# &, uu]]
(* Out[191]= 0.0833333333333 *)
Para llegar a este círculo centrado en el origen, ejecutamos las transformaciones en nuestras variables, x, y, para obtener nuevas coordenadas en términos de las mismas.
nx, ny = Inverse[vv.ww].(x, y - mean);
Entonces aquí está la ecuación de la elipse.
expr = Expand[nx^2 + ny^2] == rsqr
(* Out[201]=
0.0838086622042 - 0.000201425757605 x + 1.80374774713*10^-6 x^2 -
0.00075844948748 y + 9.11428834461*10^-7 x y +
1.71594919287*10^-6 y^2 == 0.0833333333333 *)
Comprobaremos esto pictóricamente.
ContourPlot[[email protected], x, -400, 400, y, -200, 500,
Epilog -> Red, PointSize[Medium], Point[points]]
En geometría analítica, la elipse se define como el conjunto de puntos (X,Y) del plano cartesiano que, en casos no degenerados, satisfacen la ecuación implícita
con 
y 
dónde
Permite ajustar puntos con curvas de segundo orden (que incluyen elipse).
elipse = a11*x^2 + a22*y^2 + 2*a12*x*y + 2*a13*x + 2*a23*y + a33;
coeff = a11, a22, a12, a13, a23, a33;
fitResult=FindFit[points/.x_,y_:>x,y,0,elipse,#<0&/@coeff,coeff,x,y];
elipse=elipse/.fitResult
(*a11->0.002852802341,a22->0.002611181564,a12->0.0008537122851,a13->-0.188169993,a23->-0.5733778586,a33->-2.326357837*)
Luego verifica el resultado
Show[ContourPlot[res==0,x,Min[points],Max[points],y,Min[points],Max[points]],ListPlot[points],PlotRange->All]
