Pseudo-Euclidean Vector Spaces

by Rolf Sulanke

January 27, 1999; revised May 26, 1999;
revised November 1, 2000

Summary

This notebook contains the basic definition of vector operations in the n-dimensional pseudo-euclidean space, including the euclidean case. The dimension dim and the index ind - the number of diagonal elements equal -1 in an orthogonalized basis-  are the characterising constants for the pseudo-euclidean vector spaces, which are vector spaces over the real numbers.  In special relativity theory the "world" = "space-time"  is based on the pseudo-euclidean vector space of dimension dim = 4 and index ind = 1, corresponding to the three-dimensionality of the physical space and the one-dimensionality of time. For the N-dimensional Möbius space  we have to set dim = N +2 and ind = 1 . In the n-dimensional euclidean case we have, of course, dim = n and ind = 0. In the Lie geometry of spheres one has dim = 6 and ind = 2. 
Changes: The most interesting subject treated in this notebook is the orthogonalization of vector sequences in the pseudo-euclidean case, which has been corrected and refined in the last revision. Furthermore, the procedure psfilter has been simplified: now a repeated application of psfilter  is not necessary in any cases.

Keywords

Dimension dim, index ind, pseudo-euclidean scalar product, generalized cross product, spacelike, timelike, isotropic vectors, norming, orthogonalization, renorming, psfilter, orthopairs.

Address

sulanke@mathematik.hu-berlin.de
Homepage:
http://www-irm.mathematik.hu-berlin.de/~sulanke/

The author will be very thankful for hints, corrections, comments.

Initialization

The needed packages

Before initializing the present notebook ensure that the packages
euvec.m
pseuvec.m
and the master package init.m are in your working directory. These packages are contained in the item "Spheres".
 For the present notebook the packages mspher.m, mcirc.m, spiralsf.m are not needed.

For Windows

For the commands in this subsection:
DO NOT CONFIRM THE EVALUATION OF THE INITIALIZATION CELLS !!!

If the result of the following command

Directory[]

does not coincide with the working directory, e.g. "C:\Eigene Dateien\mymath\spheres.v4", containing the needed packages, insert the Path of your working directory in the next command line and activate it first, without calling the initialization:

SetDirectory["C:\\Eigene Dateien\\mymath\\spheres4"]

If you want to use the packages euvec.m, pseuvec.m only,  carry out the next command.

Needs["pseuvec`"]

If you want  that all commands contained in the packages of the item "Spheres" are available, activate

Get["init.m"]

Now: Activite Kernel/Evaluation/Evaluate Initialization from the menu!

For Linux (and other  systems: not tested!)

For the next command:
DO NOT CONFIRM THE EVALUATION OF THE INITIALIZATION CELLS !!!

If you want to use the package euvec.m, pseuvec.m only, remove the file init.m  from the working directory and carry out the next command.

Needs["pseuvec`"]

Now: Activite Kernel/Evaluation/Evaluate Initialization from the menu!

1. Introduction. Usages

1.1 Introduction

As we wrote this notebook, we have been interested in the Moebius geometry of the 3-sphere. Thus, the default values of dimension and index are

dim  =  5; ind = 1;

Always have in mind that almost all the defined concepts  depend on dim and ind; if you want to change the default values, simply enter the commands
dim = n; ind = m;
 where n and m, 0 ≤ m ≤ n, are your preferred values. But do this before calling the Initialization, or importing the package pseuvec.m , since in the definitions of  this package dim and ind enter as global constants.

The basic concepts used in the pseudo-euclidean vector algebra are contained  in the package pseuvec.m, which automatically imports the package euvec.m containing general definitions of  n-dimensional linear algebra.

If you want information about these  loaded concepts, enter e.g.

?ide

Often you can inform yourself about the definitions by evoking or applying them:

ide[]

The standard basis is

Table[stb[i], {i,  dim}]

and the matrix of the scalar products of the basis vectors is

Table[pssp[stb[i],  stb[j]],  {i, dim},  {j, dim}] // 
      MatrixForm

Here pssp denotes the pseudo-euclidean scalar product.

?? pssp

Since the names are prepended by the ContextNames, sometimes the definitions obtained in this way are hard to read. Therefore you will find some of them in the original form in unevaluable cells.  It is recommended to look into the packages.

All the definitions contained in the packages euvec.m and pseuvec.m are listed in the next subsection 1.2.

1.2 The usages

From the package euvec.m:

Valid for arbitrary values of the dimension dim are

?nul

?stb

?vec

?randomv

?cross

?dotnorme

?esorthonorm

Options[esorthonorm]

?normed

?print

?neglect

?orthobase

?outzero

?sterproj

?invsterproj

?hyperplane

?unit

?rank

?renorm

?smoothing

?sphericalreflection

The following commands work for dimension dim = 4 only:

?ortho1

