Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

I first make a function to get a random vector on unit sphere in a swath around the equator. That is what the parameter $\gamma$ controls; if $\gamma = 1/2$, the vectors can be chosen anywhere on the sphere; if $\gamma = 0$, the vectors can only be chosen from the equatorial plane etc.

randomVectorOnUnitSphere[\[Gamma]_] := Module[{},
\[Phi] = RandomReal[{0, 2 \[Pi]}];
w = RandomReal[{1/2 - \[Gamma], 1/2 + \[Gamma]}];
nx = Cos[\[Phi]] Sin[ArcCos[2 w - 1]];
ny = Sin[\[Phi]] Sin[ArcCos[2 w - 1]];
nz = Cos[ArcCos[2 w - 1]];
Return[{nx, ny, nz}];
]

These two random vectors span a plane. I would like to find an orthonormal basis for that plane so I use the Gram-Schmidt process as so:

randMod[\[Gamma]_] := Module[{},
v1 = randomVectorOnUnitSphere[\[Gamma]];(*random unit vector 1*)
v2 = randomVectorOnUnitSphere[\[Gamma]];(*random unit vector 2*)
u1 = v1; (* start Gram-Schmidt process *)
proj = u1.v2/u1.u1 u1;
u2temp = v2 - proj;
u2 = u2temp/Sqrt[u2temp.u2temp];
Return[{u1,u2}] (*the output here is 2 orthonormal vectors that span the plane*)
];

I look at 10000 pairs and plot on unit sphere:

randTable1 = ParallelTable[randMod[1/10][[1]], {i, 1, 10000}];
randTable2 = ParallelTable[randMod[1/10][[2]], {i, 1, 10000}];
randPlot1 = 
ListPointPlot3D[{randTable1, {{0, 0, 1}}}, 
PlotStyle -> {Black, {PointSize -> 0.025, Red}}];
randPlot2 = 
ListPointPlot3D[{randTable2, {{0, 0, 1}}}, 
PlotStyle -> {Black, {PointSize -> 0.025, Red}}];
{Show[Graphics3D[Sphere[{0, 0, 0}, 1]], randPlot1], 
Show[Graphics3D[Sphere[{0, 0, 0}, 1]], randPlot2]} 

and get the output for all the $u_1$ (top) and $u_2$ (bottom):

This is $u_1$This is $u_2$ (not sure why size larger, sorry)

The $u_2$ vectors are concentrated at the equator but there is a sizable number that are spreadout over the sphere. I would expect to get something that looked just like $u_1$ instead.

The initial vectors $v_1$ and $v_2$ lie in the same 'swath' as $u_1$ so I know that it is the G-S process that is failing.

I have looked at some of the culprit points and have noticed that their initial starting vectors were near parallel. I think this is the cause of the problem. I tried increasing WorkingPrecision when I call RandomReal but that did not seem to have any noticeable effect. Is there a way I can fix this problem?

share|improve this question
2  
Why don't you expect vectors u2 away from the equatorial plane? Suppose you start with v1={1,0,0} and v2={Sqrt[1-eps^2],0,eps} where eps is small, then u2 will be {0,0,1}. –  Heike Feb 3 '12 at 21:18
add comment

1 Answer

I'm not sure why you'd expect the resulting vectors to remain in the original swath because you're forcing the angle between the two resulting vectors to be 90 deg. Pairs of vectors with cross products far from the poles of the sphere will necessarily result in one of the vectors "popping out" of the swath.

I'd also like to point out that you don't need Gram-Schmidt to create an orthonormal basis for the two vectors. Instead, you can just do the following: Given two vectors v1 and v2, define u1=v1 then construct u2=Normalize[v1 x v2]. Finally, u3 = u2 x u1. Now u1 will be aligned with v1 and u3 will be in the plane defined by v1 and v2 but it will be orthogonal to u1 and in the general direction of v2.

share|improve this answer
4  
You could also construct a basis by using Orthogonalize[{v1,v2, Cross[v1,v2]}] –  Heike Feb 3 '12 at 21:11
    
ok, I see my error. I did not visualize how the basis vectors could get to the poles but now I see. For my situation I want the basis vectors to remain in the swath. Should I delete this question and ask a new one or modify this one? –  BeauGeste Feb 3 '12 at 23:10
1  
Generally you don't delete a question unless it's seriously terrible (and this is not). I'd just ask a new question. –  David Z Feb 4 '12 at 3:23
add comment

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.