Ecuación de la esfera dada 4 puntos

Question

Halla la ecuación de la Esfera dándole los 4 puntos (3,2,1), (1,-2,-3), (2,1,3) y (-1,1,2).

La solución *fracasada* que probé es bastante directa:

Tenemos que encontrar el centro de la esfera.

Teniendo los puntos:

$$p_{1}(3,2,1),\, p_{2}(1,-2,-3),\, p_{3}(2,1,3),\, p_{4}(-1,1,2)$$

2 Se pueden crear triángulos utilizando estos puntos, llamémoslos $A$ nuestro triángulo formado por los puntos $p_{1},\,p_{2}\, and\,\, p_{3}$ Y $B$ nuestro triángulo formado por los puntos $p_{1},\, p_{3}\, and \,\,p_{4}$ .

Calcula los centroides de cada triángulo: $$CA = (2,1/3,1/3)\\ CB = (4/3,4/3,2)$$ Y también, un vector normal para cada triángulo: $$\overrightarrow{NA} = \overrightarrow{p_{1}p_{2}} \times \overrightarrow{p_{1}p_{3}}\\ \overrightarrow{NB} = \overrightarrow{p_{1}p_{3}} \times \overrightarrow{p_{1}p_{4}}$$

$$\overrightarrow{p_{1}p_{2}} = <-2,-4,-4>\\ \overrightarrow{p_{1}p_{3}} = <-1,-1,2>\\ \overrightarrow{p_{1}p_{4}} = <-4,-1,1>\\ \:\\ \overrightarrow{NA} = <-12, 8, -2>\\ \overrightarrow{NB} = <1, -7, -3>\\ $$

Con los centroides y las normales de los triángulos $A$ y $B$ podemos construir dos ecuaciones paramétricas para una línea, la primera interseca el centroide de $A$ y el otro el centroide $B$ . $$Line \enspace A\\ x = 2-12t \quad y = 1/3+8t \quad z = 1/3-2t\\ \:\\ Line \enspace B\\ x = 4/3 + s \quad y = 4/3 - 7s \quad z = 2 - 3s$$

El punto de intersección de estas líneas debería ser el centro de la esfera, desgraciadamente este sistema de ecuaciones no es linealmente dependiente, es decir, no se intersecan entre sí. ¿Cuál podría ser el problema aquí?

Answer 1

Citaría el hermoso método de W.H.Beyer para encontrar el centro y el radio de la esfera $(x-a)^2+(y-b^2)+(y-c)^2=R^2$

Answer 2

Otro método es comenzar con la ecuación de la esfera:

$$(x-u)^2+(y-v)^2+(z-w)^2=r^2$$

$(u,v,w)$ son las coordenadas del centro de la esfera y $r$ es el radio. Enchufa en los puntos dados $p_1, p_2, p_3, p_4$ para $x,y,z$ en la ecuación y se obtienen cuatro ecuaciones con variables $u,v,w,r$ . Pero estas ecuaciones contienen términos cuadráticos. Resta una ecuación de cada una de las otras tres. Se obtienen tres ecuaciones lineales y las variables $u, v, w$ . Este sistema puede resolverse de la forma habitual para encontrar el centro $(u,v,w)$ Por ejemplo, la eliminación gaussiana. Ahora, introduzca el valor calculado $u,v,w$ en una de las ecuaciones iniciales para encontrar el radio $r$ .

Nota:

Las tres ecuaciones lineales encontradas son la ecuación de los planos descritos por @Jyrki Lahtonen

El problema se puede resolver con este programa Maxima:

/\* define the four points p\[1\],...,p\[4\] \*/
p\[1\]:\[3,2,1\]; p\[2\]:\[1,-2,-3\]; p\[3\]:\[2,1,3\]; p\[4\]:\[-1,1,2\];

/\* ceq is the equation of the circle
(u,v,w) is the center, r is the radius \*/
ceq:(x-u)^2+(y-v)^2+(z-w)^2=r^2;

/\* plug in the points in the circle equation 
to get 4 equation eq1\[1\],...,eq1\[4\] \*/
for i:1 thru 4 do 
 eq1\[i\]:ev(ceq,map("=",\[x,y,z\],p\[i\]));