?orthoerg

?orthoframe

?tg4frame

?subsphere

?subspheremf

From the package pseuvec.m:

?ide

?pssp

?dual

?psCross

?pscross

?ch

?pr

?orthonorm

Options[orthonorm]

?normed

?print

?neglect

?psfilter

?normalize

?orthopair

2. Vector Operations in Euclidean and Pseudo-Euclidean Spaces

2.1 Basic definitions

Of course, we can apply all the operations of real linear vector algebra. We have the null vector:

\(\(\(
\)\(nul[]\)\)\)

The standard basis is:

Table[stb[i], {i,  dim}]

The basic concept of pseudo-euclidean geometry is the indefinite scalar product pssp. 
The standard basis {stb[i]} is orthonormal with respect to this scalar product:

MatrixForm[
      Table[pssp[stb[i],  stb[j]],  {i,  dim},  {j,  dim}]]

The diagonal of this matrix is

ide[]

The dualization of a vector

vec[x]

is given by

dual[vec[x]]

This corresponds to the dual coordinates of a vector in Riemannian geometry. The scalar product pssp is related to the usual positive definite scalar product, denoted in Mathematica by  Dot ., as follows

pssp[vec[x],  vec[y]]  ==  vec[x] . dual[vec[y]]

The definition of the pseudo-orthogonal scalar product of type (dim - ind, ind) implies the following rules:

pssp[x, y] == pssp[y, x]

pssp[x + y, u] == pssp[x, u] + pssp[y, u]

pssp[a k, l] == a pssp[k, l]

The Dot product is more general. Therefore commutativity is not valid, in general:

x . y  ==  y . x

Simplify[%]

(vec[x] + vec[a]) . vec[y]  ==  
      vec[y] . vec[x]  +  vec[y] . vec[a]

Simplify[%]

2.2 Type of a vector

A non zero vector v is called  space-, time-, or light-like (in relativity theory),  if  its scalar square pssp[x,x]  is positive, negative, or zero, respectively. In geometric context one mostly uses "isotropic" instead of "lightlike". This property of the vector can be indicated by the function ch (== character), which gives the sign of pssp[v, v].

?ch

Examples

Error messages: The nullvector

ch[nul[]]

or a vector with incorrect dimension:

ch[stb[1, 3]]

Timelike vectors

ch[stb[dim]]

ch[stb[5]*Cosh[t]  +  stb[1]*Sinh[t]]

Isotropic vectors:

ch[stb[1] Cos[Pi/4]  +  stb[dim] Sin[Pi/4]]

ch[stb[1]*a +  stb[dim]*a ]

Spacelike vectors:

ch[stb[1]]

ch[stb[1]*Cosh[t]  +  stb[5]*Sinh[t]]

ch[stb[1]*Cos[t] + stb[2]*Sin[t]]

vec[x]  -  stb[dim]*x[dim]

ch[vec[x] - stb[dim]*x[dim]]

Type of a random vector

ch[randomv[]]

v1 = randomv[]

ch[v1]

The empirical probability that such a random vector is timelike equals approximately 0.06:

\(zj\  = \ 0; 
    For[i = 1, 
      i < 10000, \ \(i++\), 
	\((v1\  = \ randomv[]; 
        If[ch[v1] == \ \(-1\), \ zj\  = \ zj + 1])\)]; N[zj/10000]\)

Remove[zj,  v1, i]

Sometimes, clearly, it is impossible to decide which character  the vector has, since it depends of the value of parameters:

ch[vec[x]]

ch[stb[1] Cos[t]  +  stb[dim] Sin[t]]

Plotting shows the character in dependence of the parameter:

Plot[
      ch[stb[1] Cos[t]  +  stb[dim] Sin[t]],  {t, 0,  
        2 Pi}]

In the next example the character does not depend on the parameter; therefore it will be evaluated:

ch[stb[1] Cos[t]  +  stb[dim - ind] Sin[t]]

2.3 The norming

Now we consider the norming of the vector

normalize[vec[x]]

Since Mathematica often has difficulties with evaluations and simplifications of  Sqrt or Abs it is recommended to use normalize only in unavoidable cases, and, if possible, only  at the end of a module or calculation.

Examples

?normalize

normalize[stb[1] + 2 stb[2]]

Here you can see the problems if you want to norm general vectors:

normalize[stb[1]*Cosh[t] + stb[2]]

pssp[%, %]

Simplify[%]

Simplify[%,  t  ∈  Reals]

In newer versions, Mathematica simplifies expressions with additional conditions, see Simplify in the help system. Sometimes the result depends on the place in the sequence of evaluations:

TrueQ[1 +  Abs[ Cosh[t]]  >=  0]

Simplify[1 +  Abs[ Cosh[t]]  >=  0,  
      t  ∈  Reals]

