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[1]:=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[1]={{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[1]:=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[1]:= Solve::tdep: The equations appear to \ involve the variables to be solved for in an essentially \ non-algebraic way. >> During evaluation of In[1]:= Solve::tdep: The equations appear to \ involve the variables to be solved for in an essentially \ non-algebraic way. >> During evaluation of In[1]:= 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[1]:= Solve::tdep: The equations appear to \ involve the variables to be solved for in an essentially \ non-algebraic way. >> During evaluation of In[1]:= General::stop: Further output of \ Solve::tdep will be suppressed during this calculation. >> During evaluation of In[1]:= 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[1]= {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[2]:=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[2]={1.,2.,3.} NSolve (but not Solve) also works using a random center for testing: In[3]:=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[3]={{-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[4]:=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[4]={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 Det@PadRight[{p,p1,p2,p3},{4,4},1]==0 Det@Append[Transpose@{p,p1,p2,p3},{1,1,1,1}]==0 Alexander Elkins