# Sphere equation given 4 points

Find the equation of the Sphere give the 4 points (3,2,1), (1,-2,-3), (2,1,3) and (-1,1,2).

The *failed* solution I tried is kinda straigh forward:

We need to find the center of the sphere.

Having the points:

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

2 Triangles can be created using these points, let's call $A$ our triangle formed by the points $p_{1},\,p_{2}\, and\,\, p_{3}$; And $B$ our triangle formed by the points $p_{1},\, p_{3}\, and \,\,p_{4}$.

Calculate the centroids of each triangle: $$CA = (2,1/3,1/3)\\ CB = (4/3,4/3,2)$$ And also, a normal vector for each triangle: $$\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>\\$$

With the centroids and normals of triangles $A$ and $B$, we can build two parametric equations for a line, the first one intersects the centroid of $A$ and the other one the centroid $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$$

The point where this lines intersect should be the center of the sphere, unfortunately this system of equations is not linearly dependent, that means that they do not intersect each other. What could be the problem here?

-
+1 for showing your work and explaining your difficulty in detail. This is essential to giving helpful answers. Many newbies fail there. Last but not least. Welcome to Math.SE! –  Jyrki Lahtonen Aug 12 '14 at 7:04
@JyrkiLahtonen Thank you! I really like the stackexchange sites. –  DanielRS Aug 12 '14 at 10:42

I would cite the beautiful method from W.H.Beyer to find the center and radius of the sphere $(x-a)^2+(y-b^2)+(y-c)^2=R^2$

-
This is what I was refering to in my answer. –  Claude Leibovici Aug 12 '14 at 9:53
Hi Claude ! Thanks for your notification. Of course, the method from W.H.Bayer in specific to the case of four given points. In case of more than four scattered points one have to use a fitting method, which is something else. Cordially, JJ. –  JJacquelin Aug 12 '14 at 10:14
For me, the problem is the same but I bet that this could take us for a loooong discussion. Cheers. By the way, I really enjoy your books. –  Claude Leibovici Aug 12 '14 at 10:16
Thank you! After thinking about it a little, my approach was incorrect. –  DanielRS Aug 12 '14 at 10:38

The centroid of a triangle is usually not at the same distance from the vertices. In the plane the center of the outer circle (= the intersection of the normals bisecting the sides) has this property.

But, because this is in 3D you might as well find the point of intersection of the planes the have $\vec{p_1p_2}$, $\vec{p_2p_3}$ and $\vec{p_3p_4}$ as their normals, and pass through the midpoints of those line segments. The points on the first plane are at the same distance from $p_1$ and $p_2$ et cetera. Thus the intersection of those three planes is what you want.

-

Hint

Using a totally different approach, you can also look at your problem using what JJacquelin (participant to MSE) proposed in his book

Pages 17 and 18 give the full approach for a spherical regression. It is quite simple and reduces to a linear system of four equations for four unknowns from which are deduced the coordinates of the center and the radius of the sphere.

Since I do not want to spoil him, I shall not give more. I applied his method to your case and it works splendid.

-

For the ones calculating by hand, the error-prone numbers are:

the center of the sphere is at $p=(24/19,-16/19,4/19)$

the squared radius is $r^2=4230/361$.

Here is my small 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);

-
Essentialiy I'm constructing the center point $p$ as a linear combination of the outward normals to the simplex $S = Conv(x_0,x_1,x_2,x_3)$, where the $x_k$ are the given points. –  andre Aug 12 '14 at 11:41

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

$(u,v,w)$ are the coordinates of the center of the sphere and $r$ is the radius. Plug into the given points $p_1, p_2, p_3, p_4$ for $x,y,z$ in the equation and you get four equations with variables $u,v,w,r$. But these equations contain quadratic terms. Subtract one equation from each of the other three. You get three linear equations and the variables $u, v, w$. Such a system can be solved in the usual way to find the center $(u,v,w)$, e.g. Gaussian elimination. Now plug the calculated $u,v,w$ into one of the initial equations to find the radius $r$.

Note:

The three linear equations found are the equation of the planes described by @Jyrki Lahtonen

The problem can be solved with this Maxima program:

/* 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];



Actually the problem can be solved using only one Maxima solve command:

(%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]]

-