Simplify[Abs[ Cosh[t]]  >=  0,  
      t  ∈  Complexes]

Simplify[Abs[t] >=  0,  t  ∈  Complexes]

TrueQ[Simplify[1 + Abs[Cosh[abcd]] >=  0]]

Mathematica does not return True, since generally t could be a complex parameter. But

TrueQ[Abs[1 +   Cosh[t]^2 ] >=  0]

TrueQ[Simplify[Abs[1 +   Cosh[t]^2 ]] >=  0]

is not recognized, too.  In the next example the norming works ignoring the isotropic case t=0:

normalize[stb[1]*Cosh[t] + stb[5]]

% /.  t  -> 0

Norming an isotropic vector, or the nullvector, returns the vector itself

normalize[stb[1] +  stb[5]]

normalize[stb[1]*a +  stb[5]*a]

normalize[stb[5]*a +  stb[2]*a  + stb[3 ]*b]

pssp[%, %]

normalize[nul[]]

2.4 The general cross product

As in euclidean geometry, the cross product of dim -1 vectors  (Subscript[a, i])  is the vector v uniquely defined by
pssp[v, x] = Det[(Subscript[a, 1]),..., (Subscript[a, n - 1]), x]	   (n = dim).
 First we remember our definition of  the euclidean cross product for arbitrary dimension, contained in euvec.m:

?cross

\(cross[b_List]\  := \ 
      Module[{d\  = \ Length[b]\  + \ 1, v}, 
		v\  = \ 
          nul[d]; 
			For[i = 1, i < d, \(i++\), 
          If[Length[b[\([i]\)]\ ] != \ d, \ 
								Return[
              Print[b[\([i]\)], \ " is not a correct argument!"]], 
								v[\([i]\)]\  = \ 
              Det[Append[b, \ stb[i, \ d]]]]]; 
			v[\([d]\)]\  = \ 
          Det[Append[b, \ stb[d, \ d]]]; 
			Return[v]]\)

The difference to Mathematica's function Cross is that cross' first argument is a list of lists, while Cross depends of dim -1 vectors.

?? Cross

Here are the generalizations to the pseudo-Euclidean case:

?? pscross

If  you want that the cross product depends of dim - 1 arguments, and not of the list of these arguments, you may use psCross:

psCross := Composition[dual,  Cross]

Examples

cross[{stb[1], stb[2],  stb[3], stb[4]}]

Cross[stb[1], stb[2],  stb[3], stb[4]]

pscross[{stb[1], stb[2],  stb[3], stb[4]}]

vt  =  %

Det[{stb[1], stb[2],  stb[3], stb[4], vec[x]}]

pssp[vt,  vec[x]]

pscross[{stb[5], stb[2],  stb[3], stb[4]}]

psCross[stb[1],  stb[2],  stb[3],  
      stb[3]*Cos[t] + stb[4]*Sin[t]]

pscross[{stb[1],  stb[2],  stb[3],  
        stb[3]*Cos[t] + stb[4]*Sin[t]}]

cross[{stb[1],  stb[2],  stb[3],  
        stb[3]*Cos[t] + stb[4]*Sin[t]}]

psCross[stb[5], stb[2],  stb[3], stb[4]]

cross[{stb[5], stb[2],  stb[3], stb[4]}]

pscross[{stb[2],  stb[3]*Cos[t] +  stb[4]*Sin[t],  
        stb[4]*Cos[t] -  stb[3]*Sin[t], stb[5]}]

Simplify[%]

cross[{vec[a], vec[b], vec[c], 
          vec[d]}]  ==  -cross[{vec[b], vec[a], vec[c], 
            vec[d]}]

Cross[vec[a], vec[b], vec[c], 
        vec[d]]  ==  -Cross[vec[b], vec[a], vec[c], vec[d]]

pscross[{vec[a], vec[b], vec[c], 
          vec[d]}]  ==  -pscross[{vec[b], vec[a], vec[c], 
            vec[d]}]

pscross[{vec[a], vec[b], vec[c], vec[d]}]  ==  -psCross[vec[b], 
          vec[a], vec[c], vec[d]]

3. Orthogonality

3.1 Definition

