From octave-sources-request at bevo dot che dot wisc dot edu Fri Jan 12 08:09:54 2001 Subject: cross.m ( correct in unitary space ) From: Rolf Fabian To: "'octave-sources UWISC'" Cc: "'melqart at free dot fr'" Date: Fri, 12 Jan 2001 15:06:21 +0100 ----------------------------------------------------------- Author: Rolf Fabian Jan. 12, 2001 ----------------------------------------------------------- This is a study on the vector (cross) product z = cross( x,y ) with complex vector arguments. In order to be correct, we request that the z is orthogonal to BOTH argument vectors and alternates sign if x,y are exchanged. The study clearly demonstrates, that Octave's V2.1.33 cross.m function returns incorrect results ( as well as the corresponding Matlab function ). A corrected cross.m -script may be found at the end of this upload. This study uses the improved/corrected 'dot.m' script uploaded yesterday Jan.11 (rev.2) for orthogonality checks. ----------------------------------------------------------- #We define two (row) test vectors, nothing special about them: :) x =[ 1 + 0i 1 - 1i 1 + 1i ] :) y =[ 3 - 1i 3 + 1i 0 + 2i ] #Now we define a function 'detcross' which calculates the #vector product coefficients as submatrix determinants, #analogous to euclidian (real) case #( may be found in many linear algebra textbooks) # | e1 e2 e3 | # | x(1) x(2) x(3) | # | y(1) y(2) y(3) | function V = detcross (A) V = zeros (1, 3); t0 = 1:3; for k = 1:3 t = t0; t (k) = []; V (k) = (-1) ^ (k + 1) * det (A (:, t)); endfor; endfunction #This function accepts the matrix 2x3 A=[x;y] #and returns the coefficients of the cross product as #a row vector. It works well for real A, i.e. real x,y #CASE STUDIES ( 'format short' used ) #We look at the output of the detcross - function #and check for orthogonality to input vectors #(1) we don't care that we've complex input x,y :)z=detcross( [x;y] ), S1=dot(z,x), S2=dot(z,y) z = 2.22e-16 - 2.00e+00i 4.00e+00 + 1.11e-16i 1.00e+00 + 5.00e+00i S1 = 10 - 6i S2 = 24 + 12i #Both scalar products S1, S2 unzero! So, z isn't orthogonal #to x or y and consequently cannot be the correct vector #product we're looking for! #(2) now we use a complex conjugated 1st argument x :)z=detcross( [conj(x);y] ), S1=dot(z,x), S2=dot(z,y) z = -6.00 + 4.00i 2.00 - 6.00i -1.00 - 1.00i S1 = 0.00e+00 - 2.22e-16i S2 = -24 + 12i #Looks better, because z is orthogonal to x now, #but z is not a correct cross-product yet, because #is not yet orthogonal to BOTH vectors x,y #(3) now we use a complex conjugated 2nd argument y :)z=detcross( [x;conj(y)] ), S1=dot(z,x), S2=dot(z,y) z = -6.00 - 4.00i 2.00 + 6.00i -1.00 + 1.00i S1 = -10 - 6i S2 = 0.00e+00 - 4.44e-16i #similar to (2), orthogonality to y OK now, #but orthogonality to x fails on the other hand #(4) now let's try complex conjugation of BOTH arguments :)z=detcross( [conj(x);conj(y)] ), S1=dot(z,x), S2=dot(z,y) z = 2.22e-16 + 2.00e+00i 4.00e+00 - 1.11e-16i 1.00e+00 - 5.00e+00i S1 = 0 S2 = 0.00e+00 - 2.22e-16i # BINGO :) z is orthogonal to both vectors x,y #an additional requirement of the correct cross product is it's #alternating sign if x,y are exchanged, # cross( x,y ) == - cross( y,x ) #let's check: :)z=detcross( [conj(y);conj(x)] ), S1=dot(z,x), S2=dot(z,y) z = -2.22e-16 - 2.00e+00i -4.00e+00 + 1.11e-16i -1.00e+00 + 5.00e+00i S1 = 0 S2 = 0.00e+00 + 2.22e-16i # OK :) #-----------------------------------------------------# # CONCLUSION #In order to calculate a correct (cross) vector product # z = cross ( x,y ) #in 3D unitary (complex) space #which has alternating sign under x,y exchange and #is ORTHOGONAL to BOTH vectors x,y we can use the same #formula(s) as in euclidian (real) space, but with #>> BOTH vectors x,y REPLACED by their << #>> complex conjugated counterparts conj(x),conj(y) << #-----------------------------------------------------# #Finally, let's have a (suprising) look at Octave's V2.1.33 #and Matlabs V5.0 output of the cross function for #our test vectors x,y given above #Octave 2.1.33 'cross.m' :)z=cross2133(x,y), S1=dot(z,x), S2=dot(z,y) z = 0 - 2i 4 + 0i 1 + 5i S1 = 10 - 6i S2 = 24 + 12i #This is case (1): neither orthogonal to x nor y :( #Matlab V5.1.0.421 PC #Many thanks to # Melqart[melqart at free dot fr] #for sending me the following result (I've no access to ML) :)z=cross(x,y) z = 0 - 2i 4 + 0i 1 + 5i #Same incorrect result as from Octave 2.1.33 script :( # #So, we've perfect Octave<->Matlab compatibility ! #But both platforms give incorrect results. #No complex conjugation is applied to any of the two argument #vectors with the consequence, that the vector product isn't #orthogonal to any of the two input vectors. -------------- correct version of 'cross.m' function ----------- ## -*- texinfo -*- ## at deftypefn {Function File} {} @var{z} = cross (@var{x}, @var{y}, @var{dim}) ##CROSS PRODUCT of two 3-dimensional vectors at var{x} and @var{y}. ## ##For real and also for complex at var{x}, @var{y} the @var{z} ##alternates sign(s) if the arguments are interchanged ##cross( at var{x},@var{y}) == - cross(@var{y},@var{x}) # ##Furthermore at var{z} is orthogonal to both input vectors: ##dot( at var{z},@var{x}) == dot(@var{z},@var{y}) == 0 ## ## at var{z} is returned as a row vector if @var{x} and @var{y} are ##both row vectors, otherwise, at var{z} is a column vector. ## ##The optional third argument invokes an error message that ##multidimensional arrays are not implemented in Octave. ## ## at example ## at group ## at var{z}=cross( @var{x}=[ 1 1-i 1+i ], @var{y}=[ 3-i 3+i 2i ] ) ## at result{} @var{z}=[ 2i 4 1-5i ] ## at end group ## at end example ## at end deftypefn ##AUTHOR Rolf Fabian Jan. 12, 2001 ## published under current GNU GENERAL PUBLIC LICENSE function z = cross (x, y, dim) if nargin != 2 if nargin == 3 error("cross: multidimensional arrays not implemented in Octave."); else usage("z = cross( x, y, dim )"); end endif [rx,cx]=size(x); [ry,cy]=size(y); if ( rx*cx == 3 && ry*cy == 3 ) z= conj( [x(2)*y(3)-x(3)*y(2); x(3)*y(1)-x(1)*y(3); x(1)*y(2)-x(2)*y(1) ]); if ( rx == 1 && ry == 1 ) z = z.'; endif else error ("cross: only 3-dimensional input vectors accepted"); endif endfunction ------------------------------------------------------------- Octave is freely available under the terms of the GNU GPL. Octave's home on the web: http://www.octave.org How to fund new projects: http://www.octave.org/funding.html Subscription information: http://www.octave.org/archive.html -------------------------------------------------------------