       Re: Sphere formula

• To: mathgroup at smc.vnet.net
• Subject: [mg109516] Re: Sphere formula
• From: "Alexander Elkins" <alexander_elkins at hotmail.com>
• Date: Mon, 3 May 2010 06:10:58 -0400 (EDT)
• References: <hr65rj\$jq6\$1@smc.vnet.net> <hrbach\$hq9\$1@smc.vnet.net> <hrjh1b\$bh1\$1@smc.vnet.net>

```"Ray Koopman" <koopman at sfu.ca> wrote in message
news:hrjh1b\$bh1\$1 at smc.vnet.net...
> On May 1, 3:51 am, "Alexander Elkins" <alexander_elk... at hotmail.com>
> wrote:
> > "S. B. Gray" <stev... at ROADRUNNER.COM> wrote in
messagenews:hrbach\$hq9\$1 at smc.vnet.net...
> >> On 4/27/2010 1:05 AM, S. B. Gray wrote:
> >>> 1. The center of a sphere through 4 points has a very nice determinant
> >>> form. (http://mathworld.wolfram.com/Sphere.html) What I want is a nice
> >>> formula for the center of a sphere through 3 points, where the center
is
> >>> in the plane of the three points. I have a formula but it's a horrible
> >>> mess of hundreds of lines, even after FullSimplify.
> >>>
> >>> 2. (Unlikely) Is there a way to get Mathematica to put a long formula
> >>> into a matrix/determinant form if there is a nice one?
> >>>
> >>> Any tips will be appreciated.
> >>>
> >>> Steve Gray
> >>
> >> Thanks to everyone who answered my question, but there is a simpler
> >> answer. I forgot the simple fact that any linear combination of two
> >> vectors lies in the plane of the two vectors.
> >>
> >> Let the three points be p1,p2,p3. Consider the linear function
> >> p=b(p2-p1)+c(p3-p1) where b,c are to be determined and p is the desired
> >> center. Now do
> >>
> >> Solve[{Norm(p-p1)==Norm(p-p2),Norm(p-p1)==Norm(p-p3)},{b,c}].
> >>
> >> This gives b,c and therefore p, which will be equidistant from p1,p2,
> >> and p3 and lie in their plane. Very simple. (I used (p-p1).(p-p1) etc.
> >> instead of Norm.)
> >>
> >> Steve Gray
> >
> > Unless I misinterpreted something in this posting, the following does
not
> > give the expected center point {1, 2, 3} using Ray Koopman' s posted
values
> > for {p1, p2, p3} as the result:
> >
> > In:=With[{p1={0.018473,-1.1359,-0.768653},
> >              p2={2.51514,5.25315,6.48158},
> >              p3={0.818313,-0.881007,-1.0825}},
> >              With[{p=b(p2-p1)+c(p3-p1)},
> >                   p/.NSolve[{(p-p1).(p-p1)==(p-p2).(p-p2),
> >                              (p-p1).(p-p1)==(p-p3).(p-p3)},{b,c}]]]
> >
> > Out={{0.797494,2.34596,2.76487}}
> > [...]
>
> The fix is simple: just change the definition of p.
>
> Block[{p1 = {0.018473,-1.1359,-0.768653},
>        p2 = {2.51514,5.25315,6.48158},
>        p3 = {0.818313,-0.881007,-1.0825}, b,c,p},
>       p = (1-b-c)p1 + b*p2 + c*p3;
>       p /. Flatten@Solve[{(p-p1).(p-p1)==(p-p2).(p-p2),
>                           (p-p1).(p-p1)==(p-p3).(p-p3)},{b,c}]]
>
> {1.,2.,3.}
As shown, Solve complains and gives no result:

In:=Block[{p1={0.018473,-1.1359,-0.768653},
p2={2.51514,5.25315,6.48158},
p3={0.818313,-0.881007,-1.0825},b,c,p},
p=(1-b-c)p1+b*p2+c*p3;
p/.Flatten@Solve[{(p-p1).(p-p1)==(p-p2).(p-p2),
(p-p1).(p-p1)==(p-p3).(p-p3)},{b,c}]]

During evaluation of In:= Solve::tdep: The equations appear to \
involve the variables to be solved for in an essentially \
non-algebraic way. >>

During evaluation of In:= Solve::tdep: The equations appear to \
involve the variables to be solved for in an essentially \
non-algebraic way. >>

During evaluation of In:= ReplaceAll::reps: {Solve[<<1>>,{b,c}]} \
is neither a list of replacement rules nor a valid dispatch table, \
and so cannot be used for replacing. >>

During evaluation of In:= Solve::tdep: The equations appear to \
involve the variables to be solved for in an essentially \
non-algebraic way. >>

During evaluation of In:= General::stop: Further output of \
Solve::tdep will be suppressed during this calculation. >>

During evaluation of In:= ReplaceAll::reps: {Solve[<<1>>,{b,c}]} \
is neither a list of replacement rules nor a valid dispatch table, \
and so cannot be used for replacing. >>

Out= {2.51514 b + 0.018473 (1 - b - c) + 0.818313 c,
5.25315 b - 1.1359 (1 - b - c) - 0.881007 c,
6.48158 b - 0.768653 (1 - b - c) -
1.0825 c} /.Solve[{(0.768653+ 6.48158 b - 0.768653 (1 - b - c) -
1.0825 c)^2 + (1.1359+ 5.25315 b - 1.1359 (1 - b - c) -
0.881007 c)^2 + (-0.018473 + 2.51514 b +
0.018473 (1 - b - c) + 0.818313 c)^2 == (-6.48158 +
6.48158 b - 0.768653 (1 - b - c) - 1.0825 c)^2 + (-5.25315 +
5.25315 b - 1.1359 (1 - b - c) - 0.881007 c)^2 + (-2.51514 +
2.51514 b + 0.018473 (1 - b - c) + 0.818313 c)^2, (0.768653+
6.48158 b - 0.768653 (1 - b - c) - 1.0825 c)^2 + (1.1359+
5.25315 b - 1.1359 (1 - b - c) - 0.881007 c)^2 + (-0.018473 +
2.51514 b + 0.018473 (1 - b - c) + 0.818313 c)^2 == (1.0825+
6.48158 b - 0.768653 (1 - b - c) - 1.0825 c)^2 + (0.881007+
5.25315 b - 1.1359 (1 - b - c) - 0.881007 c)^2 + (-0.818313 +
2.51514 b + 0.018473 (1 - b - c) + 0.818313 c)^2}, {b, c}]

If Solve is changed to NSolve the correct result is returned:

In:=Block[{p1={0.018473,-1.1359,-0.768653},
p2={2.51514,5.25315,6.48158},
p3={0.818313,-0.881007,-1.0825},b,c,p},
p=(1-b-c)p1+b*p2+c*p3;
p/.Flatten@NSolve[{(p-p1).(p-p1)==(p-p2).(p-p2),
(p-p1).(p-p1)==(p-p3).(p-p3)},{b,c}]]

Out={1.,2.,3.}

NSolve (but not Solve) also works using a random center for testing:

In:=Block[{pc=RandomReal[{-10,10},3],
r=5Orthogonalize@RandomReal[NormalDistribution[0,1],{3,3}],
b,c,p,p1,p2,p3},
{p1,p2,p3}=Table[{Cos@#,Sin@#,0}&@RandomReal[{-Pi,Pi}].r+pc,{3}];
p=(1-b-c)p1+b*p2+c*p3;
{pc,#,pc==#}&[
p/.Flatten@NSolve[{(p-p1).(p-p1)==(p-p2).(p-p2),(p-p1).(p-p1)==
(p-p3).(p-p3)},{b,c}]]]

Out={{-1.87203,5.40096,7.652},{-1.87203,5.40096,7.652},True}

If instead an equation constraining the center p to be coplanar with
{p1,p2,p3} is added, then Solve for p={x,y,z} returns the expected result:

In:=Block[{p={x,y,z},p1={0.018473,-1.1359,-0.768653},
p2={2.51514,5.25315,6.48158},
p3={0.818313,-0.881007,-1.0825}},
p/.Flatten@Solve[(p-p1).(p-p1)==(p-p2).(p-p2)==(p-p3).(p-p3)&&
Inverse[{p1,p2,p3}].{1,1,1}.p==1,p]]

Out={1.,2.,3.}

These are all equivalent implicit plane equations suitable for constraining
the center p to be coplanar with {p1,p2,p3}:

Cross[p1-p2,p1-p3].(p-p1)==0
Cross[p1-p2,p1-p3].p==Det@{p1,p2,p3}
Inverse[{p1,p2,p3}].{1,1,1}.p==1