As in euclidean geometry a list {(Subscript[v, i])} of non zero vectors is called orthogonal, if pssp[(Subscript[v, i]),(Subscript[v, j]) ] = 0 for i ≠ j; it is called orthonormal if it is orthogonal and for each i one has pssp[(Subscript[v, i]), (Subscript[v, i])] = 1, pssp[(Subscript[v, i]), (Subscript[v, i])] = -1, or pssp[(Subscript[v, i]), (Subscript[v, i])] = 0 with (Subscript[v, i])≠ o. An orthonormal sequence of vectors which does not contain isotropic vectors is linearly independent, therefore its length is less or  equal dim. By the theorem of inertia the number of timelike vectors in an orthonormal sequence is less or equal the index ind. An orthonormal sequence of maximal length dim which does not contain isotropic vectors is called an orthonormal frame; in any orthonormal frame the number of timelike vectors equals ind.
In this section we collect definitions and modules for orthogonal vectors in the pseudo-euclidean space, which we need for certain applications, especially in Möbius geometry.  The basic orthonormalization procedure orthonorm is contained in the package pseuvec.m. The Erhard Schmidt orthogonalization (for positive semidefinite scalar products) esorthonorm (defined in euvec.m) is tested in the notebook eusphere.nb, and the general orthogonalization orthonorm  is treated below.  A generalization of orthonorm for scalar products in matrix spaces and applications to Killing forms can be found in the item Lie algebras, in my homepage, or in :

http://www.mathsource.com/cgi-bin/msitem?0210-845

3.2 Pseudo-euclidean orthogonalization

3.2.1 General remarks

For definite scalar products J. M.Novak's GramSchmidt orthogonalization works well. Clearly, Novak's procedure GramSchmidt (see .../mathematica/AddOns/StandardPackages/LinearAlgebra/Orthogonalization.m) was not intended for applications to indefinite scalar products, and it did not work correctly for these. For pseudo-euclidean vector spaces we programmed a procedure orthonorm with the following properties:
Given a sequence of vectors (Subscript[v, 1]), (Subscript[v, 2]),..., (Subscript[v, k]). Construct an orthogonal sequence  (Subscript[w, 1]), (Subscript[w, 2]),...,(Subscript[w, m]), which  has the same span 
                                                               span[ (Subscript[v, 1]), (Subscript[v, 2]),..., (Subscript[v, k])] = span[ (Subscript[w, 1]), (Subscript[w, 2]),...,(Subscript[w, m])]. 
To this aim we use a function defined in euvec.m which can be useful also for other applications:

?outzero

outzero[{nul[],  stb[1]}]

3.2.2 Orthonormalization

The orthogonalization is more complicated as in the euclidean case; the following procedure is contained in pseuvec.m:

?orthonorm

Options[orthonorm]

?neglect

\(orthonorm[b_List, opts___] := 
      Module[{v, w, va, vb, l, ort = {}, i, j, f, h, sm, aus, prod, nor, pr, 
          bn, follow, iso, 
          control}, {prod, nor, pr, 
            bn} = \({innerprod, normed, print, neglect} /. {opts}\) /. 
            Options[orthonorm]; 
f[h_] := 
          smoothing[h, 10^bn]; 
sm[aus_] := 
          MapAll[f, aus]; 
w = 
          outzero[sm[b]]; 
Label[
          follow]; 
l = Length[w]; 
If[
          l == 0, Goto[control]]; 
For[i = 1, 
          i ≤ l, \(i++\), 
          If[Not[Simplify[sm[prod[w[\([i]\)], w[\([i]\)]]]] == 0], 
            AppendTo[ort, w[\([i]\)]]; 
If[
              pr == True, Print[w[\([i]\)]]]; 
v = 
              Simplify[
                sm[Table[
                    w[\([j]\)] - 
                      w[\([i]\)]*
                        prod[w[\([i]\)], w[\([j]\)]]/
                          prod[w[\([i]\)], w[\([i]\)]], {j, 
                      l}]]]; 
w = outzero[v]; \ 
            Goto[follow], Null, 
            AppendTo[ort, w[\([i]\)]]; 
If[
              pr == True, Print[w[\([i]\)]]]; 
v = 
              Simplify[
                sm[Table[
                    w[\([j]\)] - 
                      w[\([i]\)]*
                        prod[w[\([i]\)], w[\([j]\)]]/
                          prod[w[\([i]\)], w[\([i]\)]], {j, 
                      l}]]]; 
w = outzero[v]; \ 
            Goto[follow]]]; 
Label[iso]; 
        l = Length[w]; 
If[l == 0, 
          Goto[control]]; 
v = w; 
If[
          l == 1, AppendTo[ort, v[\([1]\)]]; \ 
          If[pr == True, Print[v[\([1]\)]]]; 
Goto[
            control]]; 
For[i = 1, 
          i ≤ l, \(i++\), 
          For[j = i + 1, j ≤ l, \(j++\), 
            If[Not[Simplify[sm[prod[v[\([i]\)], v[\([j]\)]]]] == 0], 
              va = v[\([j]\)] - v[\([i]\)]; 
vb = 
                v[\([j]\)] + v[\([i]\)]; \ 				v[\([i]\)] = va; 
              v[\([j]\)] = vb; 
w = 
                outzero[Simplify[sm[v]]]; Goto[follow]]]]; \ 			ort = 
          Join[ort, v]; 
Label[
          control]; 
l = Length[ort]; 
If[
          l == 0, Return[ort]]; 
If[
          nor == True, 
          ort = sm[
              Table[Simplify[normalize[ort[\([i]\)], prod]], {i, 
                  Length[ort]}]]]; 
For[i = 1, 
          i ≤ l - 1, \(i++\), 
          For[j = i + 1, j ≤ l, \(j++\), 
            If[Not[Simplify[sm[prod[ort[\([i]\)], ort[\([j]\)]]]] == 
                  0], Return[
                Print["I cannot orthogonalize the given sequence!"]], 
              Null, Return[
                Print["I cannot orthogonalize the given sequence!"]]]]]; \ 	Return[ort]]\)