/\* display this 4 equations eq1\*/
listarray(eq1);

/\* subtract  the fourth equation for each of the first 
three equation to get three linear 
equations eq2\[1\],eq2\[2\],eq2\[3\] in u,v,w \*/

for i:1 thru 3 do 
 eq2\[i\]:ev(eq1\[i\]-eq1\[4\],expand);

/\* display these three equations 
eq2\[1\],eq2\[2\],eq2\[3\] \*/
listarray(eq2);

/\* solve this system oflinear euations to get u,v,w \*/
ss:solve(listarray(eq2),\[u,v,w\]);

/\* plugin the solutions in the fourth 
equation eq1\[4\] and calculate r^2 \*/
solve(eq1\[4\],r^2),ss\[1\];

En realidad, el problema se puede resolver con un solo Maxima solve comando:

(%i1) display2d:false$
(%i2) solve(makelist(ev((z-w)^2+(y-v)^2+(x-u)^2=r^2,
map("=",\[x,y,z\],p)),p,\[\[3,2,1\],\[1,-2,-3\],\[2,1,3\],\[-1,1,2\]\]),
\[u,v,w,r\]);
(%o2) \[\[u = 24/19,v = -16/19,w = 4/19,r = -3\*sqrt(470)/19\],
       \[u = 24/19,v = -16/19,w = 4/19,r = 3\*sqrt(470)/19\]\]

Answer 3

El centro de un triángulo no suele estar a la misma distancia de los vértices. En el plano, el centro del círculo exterior (= la intersección de las normales que bisecan los lados) tiene esta propiedad.

Pero, como esto es en 3D también podrías encontrar el punto de intersección de los planos que tienen $\vec{p_1p_2}$ , $\vec{p_2p_3}$ y $\vec{p_3p_4}$ como sus normales, y pasan por los puntos medios de esos segmentos de línea. Los puntos del primer plano están a la misma distancia de $p_1$ y $p_2$ etc. Por lo tanto, la intersección de esos tres planos es lo que quieres.

Answer 4

Sugerencia

Utilizando un enfoque totalmente diferente, también puede analizar su problema utilizando lo que JJacquelin (participante de MSE) propuso en su libro

http://fr.scribd.com/doc/14819165/Regressions-coniques-quadriques-circulaire-spherique

En las páginas 17 y 18 se ofrece el enfoque completo para una regresión esférica. Es bastante sencillo y se reduce a un sistema lineal de cuatro ecuaciones para cuatro incógnitas de las que se deducen las coordenadas del centro y el radio de la esfera.

Como no quiero estropearlo, no daré más. He aplicado su método a tu caso y funciona espléndidamente.

Answer 5

Para los que calculan a mano, los números propensos a errores son:

el centro de la esfera está en $p=(24/19,-16/19,4/19)$

el radio al cuadrado es $r^2=4230/361$ .

Aquí está mi pequeño GNU Maxima script

---

display2d : false;

/*
* purpose:
*   given points x_k,
*   calculate circumsphere center p and squared radius r^2
**/

my_A(d,k,x) := block(

    [res,i,j],

    res : zeromatrix(d,d),
    for j : 0 thru k-1 do block(
        for i : 1 thru d do block(
            res[i,j+1] : x[j+1][i]-x[k+1][i]
        )
    ),
    for j : k+1 thru d do block(
        for i : 1 thru d do block(
            res[i,j] : x[j+1][i]-x[k+1][i]
        )
    ),

    return(res)

);

my_norm2(d,x) := apply("+",makelist(x[k]^2,k,1,d));
my_dist2(d,x,y) := apply("+",makelist((x[k]-y[k])^2,k,1,d));

x : [[3,2,1],[1,-2,-3],[2,1,3],[-1,1,2]];

d : length(x)-1;
A : makelist(my_A(d,k,x),k,0,d);
AinvT : makelist(transpose(invert(A[k+1])),k,0,d);
xnorm2 : makelist(my_norm2(d,x[k+1]),k,0,d);

z1 : transpose(makelist(1,k,1,d));
p : -1/2*apply("+",makelist(xnorm2[k+1]*(AinvT[k+1].z1),k,0,d));
r2 : makelist(my_dist2(d,transpose(p)[1],x[k+1]),k,0,d);