3.2.3  WARNINGS:

1. In certain cases the pseudo-euclidean orthonormalization procedure orthonorm cannot  be carried out. Then instead of the orthonormed matrix appears the message  "I cannot orthogonalize the given sequence!"  Sometimes the norming is the reason; then try the orthogonalization with the option normed -> False. Besides of the parameter dependence the reason can be a numerical badly conditioned matrix b (e.g. contains vectors with very small, or very different in magnitude, components). In this case one can try applying orthonorm  on the renormed matrix renorm[b].

2. In the pseudo-euclidean case the isotropic vectors are orthogonal to itself. By our definition, they are considered as being normed. Therefore in the resulting orthogonal or orthonormal sequence there may appear linearly dependent isotropic vectors.
3. In the pseudo-euclidean vector spaces exist vector sequences which may not be orthogonalized preserving its flag structure. In difference to Erhard Schmidt's orthogonalization the procedure orthonorm does not preserve the flag structure of the given sequence in general.

3.3 Examples, tests

3.3.1 Simple, parameterfree vectors

orthonorm[{nul[]}]

Isotropic vectors are treated correctly; here the orthonormal basis of the 2-space is non isotropic:

orthonorm[{stb[1] + stb[5], stb[5]}]

orthonorm[{stb[1] + stb[5], stb[2] + stb[5]}]

orthonorm[{stb[1] + stb[5], stb[1] - stb[5]}]

Sometimes appear orthonormal linearly dependent isotropic vectors:

orthonorm[{stb[2], 
        stb[1] + stb[5], (stb[1]  +  stb[5])*3}]

orthonorm[{stb[1], stb[5]  - stb[2], nul[],  
        stb[2]*5 +  stb[5]}]

Table[pssp[%[[i]], %[[j]]], {i, 3}, {j, 3}]

orthonorm[{stb[1], (stb[5]  - stb[2])*25,  
        stb[1]*5 +  stb[4]}]

orthonorm[{stb[1] + stb[4], stb[1]  -  stb[5], 
        stb[1]*3  - stb[5]}]

Table[pssp[%[[i]], %[[j]]], {i, 3}, {j, 3}]

Clearly,  by numerical reason,  one may expect an approximative result only, if during the calculation appear floating point numbers.

orthonorm[{stb[1] + stb[4], stb[1]  -  stb[5], 
        stb[1]*3.1  - stb[5]}]

Table[pssp[%[[i]], %[[j]]], {i, 3}, {j, 3}]

Here is a test with a random floating point number:

t1  =  10*Random[];
avt  = 
        orthonorm[{ stb[3], stb[5]  +  stb[2],  nul[],  
            stb[5] Cosh[t1]  +  5*stb[1] Sinh[t1], 
            stb[5]  +  stb[2]}];
\(	smoothing[
      Table[pssp[avt[\([i]\)], avt[\([j]\)]], {i, Length[avt]}, {j, 
          Length[avt]}]]\)

Repeat the evaluation of the next cell several times!

t1  =  10*Random[];
avt  = 
        orthonorm[{ stb[3], stb[5]  +  stb[2],  nul[],  
            stb[5] Cosh[t1]  +  5*stb[1] Sinh[t1], 
            stb[5]  +  stb[2]}];
\(	Table[pssp[avt[\([i]\)], avt[\([j]\)]], {i, 3}, {j, 3}]\)

Remove[avt]

3.3.2 Vectors depending on parameters

The following examples show that orthonorm worked correctly also on sequences depending on parameters.

orthonorm[{stb[
          1], (stb[5]  - stb[2])*25,  (stb[1]*5 +  stb[4])*
          a}, 
print -> True]

orthonorm[{stb[
          1], (stb[5]  - stb[2])*25,  (stb[1]*5 +  stb[4])*
          a}, 
print -> True,  
      normed  ->  False]

orthonorm[{stb[1] *c +  stb[5]*3, stb[1]  +  stb[5]}, 
      print -> True,  normed  ->  False]

orthonorm[{stb[1] *c +  stb[5]*3, stb[1]  +  stb[5]}, 
      print -> True]

Table[pssp[%[[i]], %[[j]]], {i, 2}, {j, 2}]

Simplify[Together[%]]

te[c_] =  %

te[c]

te[3]

Limit[te[c],  c  ->  3]

Limit[te[c],  c  -> 3,  
      Direction  ->  -1]

Limit[te[c],  c  ->  3,  Direction  ->  1]

The evaluation worked for general c, in spite of the fact that at c = 3  the character  of the constructed vectors changes. Moreover, for c= 3 the 2-space degenerates to an isotropic line.

\(orthonorm[{nul[], \ 3*stb[2], \ nul[], \ 
        stb[1] Cosh[t]\  + \ 2*stb[5] Sinh[t], 
\ \ \ \ stb[1]\  + \ 
          stb[2]}, \ print\  -> \ True]\)

Table[
      pssp[%[[i]], %[[j]]], {i, Length[%]}, {j, 
        Length[%]}]

ftest[t_]  =  Simplify[Together[%]]

ftest[0]

Limit[ftest[t],  t  ->  0,  Direction  -> 1]

Limit[ftest[t],  t  ->  0,  
      Direction  -> -1]

ftest[1. ]

ftest[0.1]

\(orthonorm[{nul[], \ 3*stb[2], \ nul[], \ 
        stb[2] Cosh[t]\  + \ stb[5] Sinh[t], 
\ \ \ \ stb[1]\  + \ 
          stb[2]}]\)

orthonorm[{ stb[1] Cosh[t]  +  stb[5] Sinh[t], 
        stb[1]  +   stb[2], stb[2],  nul[], 
        stb[1]  -  stb[2]}]

Table[
      pssp[%[[i]], %[[j]]], {i, Length[%]}, {j, 
        Length[%]}]

Simplify[Together[%]]

\(orthonorm[{nul[], \ 3*stb[2], \ nul[], \ 
        stb[1] Cosh[t]\  + \ stb[5] Sinh[t], 
\ \ \ \ stb[3]\  + \ 
          stb[2]}]\)

Remove[te,  ftest]

3.3.3 Index ind = 2

If the index is greater than 1, you can meet also mutual orthogonal, isotropic,  linearly independent vectors. If you wish to consider such examples, first quit the kernel! Change the first command line in subsection 1.1 according to yor wishes, and repeat the  initialization! If you have chosen ind ≠ 2,  ignore the comments in this subsection.

ch[stb[4] + stb[2]]

pssp[stb[4], stb[4]]

Here we have two isotropic orthogonal vectors spanning a totally isotropic 2-dimensional subspace:

orthonorm[{ stb[1]  +  stb[5], stb[4]  + stb[2]}]

These vectors are orthogonal and have norm 0, therefore they have been returned correctly. Moreover, they are linearly independent.

orthonorm[{ stb[1]  +  stb[5], stb[1]*2  +  stb[3], 
        stb[4]  + stb[2]}]

Table[pssp[%[[i]], %[[j]]], {i, 3}, {j, 3}]

The next two isotropic vectors span a 2-dimensional pseudoeuclidean subspace:

orthonorm[{ stb[1]  +  stb[5], stb[5]  + stb[2]}, 
      print  ->  True]

3.3.4 Tests with random vectors (case ind = 2)

ind

atest  = Table[randomv[5], {i, 5}];
ot  = orthonorm[atest]; Print[ot // MatrixForm];
Table[
      smoothing[pssp[ot[[i]], ot[[j]]]], {i, 5}, {j, 5}]

atest  = 
        PrependTo[Table[randomv[5], {i, 7}], 
          stb[1]  +  stb[5]];
ot  = orthonorm[atest]; Print[ot // MatrixForm];
Table[
      smoothing[pssp[ot[[i]], ot[[j]]]], {i, Length[ot]}, {j, 
        Length[ot]}]

atest  = Table[randomv[5], {i, 3}];
atest  =  
        Join[atest, {stb[1]  +  stb[5], stb[2]  +  stb[4]}];
ot  = orthonorm[atest]; Print[ot // MatrixForm];
Table[
      smoothing[pssp[ot[[i]], ot[[j]]]], {i, 5}, {j, 5}]

atest  = Table[randomv[5], {i, 7}];
atest  =  
        Join[{stb[1]  +  stb[5], stb[2]  +  stb[4]}, atest];
ot  = orthonorm[atest]; Print[ot // MatrixForm];
Table[
      smoothing[pssp[ot[[i]], ot[[j]]]], {i, 5}, {j, 5}]

atest  = Table[randomv[5], {i, 7}];
ot  = orthonorm[atest, print  ->  True]; 
    Print[ot // MatrixForm];
Table[
      smoothing[pssp[ot[[i]], ot[[j]]]], {i, Length[ot]}, {j, 
        Length[ot]}]

Remove[atest,  ot]

3.3.5 Renorming of badly conditioned matrices

PRESERVE ind = 2, dim = 5 !
The comments are valid in this case only. For ind =1, dim = 5 the sequence ot will be orthonormed with approximative correctness, test it later!

The next matrix ot has vectors with coordinates of very different absolute values. Applying orthonorm ,  we  do not get orthogonal vectors:

ot  = {{0.675496095887281988`, 0.838111706002798051`, 
          0.0599127084283752431`, 0.0726205437662291686`, 
          0.396285246447809988`}, {2.90602662075336581`, -0.838111706002798761`, -0.0599127084283752875`, -0.0726205437662292307`, 
          3.18523747019283787`}, {-9.19936356392773468`*^7, 
          3.3554431999999994`*^7, -6.36335108454432329`*^-10, 
          3.35544319999999984`*^7, -9.19936356392773468`*^7}}

The three vectors are "almost" orthogonal:

Table[
        smoothing[pssp[ot[[i]], ot[[j]]]], {i, Length[ot]}, {j, 
          Length[ot]}] // MatrixForm

at  = orthonorm[ot]

We try it without norming:

at  = orthonorm[ot,  normed  ->  False]

If one renorms the vectors of ot, one preserves their span, and hopefully, can find  a nearly orthonormal basis:

?renorm

\(renorm[b_?MatrixQ] := \ 
      Module[{bo\  = \ outzero[b], i\ }, 
	
        Table[bo[\([i]\)]/Sqrt[bo[\([i]\)] . bo[\([i]\)]], 
          \ {i, \ Length[bo]}]]\)

otre  =  renorm[ot]

orthonorm[otre,   print  ->  True,  
      normed  ->  False]

at  =  %;

\(Table[
        smoothing[pssp[at[\([i]\)], at[\([j]\)]]], {i, Length[at]}, {j, 
          Length[at]}] // 
\ \ MatrixForm\)

Here  we succeeded,  but if we want to get normed vectors, the procedure does not work:.

orthonorm[otre,   print  ->  True]

The reason is that the last vector may not be normalized correctly:

atn =  Table[normalize[at[[i]]], {i, 3}]

The scalar products after normalization are incorrect; therefore the procedure did not return the found vectors:

\(Table[
        smoothing[pssp[atn[\([i]\)], atn[\([j]\)]]], {i, Length[atn]}, {j, 
          Length[atn]}] // 
\ \ MatrixForm\)

The third vector of at seems to be isotropic; but copying it we get its 'inner hidden' value:

badv  = {-0.6642969596426436`, 0.24230053530593257`, 0, 
        0.24230053530593232`, -0.6642969596426437`}

pssp[badv, badv]

normalize[badv]

badnormv  =  %

pssp[badnormv, badnormv]

The hidden value of badnormv is

{-6.304594370522744`*^7, 2.2995838964643177`*^7, 0, 
      2.299583896464315`*^7, -6.304594370522745`*^7}

Here  we succeeded,  but in the next example we have certain difficulties.

ot  = {{0.675496095887281988`, 0.838111706002798051`, 
            0.0599127084283752431`, 0.0726205437662291686`, 
            0.396285246447809988`}, {2.90602662075336581`, -0.838111706002798761`, -0.0599127084283752875`, -0.0726205437662292307`, 
            3.18523747019283787`}, {-9.19936356392773468`*^7, 
            3.3554431999999994`*^7, -6.36335108454432329`*^-10, 
            3.35544319999999984`*^7, -9.19936356392773468`*^7}, {6.82635841776287133`*^7, -2.48989593404715137`*^7, 
            9.9150403743008848`*^-9, -2.48989593404715048`*^7, 
            6.82635841776287133`*^7}, {-1.82640220915799815`*^7, 
            6.66175309911068413`*^6, 2.6549032889775681`*^-7, 
            6.66175309911060686`*^6, -1.82640220915799815`*^7}};

This matrix is bad conditioned too, but not nearly orthogonal:

Table[
        smoothing[pssp[ot[[i]], ot[[j]]]], {i, Length[ot]}, {j, 
          Length[ot]}] // MatrixForm

orthonorm[ot]

otre  =  renorm[ot]

The matrix  otre is badly conditioned again: The last three vectors are pairwise almost linearly dependent, and we have

Det[otre]

Table[pssp[otre[[i]], otre[[j]]], {i, 5}, {j, 5}]

With the default norming we do not succeed:

at1  =  orthonorm[otre,  print  ->  True]

but without norming we get a reasonable result:

at1  =  
      orthonorm[otre,  print  ->  True, 
        normed  ->  False]

Table[
        smoothing[pssp[at1[[i]], at1[[j]]]], {i, Length[at1]}, {j, 
          Length[at1]}] // MatrixForm

That looks wrong: in the pseudo-euclidean vector space of index 2 there may not exist three mutual orthogonal linearly independent isotropic vectors! The numerical deviation from orthonormality one can see evaluating the scalar products not using smoothing:

Table[
      pssp[at1[[i]], at1[[j]]], {i, Length[at1]}, {j, 
        Length[at1]}]

Thus, orthonorm yields very small, almost orthogonal vectors, and the last three vectors are almost isotropic; the five vectors are "almost linearly dependent":

Det[at1]

3.4  psfilter

We now go back to index 1: Quit kernel, set dim = 5, ind = 1 in subsection 1.1,  and repeat the Initialization!

{dim ,  ind}

Remember that in a space of index ind  any orthogonal sequence of linearly independent vectors may contain not more than ind non spacelike vectors.

Sometimes it is necessary to place the spacelike vectors of a sequence b at the first places. The following procedure  realizes such a permutation:

?psfilter

\(psfilter[b_List]\  := \ 
      Module[{i, \ lb\  = \ Length[b], \ bn = \ {}, qn\  = \ {}}, 
		
        For[i = 1, i <= lb, \(i++\), 
			
          If[Length[b[\([i]\)]\ ]\  != \ dim, \ 
            Return[\ Print["Incorrect argument ", b[\([i]\)], 
                "!!!"]]]]; 
		
        For[i = 1, \ i <= lb, \ \(i++\), 
			
          If[ch[b[\([i]\)]\ ] == \ 1, \ AppendTo[bn, \ b[\([i]\)]], 
            AppendTo[qn, b[\([i]\)]], AppendTo[qn, b[\([i]\)]]]]; 
				
        Return[Join[bn, qn]]]\)

\(\(\(		\)\(psfilter[{a}]\)\)\)

\(\(\(		\)\(psfilter[{nul[], \ stb[5], \ stb[1]}]\)\)\)

psfilter  yields the expected result also for non orthogonal sequences:

\(\(\(		\)\(psfilter[{stb[2], stb[5]*Cos[t], \ stb[5], \ 
        stb[1]}]\)\)\)

psfilter[%]

psfilter[{ stb[1]  +  stb[5], 
        stb[5]  + stb[2],  {Cos[t],  Sin[t],  1, 0, 1}}]

3.5 Orthogonal pairs

Now we generate a pair of orthonormal spacelike vectors starting from a given spacelike one; that we need for the geometry of circles on the 3-sphere.The definition is contained in the  package pseuvec.m.

?orthopair

\(orthopair[l_List] := 
      Module[{a, l1 = l, m1 = nul[], b = dim - ind, j, k}, 
          a = Simplify[pssp[l, l]]; 
If[a ≤ 0, 
            Print["Vector ", l, " is not spacelike!"]; Return[], 
            Null, Print["Maybe, Vector ", l, " is not spacelike?"]; 
            Return[]]; \ 
          If[b < 2, Print["Orthogonal spacelike pairs don't exist!"]; 
            Return[]]; \ 	\ l1 = l1/Sqrt[a]; \ 	For[k = 1, 
            k ≤ b, \(k++\), 
            If[l1[\([k]\)] == 0, m1 = stb[k]; \ 
              Return[{l1, m1}]]]; \ 	\ For[j = 1, j ≤ b, \(j++\), 
            If[OddQ[j] && j < b, 
              m1[\([j]\)] = \(-l1[\([j + 1]\)]\)]; \ \ 		If[EvenQ[j], 
              m1[\([j]\)] = l1[\([j - 1]\)]]]; \ 		m1 = 
            m1/Sqrt[pssp[m1, m1]]; 
\ {l1, m1}] /; 
        Length[l] == dim\)

The following function controls orthonormality:

control[v_, w_] := 
      Simplify[(((pssp[v, v] - 1))^2) + pssp[v, w]^2 + 
          (((pssp[w, w] - 1))^2)]

?control

If all is correct the result mus be 0, or numerically very small. If the result depends on a parameter, integrate!

Tests

orthopair[stb[1]]

orthopair[stb[2]]

orthopair[stb[3]]

orthopair[stb[4]]

orthopair[stb[5]]

orthopair[{Cos[t], Sin[t], Cos[s], Sin[s], 0}/(Sqrt[2])]

Simplify[%]

The argument of orthopair must be a spacelike vector, but it needs not to be a unit vector:

orthopair[{Cos[t], Sin[t], Cos[s], Sin[s], 0}]

Simplify[%]

Simplify[orthopair[{1, Cos[t], Sin[t], Cos[s], Sin[s]}]]

In spite of the fact that the argument is spacelike for all values of s and t!

pssp[{1, Cos[t], Sin[t], Cos[s], Sin[s]}, {1, Cos[t], Sin[t], Cos[s], 
        Sin[s]}]

orthopair[ 3*randomv[ ]]

Table[
      Simplify[pssp[%[[i]],  %[[j]]]], {i, 2}, {j, 2}]

control[%%[[1]],  %%[[2]]]