Álgebra Lineal
Colin Fahey
1. Software
LinearAlgebra.zip
Álgebra lineal código fuente (C#)
19910 bytes
MD5: 11d8c8035cac30ba543e5e0b72ee9767
2. Introducción
Este artículo describe los vectores y matrices en (d) dimensiones espaciales.
3. (d) dimensiones espaciales: atributos
3.1 Array
“serie:” Una colección de variables de tal forma que cada variable tiene un nombre único, y de tal manera que los nombres se le puede asignar un orden.
Integer valores pueden ser utilizados como nombres para las variables en una matriz.
Por ejemplo, si un conjunto contiene (d) variables, entonces la { 0, 1, 2, ..., (d-1) } valores pueden ser los nombres asignados a las variables en la matriz.
3.2 (d)-dimensional Vector
“(d)-dimensional vector:” Una serie de variables (d).
“vector componente:” Una variable dentro de un vector.
3.3 (d)-dimensional vector espacial
“una dimensión espacial:” El conjunto completo de valores que pueden almacenarse de una variable.
“(d) dimensiones del espacio:” El conjunto completo de combinaciones de valores que pueden ser almacenados por un conjunto de variables (d).
La definición formal de un “espacio vectorial:”
Que (T) ser un tipo básico (por ejemplo: número real, entero, número complejo, número racional, etc.)
Cualquier variable de un tipo básico se llama un “escalar.”
A “(T) de tipo (d)-dimensional vector espacial” es el conjunto (S) de (d) dimensiones vectores tener dos operaciones, además (+) vector, y la multiplicación escalar (*), que cumpla las condiciones siguientes:
(1) Si (v) y (w) son dos vectores en (S), entonces (v + w) es también un vector en (S);
(2) Si (u), (v), y (w) son tres vectores en (S), entonces (u + v) + w = u + (v + w);
[aditivo conmutatividad]
(3) Si (v) y (w) son dos vectores en (S), entonces (v + w) = (w + v);
[aditivo asociatividad]
(4) Hay un “vector cero,” (0), en (S), de tal manera que para cualquier vector (v) en (S), (v + (0)) = v;
[aditivo identidad]
(5) Si es (c) escalar cualquier tipo de (T), y (v) es cualquier vector en (S), entonces el producto (c * v) es un vector en (S);
(6) Si (a), (b), y (c) son escalares de cualquier tipo (T), y (v) y (w) son vectores en cualquier (S), entonces (a + b) * v = a*v + b*v, y c*(v + w) = c*v + c*w;
[multiplicativo distributivity]
(7) Si (a) y (b) son escalares de cualquier tipo (T), y (v) es cualquier vector en (S), entonces (a*b)*v = a*(b*v);
(8) Si “1” es un tipo de escalar (T) tal que (1*1)=1, y (v) es cualquier vector en (S), entonces (1*v) = v;
(9) Para cada vector (v) en (S), el vector (-1)*v = -v cumple v + (-v) = (0);
[aditivo inverso]
3.4 (d)-dimensional vector código
El código de abajo muestra cómo un (d)-dimensional vector, con 64 bits de punto flotante componentes, se pueden aplicar.
Un array es suficiente para representar un vector.
El código que aparece a continuación tiene un conjunto que figura en una clase, sólo por conveniencia.
El código no está destinado a ser eficiente.
Una estructura (por ejemplo: “struct”, un tipo de valor) en representación de los vectores de un número fijo de dimensiones (ejemplos: 3 o 4) es probable que sea mucho más eficiente que el general (d)-dimensional vector de clase se muestra aquí.
Aunque el código que aparece a continuación define un vector con punto flotante componentes, este documento también hace uso de vectores con componentes entero.
El código que aparece a continuación pueden ser fácilmente modificados para aplicar los vectores con componentes entero.
using System;
using System.Collections.Generic;
using System.Text;
// Multidimensional vector with 64-bit floating-point components:
public class VectorF64
{
private double[] components;
public int Dimensions()
{
if (null == this.components)
{
return (0);
}
return (this.components.Length);
}
public VectorF64()
{
this.components = null;
}
public VectorF64( params double[] paramValues )
{
this.components = null;
if (null == paramValues)
{
return;
}
int dimensions = paramValues.Length;
this.components = new double[dimensions];
for (int i = 0; i < dimensions; i++)
{
this.components[i] = paramValues[i];
}
}
public VectorF64( VectorF64 other )
{
this.components = null;
if (null == other)
{
return;
}
if (null == other.components)
{
return;
}
int dimensions = other.Dimensions();
this.components = new double[dimensions];
for (int i = 0; i < dimensions; i++)
{
this.components[i] = other.components[i];
}
}
public double this[int index]
{
get
{
if (null == this.components)
{
return (0.0);
}
if
(
(index >= 0)
&& (index < this.components.Length)
)
{
return (this.components[index]);
}
return (0.0);
}
set
{
if (null == this.components)
{
return;
}
if
(
(index >= 0)
&& (index < this.components.Length)
)
{
this.components[index] = value;
}
}
}
public void Write( int precision )
{
if (null == this.components)
{
Log.Write( String.Empty + '(' + ' ' + ')' );
return;
}
int dimensions = this.Dimensions();
if (0 == dimensions)
{
Log.Write( String.Empty + '(' + ' ' + ')' );
return;
}
// Determine the largest component width in characters
// so that we can make all components an equal width.
int largestComponentWidth = 1;
for (int i = 0; i < dimensions; i++)
{
// { index [,minwidth] [:typeCode[precision]] }
// (minwidth<0) means left-justify
String text = String.Format( String.Empty + '{' + '0' + ':'
+ 'g' + precision + '}', this[i] );
if (text.Length > largestComponentWidth)
{
largestComponentWidth = text.Length;
}
}
Log.Write( '(' );
for (int i = 0; i < dimensions; i++)
{
Log.Write( ' ' );
String text =
String.Format( String.Empty + '{' + '0' + ','
+ largestComponentWidth + ':' + 'g' + precision + '}',
this[i] );
Log.Write( text );
if ((i + 1) < dimensions)
{
Log.Write( ',' );
}
else
{
Log.Write( ' ' );
}
}
Log.Write( ')' );
}
public void WriteLine( int precision )
{
this.Write( precision );
Log.WriteLine();
}
public void WriteLine()
{
const int defaultPrecision = 8;
this.Write( defaultPrecision );
Log.WriteLine();
}
// . . .
public static void Test()
{
// A 3-dimensional vector with 64-bit floating-point components:
VectorF64 v3 = new VectorF64( 0.0, 1.0, 2.0 );
v3.WriteLine(); // ( 0, 1, 2 )
// A 4-dimensional vector with 64-bit floating-point components:
VectorF64 v4 = new VectorF64( 0.0, 1.0, 2.0, 3.0 );
v4.WriteLine(); // ( 0, 1, 2, 3 )
// . . .
}
}
“(d)-dimensional vector cero:” un vector con todos los componentes (d) igual a cero.
public class VectorF64
{
// . . .
public static VectorF64 Zero( int dimensions )
{
VectorF64 zero = new VectorF64();
zero.components = new double[dimensions];
for (int i = 0; i < dimensions; i++)
{
zero[i] = 0.0;
}
return (zero);
}
// . . .
public static void Test()
{
// . . .
// An 8-dimensional vector with all 8 64-bit floating-point
// components set to zero:
VectorF64 z = VectorF64.Zero( 8 );
z.WriteLine(); // ( 0, 0, 0, 0, 0, 0, 0, 0 )
// . . .
}
}
3.5 (d) dimensiones del espacio: puntos
“(d) dimensiones del espacio:” Una serie de (d) variables con valores específicos “(coordinar valores).”
“(d) dimensiones espaciales origen:” Una serie de variables (d) con todos los valores igual a cero.
3.6 (d) dimensiones vectores: no relativo y relativo
A (d)-dimensional vector que “no” es “pariente” (d) es un vector-dimensional que directamente representa un estado o configuración.
A (d)-dimensional vector que es “relativo” (d) es un vector-dimensional que representa los cambios a un conjunto de componentes.
Un familiar vector puede representar la diferencia entre los dos Estados no relativo vectores.
Habida cuenta de un pariente de vectores, la determinación de un no pariente estado o la configuración usando vectores relativa que exige la combinación de vectores que en relación con un no pariente vector.
No relativo vectores y relativos son los dos vectores de vectores.
Si un vector no es pariente o familiar debe ser especificado cuando los vectores se define.
Si un (d)-dimensional vector se interpreta como no pariente, entonces la (d)-dimensional vector puede representar un punto en (d) dimensiones espaciales.
3.7 (d)-dimensional vector suma, resta, y la ampliación
Vector suma, resta, y la ampliación:
public class VectorF64
{
// . . .
public static VectorF64 operator +( VectorF64 a, VectorF64 b )
{
if ((null == a) || (null == b))
{
return (new VectorF64()); // Vector not specified.
}
if ((null == a.components) || (null == b.components))
{
return (new VectorF64()); // Vector is empty.
}
if (a.Dimensions() != b.Dimensions())
{
return (new VectorF64()); // Vectors not the same size.
}
int dimensions = a.Dimensions();
VectorF64 result = VectorF64.Zero( dimensions );
for (int i = 0; i < dimensions; i++)
{
result[i] = a[i] + b[i];
}
return (result);
}
public static VectorF64 operator -( VectorF64 a, VectorF64 b )
{
if ((null == a) || (null == b))
{
return (new VectorF64()); // Vector not specified.
}
if ((null == a.components) || (null == b.components))
{
return (new VectorF64()); // Vector is empty.
}
if (a.Dimensions() != b.Dimensions())
{
return (new VectorF64()); // Vectors not the same size.
}
int dimensions = a.Dimensions();
VectorF64 result = VectorF64.Zero( dimensions );
for (int i = 0; i < dimensions; i++)
{
result[i] = a[i] - b[i];
}
return (result);
}
public static VectorF64 operator -( VectorF64 a )
{
if (null == a)
{
return (new VectorF64()); // Vector not specified.
}
if (null == a.components)
{
return (new VectorF64()); // Vector is empty.
}
int dimensions = a.Dimensions();
VectorF64 result = VectorF64.Zero( dimensions );
for (int i = 0; i < dimensions; i++)
{
result[i] = (-( a[i] ));
}
return (result);
}
public static VectorF64 operator *( double scale, VectorF64 a )
{
if (null == a)
{
return (new VectorF64()); // Vector not specified.
}
if (null == a.components)
{
return (new VectorF64()); // Vector is empty.
}
int dimensions = a.Dimensions();
VectorF64 result = VectorF64.Zero( dimensions );
for (int i = 0; i < dimensions; i++)
{
result[i] = scale * a[i];
}
return (result);
}
public static VectorF64 operator *( VectorF64 a, double scale )
{
if (null == a)
{
return (new VectorF64()); // Vector not specified.
}
if (null == a.components)
{
return (new VectorF64()); // Vector is empty.
}
int dimensions = a.Dimensions();
VectorF64 result = VectorF64.Zero( dimensions );
for (int i = 0; i < dimensions; i++)
{
result[i] = scale * a[i];
}
return (result);
}
// . . .
public static void Test()
{
// . . .
// Examples of vector addition, subtraction, and scaling:
VectorF64 a = new VectorF64( 0.0, 1.0, 2.0, 3.0 );
a.WriteLine(); // ( 0, 1, 2, 3 )
VectorF64 b = new VectorF64( 3.0, 2.0, 1.0, 0.0 );
b.WriteLine(); // ( 3, 2, 1, 0 )
VectorF64 c = new VectorF64();
c.WriteLine(); // ( )
c = a + b;
c.WriteLine(); // ( 3, 3, 3, 3 )
c = a - b;
c.WriteLine(); // ( -3, -1, 1, 3 )
c = -b;
c.WriteLine(); // ( -3, -2, -1, 0 )
c = 3.0 * a;
c.WriteLine(); // ( 0, 3, 6, 9 )
// . . .
}
}
3.8 (d)-dimensional base vectores
La definición formal de una “base” de un espacio vectorial:
Que (T) ser un tipo básico (por ejemplo: número real, entero, número complejo, número racional, etc.)
Cualquier variable de un tipo básico se llama un “escalar.”
Que (V) ser un “(T) de tipo (d)-dimensional espacios vectoriales.”
Si el no-cero vectores { u1, u2, ..., ud } en (V) son tales que cada vector (v) en (V) puede escribirse como una “combinación lineal” de los vectores, v = c1*u1 + c2*u2 + ... + cd*ud, donde { c1, c2, ..., cd } son escalares de tipo (T), entonces es (V) “abarcó” { u1, u2, ..., ud } por vectores.
Cualquier conjunto de no-cero vectores { u1, u2, ..., ud } que “abarcan” espacios vectoriales (V) es nombrado una “base” de (V).
Una simple “base” de una (d)-dimensional vector espacial es el conjunto de (d) distintas dimensiones (d) vectores, cada una con un elemento igual a uno y todos los demás componentes iguales a cero.
Esa base son vectores “ortonormales,” lo que significa que son mutuamente perpendicular “(ortogonal)” y que cada vector tiene unidad de longitud.
Cada uno de esos vectores es un vector unitario paralelo a uno de los ejes de coordenadas (d).
Expresando un vector arbitrario como una “combinación lineal” de estos vectores base es directa; cada uno de los componentes de la arbitrariedad vector se multiplica por la correspondiente base de vectores, y estos productos se suman a la forma arbitraria de vectores.
El código de abajo muestra cómo un vector puede ser expresado como una “combinación lineal” de vectores base.
El código que aparece a continuación define arbitrariamente ortonormales un conjunto de vectores base.
public class VectorF64
{
// . . .
public static VectorF64 BasisVector( int dimensions, int componentIndex )
{
if (dimensions < 0)
{
// Invalid number of dimensions specified.
return (new VectorF64());
}
VectorF64 basisVector = VectorF64.Zero( dimensions );
if ((componentIndex >= 0) && (componentIndex < dimensions))
{
basisVector[ componentIndex ] = 1.0;
}
return (basisVector);
}
// . . .
public static void Test()
{
// . . .
// The 4 basis vectors of 4-dimensional space,
// each with 4 64-bit floating-point components:
VectorF64 b0 = VectorF64.BasisVector( 4, 0 );
b0.WriteLine(); // ( 1, 0, 0, 0 )
VectorF64 b1 = VectorF64.BasisVector( 4, 1 );
b1.WriteLine(); // ( 0, 1, 0, 0 )
VectorF64 b2 = VectorF64.BasisVector( 4, 2 );
b2.WriteLine(); // ( 0, 0, 1, 0 )
VectorF64 b3 = VectorF64.BasisVector( 4, 3 );
b3.WriteLine(); // ( 0, 0, 0, 1 )
// . . .
}
}
Cualquier vector en (d) dimensiones del espacio se puede expresar como la suma de los productos de base y los números de vectores:
public class VectorF64
{
// . . .
public static void Test()
{
// . . .
// The following two vectors are equivalent:
// A 4-dimensional vector with 64-bit floating-point components:
VectorF64 va = new VectorF64( 0.1, 1.1, 2.2, 3.3 );
va.WriteLine(); // ( 0.1, 1.1, 2.2, 3.3 )
// A 4-dimensional vector formed by scaling the 4 independent
// basis vectors for 4-dimensional space:
VectorF64 vb =
0.1 * VectorF64.BasisVector( 4, 0 )
+ 1.1 * VectorF64.BasisVector( 4, 1 )
+ 2.2 * VectorF64.BasisVector( 4, 2 )
+ 3.3 * VectorF64.BasisVector( 4, 3 );
vb.WriteLine(); // ( 0.1, 1.1, 2.2, 3.3 )
// . . .
}
}
3.9 (d) dimensiones del espacio: la distancia entre los puntos
Que (P) ser un (d) vector-dimensional que representa un punto en (d) dimensiones espaciales.
Que (Q) ser un (d) vector-dimensional que representa un punto en (d) dimensiones espaciales.
Que (R) ser un (d)-dimensional vector que representa el cambio en las coordenadas (d) a punto de llegar a (P) punto (Q); R = (Q - P).
En una dimensión espacial, P = ( p0 ), Q = ( q0 ), y R = (Q - P) = ( q0 - p0 ).
La distancia entre los dos puntos es la siguiente: Abs( q0 - p0 ).
En dos dimensiones espaciales, P = ( p0, p1 ), Q = ( q0, q1 ), y R = (Q - P) = ( q0-p0, q1-p1 ).
Interpretación de los dos desplazamientos perpendiculares como la perpendicular lados de un triángulo derecho-, la distancia entre los puntos corresponde a la longitud de la hipotenusa de dicho triángulo.
La fórmula Pythagorean, (a*a) + (b*b) = (c*c), donde (a) y (b) son las longitudes de los lados perpendiculares de un triángulo derecho-, y (c) es la longitud de la hipotenusa (skew lado), puede ser usado para determinar la distancia entre los dos puntos: Sqrt( Sq(q0-p0) + Sq(q1-p1) ).
En el espacio tridimensional, P = ( p0, p1, p2 ), Q = ( q0, q1, q2 ), y R = (Q - P) = ( q0-p0, q1-p1, q2-p2 ).
Interpretar el desplazamiento ( q0-p0, q1-p1, 0 ) como la perpendicular lados de un triángulo derecho-, y utilizando la fórmula Pythagorean, da la distancia entre el punto (P) y el punto ( q0, q1, p2 ): d01 = Sqrt( Sq(q0-p0) + Sq(q1-p1) ).
El desplazamiento ( q0-p0, q1-p1, 0 ) es perpendicular al desplazamiento ( 0, 0, q2-p2 ), y de otro derecho-puede ser triángulo formado, y la Pythagorean fórmula puede ser utilizado de nuevo.
De este modo, la distancia desde el punto (P) a punto (Q) viene dada por: Sqrt( Sq(d01) + Sq(q2-p2) ) = Sqrt( (Sq(q0-p0) + Sq(q1-p1)) + Sq(q2-p2) ) = Sqrt( Sq(q0-p0) + Sq(q1-p1) + Sq(q2-p2) ).
El método de extender la fórmula de distancia de dos dimensiones espacio a espacio tridimensional puede aplicarse varias veces para finalmente determinar la distancia fórmula para (d) dimensiones espaciales: Sqrt( (Sq(q0-p0) + Sq(q1-p1)) + Sq(q2-p2) + ... + Sq(qd-pd) ).
El código que aparece a continuación define una función llamada “Largo” que calcula la longitud de un (d)-dimensional vector.
Cuando un vector representa el desplazamiento entre dos puntos en (d) dimensiones espaciales, la longitud del vector representa la distancia entre esos dos puntos.
public class VectorF64
{
// . . .
public double Length()
{
if (null == this.components)
{
return (0.0); // Vector empty.
}
int dimensions = this.Dimensions();
double sumOfSquares = 0.0;
for (int i = 0; i < dimensions; i++)
{
sumOfSquares += (this[i] * this[i]);
}
double length = Math.Sqrt( sumOfSquares );
return (length);
}
public static double Length( VectorF64 a )
{
if (null == a)
{
return (0.0); // Vector not specified.
}
if (null == a.components)
{
return (0.0); // Vector is empty.
}
return (a.Length());
}
// . . .
public static void Test()
{
// . . .
// Example of vector length:
// A 6-dimensional vector representing a point (p):
VectorF64 p = new VectorF64( 0.0, 1.0, 2.0, 3.0, 4.0, 5.0 );
p.WriteLine(); // ( 0, 1, 2, 3, 4, 5 )
// A 6-dimensional vector representing a point (q):
VectorF64 q = new VectorF64( -5.0, 4.0, -3.0, 2.0, -1.0, 0.0 );
q.WriteLine(); // ( -5, 4, -3, 2, -1, 0 )
// A 6-dimensional vector representing the displacement
// from point (p) to point (q):
VectorF64 r = q - p;
r.WriteLine(); // ( -5, 3, -5, -1, -5, -5 )
// The distance between point (p) and point (q)
// in 6-dimensional space:
double distance = r.Length();
Log.WriteLine( distance ); // 10.4880884817015
// . . .
}
}
3.10 (d) dimensiones vectores: producto escalar
El “producto escalar” convierte (d) dos dimensiones vectores a un número.
El código que aparece a continuación calcula un producto escalar de dos vectores:
public class VectorF64
{
// . . .
public static double Dot( VectorF64 a, VectorF64 b )
{
if ((null == a) || (null == b))
{
return (0.0); // Vector not specified.
}
if ((null == a.components) || (null == b.components))
{
return (0.0); // Vector is empty.
}
if (a.Dimensions() != b.Dimensions())
{
return (0.0); // Vectors not the same size.
}
int dimensions = a.Dimensions();
double dotProduct = 0.0;
for (int i = 0; i < dimensions; i++)
{
dotProduct += (a[i] * b[i]);
}
return (dotProduct);
}
// . . .
}
Para cualquier vector, (A), Length(A) = Sqrt(Dot(A,A)).
3.11 (d) dimensiones vectores: definición de “paralelas”
Vectores (A) y (B) son “paralelos” si todos se dan las circunstancias siguientes:
(1) A.Length() > 0;
(2) B.Length() > 0;
(3) Abs(Dot(A,B)) = A.Length()*B.Length().
El código siguiente determina si un par de vectores son paralelos (posiblemente anti-alineados).
Números de punto flotante puede acumular errores en la parte fraccional debido a la limitada precisión, y, por tanto, el código debería incluir la no-cero tolerancias al comparar los números de punto flotante.
El código de ejemplo incluye valores de tolerancia, pero el ejemplo valores de tolerancia podría no ser apropiado para algunas tareas.
public class VectorF64
{
// . . .
public static bool Parallel( VectorF64 a, VectorF64 b )
{
// Smallest normalized float : 1.1754943e-38
// Smallest normalized double: 2.2250738585072020e-308
double nonZeroThreshold = 1.0e-38; // conservative for double
// double: (52+1)-bit mantissa; log10(2^53)=15.95 decimal digits
double fractionalDifferenceThreshold = 1.0e-14; // conservative
if ((null == a) || (null == b))
{
return (false); // Vector is not specified.
}
if ((null == a.components) || (null == b.components))
{
return (false); // Vector is empty.
}
if (a.Dimensions() != b.Dimensions())
{
return (false); // Vectors not the same size.
}
double lengthA = a.Length();
if (lengthA <= nonZeroThreshold)
{
return (false);
}
double lengthB = b.Length();
if (lengthB <= nonZeroThreshold)
{
return (false);
}
double oneImpliesParallel =
Math.Abs( VectorF64.Dot( a, b ) ) / (lengthA * lengthB);
double absoluteDifferenceFromOne =
Math.Abs( oneImpliesParallel - 1.0 );
if (absoluteDifferenceFromOne <= fractionalDifferenceThreshold)
{
return (true);
}
return (false);
}
// . . .
public static void Test()
{
// . . .
// Example of testing for parallel vectors:
// A 6-dimensional vector:
VectorF64 vf = new VectorF64( 0.0, 1.0, 2.0, 3.0, 4.0, 5.0 );
vf.WriteLine(); // ( 0, 1, 2, 3, 4, 5 )
// A 6-dimensional vector:
VectorF64 vg = new VectorF64( 0.0, -2.0, -4.0, -6.0, -8.0, -10.0 );
vg.WriteLine(); // ( 0, -2, -4, -6, -8, -10 )
// Determine if the specified vectors are parallel
// (or "anti-aligned"):
bool parallel = VectorF64.Parallel( vf, vg );
Log.WriteLine( parallel ); // True
// Add a non-negligible displacement to a component of a vector:
vf[0] += 1.0e-5;
vf.WriteLine(); // ( 1E-05, 1, 2, 3, 4, 5 )
// Determine if the specified vectors are parallel
// (or "anti-aligned"):
parallel = VectorF64.Parallel( vf, vg );
Log.WriteLine( parallel ); // False
// . . .
}
}
3.12 (d) dimensiones vectores: definición de la “perpendicular”
Vectores (A) y (B) son “perpendiculares” si todos se dan las circunstancias siguientes:
(1) A.Length() > 0;
(2) B.Length() > 0;
(3) Abs(Dot(A,B)) = 0.
El código siguiente determina si un par de vectores son perpendiculares. Números de punto flotante puede acumular errores en la parte fraccional debido a la limitada precisión, y, por tanto, el código debería incluir la no-cero tolerancias al comparar los números de punto flotante.
El código de ejemplo incluye valores de tolerancia, pero el ejemplo valores de tolerancia podría no ser apropiado para algunas tareas.
public class VectorF64
{
// . . .
public static bool Perpendicular( VectorF64 a, VectorF64 b )
{
// Smallest normalized float : 1.1754943e-38
// Smallest normalized double: 2.2250738585072020e-308
double nonZeroThreshold = 1.0e-38; // conservative for double
// double: (52+1)-bit mantissa; log10(2^53)=15.95 decimal digits
double fractionalDifferenceThreshold = 1.0e-14; // conservative
if ((null == a) || (null == b))
{
return (false); // Vector is not specified.
}
if ((null == a.components) || (null == b.components))
{
return (false); // Vector is empty.
}
if (a.Dimensions() != b.Dimensions())
{
return (false); // Vectors not the same size.
}
double lengthA = a.Length();
if (lengthA <= nonZeroThreshold)
{
return (false);
}
double lengthB = b.Length();
if (lengthB <= nonZeroThreshold)
{
return (false);
}
double zeroImpliesPerpendicular =
Math.Abs( VectorF64.Dot( a, b ) ) / (lengthA * lengthB);
double absoluteDifferenceFromZero =
Math.Abs( zeroImpliesPerpendicular );
if (absoluteDifferenceFromZero <= fractionalDifferenceThreshold)
{
return (true);
}
return (false);
}
// . . .
public static void Test()
{
// . . .
// Example of testing for perpendicular vectors:
// A 6-dimensional vector:
VectorF64 vf2 = new VectorF64( 0.0, 1.0, 2.0, 0.0, 4.0, 5.0 );
vf2.WriteLine(); // ( 0, 1, 2, 0, 4, 5 )
// A 6-dimensional vector:
VectorF64 vg2 = new VectorF64( 10.0, 0.0, 0.0, -5.0, 0.0, 0.0 );
vg2.WriteLine(); // ( 10, 0, 0, -5, 0, 0 )
// Determine if the specified vectors are perpendicular
bool perpendicular = VectorF64.Perpendicular( vf2, vg2 );
Log.WriteLine( perpendicular ); // True
// Add a non-negligible displacement to a component of a vector:
vf2[0] += 1.0e-13;
vf2.WriteLine(); // ( 1E-13, 1, 2, 0, 4, 5 )
// Determine if the specified vectors are perpendicular
perpendicular = VectorF64.Perpendicular( vf2, vg2 );
Log.WriteLine( perpendicular ); // False
// . . .
}
}
3.13 Matrices
“matriz:” Una colección de variables de tal forma que cada variable tiene una combinación única de una “fila” nombre y un nombre de “columna.”
“entrada:” Una variable dentro de una matriz.
Integer valores pueden ser utilizados como nombres de “fila” y “columna de” nombres para las variables en una matriz.
Por ejemplo, si una matriz (totalRows) ha (totalColumns) filas y columnas, entonces la { 0, 1, ..., (totalRows-1) } valores pueden ser los nombres asignados a las filas, y los valores pueden ser { 0, 1, ..., (totalColumns-1) } los nombres asignados a las columnas.
Por lo tanto, un variables en una matriz puede ser especificado por la especificación de un par de enteros, ( row, column ), lo que indica la combinación de fila y columna índices correspondientes a la variable concreta.
El tamaño de una matriz se especifica como “(totalRows) * (totalColumns)” (o “(totalRows) de (totalColumns)).”
Este fin de dimensiones es el mismo que el orden de dimensiones utilizado para especificar las entradas en la matriz “(( row, column )).”
[Esta convención es algo lamentable, porque para muchos de dos dimensiones usos (ejemplos: imágenes, gráficos, etc) la convención habitual es especificar las dimensiones como “width * height” y coordina como “( horizontal, vertical )” (o “( x, y )).”
Esto es lo contrario de la orden de las dimensiones y las coordenadas utilizadas para describir matrices y sus entradas. ]
Una matriz con (totalRows) igual a (totalColumns) se llama “cuadrado,” de lo contrario, la matriz se denomina “rectangular.”
Una matriz puede ser considerado como que contiene un conjunto de “vectores fila,” donde las variables en cada fila se interpretan como pertenecientes a un vector.
Una matriz puede también ser considerado como que contiene un conjunto de “vectores columna,” donde las variables en cada columna se interpretan como pertenecientes a un vector.
Las matrices pueden representar una gran variedad de relaciones matemáticas.
El significado de una matriz, y las operaciones que podrían ser apropiadas para su procesamiento entradas de una matriz, depende del contexto.
Sin embargo, hay reglas básicas de aritmética matriz que son relevantes para muchos contextos, y estas normas básicas se definirán en una próxima sección.
Un array y un valor (totalColumns) son suficientes para representar una matriz.
El conjunto puede tener (totalRows * totalColumns) variables, y la entrada en ( row, column ) puede corresponder a la variable array a Índice ((totalColumns * row) + column).
El código que aparece a continuación define una matriz, con 64 bits de punto flotante entradas.
Un array y un valor (totalColumns) son suficientes para representar una matriz.
El código que aparece a continuación tiene un abanico y un (totalColumns) valor que figura en una clase, sólo por conveniencia.
El código que aparece a continuación no pretende ser eficiente.
Una estructura (por ejemplo: “struct”, un tipo de valor) que representan a las matrices de dimensiones fijas particular (ejemplos: 2*2, 3*3, o 4*4) es probable que sea mucho más eficiente que el general (totalRows * totalColumns) clase se muestra aquí.
Aunque el código que aparece a continuación define con una matriz de punto flotante entradas, este artículo también hace uso de matrices con entradas entero.
El código que aparece a continuación se puede modificar fácilmente para aplicar las matrices con entradas entero.
using System;
using System.Collections.Generic;
using System.Text;
// Matrix with 64-bit floating-point entries:
public class MatrixF64
{
private int totalRows;
private int totalColumns;
private double[] entries;
// matrix[ row, column ] = entries[ (totalColumns * row) + column ]
public int Columns()
{
return (this.totalColumns);
}
public int Rows()
{
return (this.totalRows);
}
public MatrixF64()
{
this.totalRows = 0;
this.totalColumns = 0;
this.entries = null;
}
public MatrixF64( int rows, int columns, params double[] paramValues )
{
this.totalRows = 0;
this.totalColumns = 0;
this.entries = null;
if (rows <= 0)
{
return; // Specified an invalid row count.
}
if (columns <= 0)
{
return; // Specified an invalid column count.
}
int totalEntries = (rows * columns);
this.totalRows = rows;
this.totalColumns = columns;
this.entries = new double[totalEntries];
for (int k = 0; k < totalEntries; k++)
{
this.entries[k] = 0.0;
}
if (null == paramValues)
{
return; // No entries specified.
}
int totalParamValues = paramValues.Length;
if (totalParamValues != totalEntries)
{
// The number of specified entries does not match the
// total number of matrix entries.
return;
}
for (int k = 0; k < totalParamValues; k++)
{
this.entries[k] = paramValues[k];
}
}
public MatrixF64( MatrixF64 other )
{
this.totalRows = 0;
this.totalColumns = 0;
this.entries = null;
if (null == other)
{
return;
}
if
(
(null == other.entries)
|| (0 == other.totalRows)
|| (0 == other.totalColumns)
)
{
return;
}
int totalEntries = (other.totalRows * other.totalColumns);
if (other.entries.Length != totalEntries)
{
return;
}
this.totalRows = other.totalRows;
this.totalColumns = other.totalColumns;
this.entries = new double[totalEntries];
for (int k = 0; k < totalEntries; k++)
{
this.entries[k] = other.entries[k];
}
}
public double this[int row, int column]
{
get
{
if
(
(this.totalRows > 0)
&& (this.totalColumns > 0)
&& (null != this.entries)
&& (row >= 0)
&& (row < this.totalRows)
&& (column >= 0)
&& (column < this.totalColumns)
)
{
int k = (this.totalColumns * row) + column;
if ((k >= 0) && (k < this.entries.Length))
{
return (this.entries[k]);
}
}
return (0.0);
}
set
{
if
(
(this.totalRows > 0)
&& (this.totalColumns > 0)
&& (null != this.entries)
&& (row >= 0)
&& (row < this.totalRows)
&& (column >= 0)
&& (column < this.totalColumns)
)
{
int k = (this.totalColumns * row) + column;
if ((k >= 0) && (k < this.entries.Length))
{
this.entries[k] = value;
}
}
}
}
public void WriteLine( int precision )
{
if
(
(this.totalRows <= 0)
|| (this.totalColumns <= 0)
|| (null == this.entries)
)
{
Log.WriteLine( String.Empty + '[' + ' ' + ']' );
return;
}
// Determine the largest entry width in characters
// so that we can make all columns an equal width.
int largestEntryWidth = 1;
for (int i = 0; i < this.totalRows; i++)
{
for (int j = 0; j < this.totalColumns; j++)
{
// { index [,minwidth] [:typeCode[precision]] }
// (minwidth<0) means left-justify
String text = String.Format( String.Empty
+ '{' + '0' + ':' + 'g' + precision + '}', this[i,j] );
if (text.Length > largestEntryWidth)
{
largestEntryWidth = text.Length;
}
}
}
// Print each row of the matrix.
for (int i = 0; i < this.totalRows; i++)
{
Log.Write( '[' );
for (int j = 0; j < this.totalColumns; j++)
{
Log.Write( ' ' );
String text = String.Format( String.Empty
+ '{' + '0' + ',' + largestEntryWidth + ':' + 'g'
+ precision + '}', this[i,j] );
Log.Write( text );
if ((j + 1) < this.totalColumns)
{
Log.Write( ',' );
}
else
{
Log.Write( ' ' );
}
}
Log.WriteLine( ']' );
}
}
public void WriteLine()
{
const int defaultPrecision = 8;
this.WriteLine( defaultPrecision );
}
// . . .
public static void Test()
{
// A 5x3 matrix (5 rows x 3 columns) with
// 64-bit floating-point entries:
MatrixF64 m =
new MatrixF64
(
5, 3,
0.0, 1.0, 2.0, // row 0
3.0, 4.0, 5.0, // row 1
6.0, 7.0, 8.0, // row 2
9.0, 10.0, 11.0, // row 3
12.0, 13.0, 14.0 // row 4
);
m.WriteLine();
// [ 0, 1, 2 ] // row 0
// [ 3, 4, 5 ] // row 1
// [ 6, 7, 8 ] // row 2
// [ 9, 10, 11 ] // row 3
// [ 12, 13, 14 ] // row 4
// . . .
}
}
“cero matriz:” una matriz con todas las entradas iguales a cero.
public class MatrixF64
{
// . . .
public static MatrixF64 Zero( int rows, int columns )
{
MatrixF64 zero = new MatrixF64( rows, columns );
// The constructor above sets all entries to zero.
return (zero);
}
// . . .
public static void Test()
{
// . . .
// An 8x2 matrix (8 rows x 2 columns) with all 16
// 64-bit floating-point entries set to zero:
MatrixF64 z = MatrixF64.Zero( 8, 2 );
z.WriteLine();
// [ 0, 0 ] // row 0
// [ 0, 0 ] // row 1
// [ 0, 0 ] // row 2
// [ 0, 0 ] // row 3
// [ 0, 0 ] // row 4
// [ 0, 0 ] // row 5
// [ 0, 0 ] // row 6
// [ 0, 0 ] // row 7
// . . .
}
}
3.14 Matriz de suma, resta, multiplicación y
Matriz de suma, resta y multiplicación.
public class MatrixF64
{
// . . .
public static MatrixF64 operator +( MatrixF64 a, MatrixF64 b )
{
if ((null == a) || (null == b))
{
return (null);
}
MatrixF64 result = MatrixF64.Zero( a.totalRows, a.totalColumns );
if
(
(null == a.entries)
|| (null == b.entries)
|| (a.totalRows != b.totalRows)
|| (a.totalColumns != b.totalColumns)
)
{
return (result);
}
for (int i = 0; i < a.totalRows; i++)
{
for (int j = 0; j < a.totalColumns; j++)
{
result[i, j] = a[i, j] + b[i, j];
}
}
return (result);
}
public static MatrixF64 operator -( MatrixF64 a, MatrixF64 b )
{
if ((null == a) || (null == b))
{
return (null);
}
MatrixF64 result = MatrixF64.Zero( a.totalRows, a.totalColumns );
if
(
(null == a.entries)
|| (null == b.entries)
|| (a.totalRows != b.totalRows)
|| (a.totalColumns != b.totalColumns)
)
{
return (result);
}
for (int i = 0; i < a.totalRows; i++)
{
for (int j = 0; j < a.totalColumns; j++)
{
result[i, j] = a[i, j] - b[i, j];
}
}
return (result);
}
public static MatrixF64 operator -( MatrixF64 a )
{
if (null == a)
{
return (null);
}
MatrixF64 result = MatrixF64.Zero( a.totalRows, a.totalColumns );
if (null == a.entries)
{
return (result);
}
for (int i = 0; i < a.totalRows; i++)
{
for (int j = 0; j < a.totalColumns; j++)
{
result[i, j] = (-( a[i, j] ));
}
}
return (result);
}
public static MatrixF64 operator *( double scale, MatrixF64 a )
{
if (null == a)
{
return (null);
}
MatrixF64 result = MatrixF64.Zero( a.totalRows, a.totalColumns );
if (null == a.entries)
{
return (result);
}
for (int i = 0; i < a.totalRows; i++)
{
for (int j = 0; j < a.totalColumns; j++)
{
result[i, j] = scale * a[i, j];
}
}
return (result);
}
public static MatrixF64 operator *( MatrixF64 a, double scale )
{
if (null == a)
{
return (null);
}
MatrixF64 result = MatrixF64.Zero( a.totalRows, a.totalColumns );
if (null == a.entries)
{
return (result);
}
for (int i = 0; i < a.totalRows; i++)
{
for (int j = 0; j < a.totalColumns; j++)
{
result[i, j] = scale * a[i, j];
}
}
return (result);
}
public static MatrixF64 operator *( MatrixF64 a, MatrixF64 b )
{
if ((null == a) || (null == b))
{
return (null);
}
MatrixF64 result = MatrixF64.Zero( a.totalRows, b.totalColumns );
// The number of columns of the first matrix must be equal
// to the number of rows of the second matrix.
if
(
(null == a.entries)
|| (null == b.entries)
|| (a.totalColumns != b.totalRows)
)
{
return (result);
}
// Entry (i,j) of the result is equal to the dot product
// of row (i) of (a) and column (j) of (b).
for (int i = 0; i < a.totalRows; i++)
{
for (int j = 0; j < b.totalColumns; j++)
{
double dotProduct = 0.0;
for (int k = 0; k < a.totalColumns; k++)
{
dotProduct += (a[i, k] * b[k, j]);
}
result[i, j] = dotProduct;
}
}
return (result);
}
// . . .
public static void Test()
{
// . . .
// Examples of matrix addition, subtraction, and multiplication:
MatrixF64 a =
new MatrixF64
(
3, 5,
0.0, 1.0, 2.0, 3.0, 4.0, // row 0
5.0, 6.0, 7.0, 8.0, 9.0, // row 1
10.0, 11.0, 12.0, 13.0, 14.0 // row 2
);
MatrixF64 b =
new MatrixF64
(
3, 5,
4.0, -3.0, 2.0, -1.0, 0.0, // row 0
-9.0, 8.0, -7.0, 6.0, -5.0, // row 1
14.0, -13.0, 12.0, -11.0, 10.0 // row 2
);
MatrixF64 c =
new MatrixF64
(
5, 3,
0.0, 1.0, 2.0, // row 0
3.0, 4.0, 5.0, // row 1
6.0, 7.0, 8.0, // row 2
9.0, 10.0, 11.0, // row 3
12.0, 13.0, 14.0 // row 4
);
MatrixF64 result = new MatrixF64();
result = a + b;
result.WriteLine();
// [ 4, -2, 4, 2, 4 ]
// [ -4, 14, 0, 14, 4 ]
// [ 24, -2, 24, 2, 24 ]
result = a - b;
result.WriteLine();
// [ -4, 4, 0, 4, 4 ]
// [ 14, -2, 14, 2, 14 ]
// [ -4, 24, 0, 24, 4 ]
result = -a;
result.WriteLine();
// [ 0, -1, -2, -3, -4 ]
// [ -5, -6, -7, -8, -9 ]
// [ -10, -11, -12, -13, -14 ]
result = 3.0 * a;
result.WriteLine();
// [ 0, 3, 6, 9, 12 ]
// [ 15, 18, 21, 24, 27 ]
// [ 30, 33, 36, 39, 42 ]
result = a * c; // (3x5) * (5x3) = (3x3)
result.WriteLine();
// [ 90, 100, 110 ]
// [ 240, 275, 310 ]
// [ 390, 450, 510 ]
result = c * a; // (5x3) * (3x5) = (5x5)
result.WriteLine();
// [ 25, 28, 31, 34, 37 ]
// [ 70, 82, 94, 106, 118 ]
// [ 115, 136, 157, 178, 199 ]
// [ 160, 190, 220, 250, 280 ]
// [ 205, 244, 283, 322, 361 ]
// . . .
}
}
“matriz identidad:” Una matriz cuadrada (un total de filas es igual a total columnas) con entradas a (1) igual a la diagonal (entradas con un índice de fila igual a su índice de columna), y con todas las demás entradas a (0) igual.
Si una matriz cuadrada, (M), se multiplica por una “matriz de identidad,” (I), del mismo número de filas (y columnas), el producto es igual a (M).
(I) multiplicando por (M) también produce un producto igual a (M).
Por lo tanto, la “matriz identidad” es similar al número “1” para la multiplicación de números (escalares).
El código siguiente crea una matriz identidad con un número determinado de filas.
public class MatrixF64
{
// . . .
public static MatrixF64 Identity( int rows )
{
MatrixF64 identity = MatrixF64.Zero( rows, rows );
for (int i = 0; i < rows; i++)
{
identity[i, i] = 1.0;
}
return (identity);
}
// . . .
public static void Test()
{
// . . .
// Examples of multiplying by an identity matrix:
MatrixF64 identity3x3 = MatrixF64.Identity( 3 );
// [ 1, 0, 0 ]
// [ 0, 1, 0 ]
// [ 0, 0, 1 ]
MatrixF64 s3x3 =
new MatrixF64
(
3, 3,
0.0, 1.0, 2.0, // row 0
3.0, 4.0, 5.0, // row 1
6.0, 7.0, 8.0 // row 2
);
result = s3x3 * identity3x3;
result.WriteLine();
// [ 0, 1, 2 ]
// [ 3, 4, 5 ]
// [ 6, 7, 8 ]
// (the same as s3x3)
result = identity3x3 * s3x3;
result.WriteLine();
// [ 0, 1, 2 ]
// [ 3, 4, 5 ]
// [ 6, 7, 8 ]
// (the same as s3x3)
// . . .
}
}
3.15 Matrix: LU factoring; atrás de sustitución
“Matrix LU factoring:” un procedimiento que convierte una matriz cuadrada, (M), a dos nuevas matrices cuadradas, (L) y (U), que tengan el mismo tamaño que (M), de tal manera que (L) * (U) = (M), y de tal forma que la matriz (U) es “triangular superior” (todas las entradas por debajo de la diagonal son cero), y la matriz es “menor” (L) “triangular” (todas las entradas por encima de la diagonal son iguales a cero).
Matrix LU factoring puede ser utilizado como parte de un procedimiento para resolver un sistema de ecuaciones, o para encontrar la inversa de una matriz, o para encontrar el determinante de una matriz.
Matrix LU factoring es una alternativa a la eliminación “Gauss-Jordan procedimiento.
Gauss-Jordan eliminación requiere un sistema de ecuaciones (A)*(x)=(b), mientras que LU factoring sólo requiere una matriz (A).
Por otra parte, tras la determinación de la LU factoring de una matriz (A), es muy fácil determinar (x) dado ninguna (b).
El procedimiento para resolver el vector (x) en (A)*(x)=(b), habida cuenta de (b) y la LU factoring (L)*(U)=(A), implica un procedimiento que será llamado “de sustitución de volver” aquí.
El código que aparece a continuación incluye un procedimiento para calcular los factores de LU cualquier plaza, no singular matriz.
El código que aparece a continuación incluye un procedimiento para hacer back-LU sustitución.
Precaución: El código de computadora a continuación calcula la L y U factores de una fila de permutada versión de una determinada matriz.
Por otra parte, el L y U matriz de resultados se combinan en una sola matriz de salida.
public class MatrixF64
{
// . . .
public static void FindLUFactorsOfARowPermutedVersionOfAnOriginalMatrix
(
MatrixF64 originalMatrix,
ref MatrixF64 LUFactorsMatrix,
ref int[] rowPermutationOfOriginalMatrix,
ref bool rowExchangeParityIsOdd
)
{
// ''LU factoring'' involves factoring a square matrix (M) in to a
// lower-triangular matrix (L), and an upper-triangular matrix (U),
// such that (L)*(U) = (M).
//
// However, the specified matrix (originalMatrix) might not be
// optimal for direct LU factoring, due to the locations of extreme
// values in the matrix. Therefore, this function instead implicitly
// factors a row-permuted version of originalMatrix, where the
// permutation of rows is determined by the values in the specified
// matrix. The specific row permutation selected by this function is
// returned in an array of row indices (rowPermutationOfOriginalMatrix).
// The resulting (L) and (U) factors are merged in to a single
// output matrix (LUFactorsMatrix). Therefore:
//
// (L part of LUFactorsMatrix) * (U part of LUFactorsMatrix)
// = (originalMatrix with rows permuted according to
// rowPermutationOfOriginalMatrix).
//
// Although factoring a row-permuted version of originalMatrix
// makes it more difficult to interpret the results of this function,
// the resulting LU factoring is likely to be more accurate. In a
// sense, this function indicates (via rowPermutationOfOriginalMatrix)
// how to permute the rows of originalMatrix to be able to produce
// the most accurate LU factoring directly, and this function produces
// that factoring (via the LUFactorsMatrix output).
//
// When using the (L) and (U) matrices (merged in to LUFactorsMatrix)
// to solve a system of equations, it is necessary to permute the rows
// of the solution vector of the system of equations according to
// rowPermutation.
//
// The output matrix (LUFactorsMatrix) is a combination of
// the (L) and (U) matrices:
//
// To form the (L) matrix from (LUFactorsMatrix):
// assume (0) for entries above the diagonal,
// assume (1) for entries on the diagonal,
// and use (LUFactorsMatrix) entries below the diagonal.
//
// To form the (U) matrix from (LUFactorsMatrix):
// assume (0) for entries below the diagonal,
// and use (LUFactorsMatrix) entries on and above the diagonal.
//
// The output array (rowPermutationOfOriginalMatrix) indicates how the
// rows of the original matrix have implicitly (not actually) been
// permuted.
//
// The output boolean (rowExchangeParityIsOdd) indicates if the parity
// of the permutation of rows is odd.
//
// This implementation is partly based on pseudo-code appearing in
// ''Introduction to Algorithms'' (Cormen; Lieserson; Rivest;
// 24th printing, 2000; MIT Press; ISBN 0-262-03141-8).
// See the section named ''Overview of LUP decomposition'', in
// chapter 31 (''Matrix Operations''). The implementation here follows
// the ''LUP-Decomposition'' pseudo-code appearing in a section named
// ''Computing an LUP decomposition''.
//
// This implementation also uses ideas found in the implementation of
// ''ludcmp'' in the book ''Numerical recipes in C : The art of
// scientific computing'' (Press; Teukolsky; Vetterling; Flannery;
// second edition; 1996 reprinting; Cambridge University Press).
// The overall loop structure here resembles that of ''ludcmp''.
// The idea of determining and caching row scale factors in advance of
// finding pivots is used here. The method in ''ludcmp'' to avoid zero
// pivot values has been modified here to handle excessively-small
// pivot values, too. HOWEVER, the row permutation indices produced by
// the following implementation is a true permutation of
// {0,1,...,(totalRows-1)}, whereas ''ludcmp'' (and ''lubksb'')
// interpret their ''permutation indices'' as swaps to perform
// (roughly, for each (i), swap b[ i ] and b[ indx[i] ]). So, the
// following implementation of LU-factoring is not directly compatible
// with ''ludcmp'' or ''lubksb''. Also, the following procedure
// does not destroy the original matrix (unlike ''ludcmp'').
double singularMatrixIfMaxRowElementIsLessThanThis = (1.0e-19);
double forcePivotsToHaveAtLeastThisAbsoluteValue = (1.0e-19);
LUFactorsMatrix = null;
rowPermutationOfOriginalMatrix = null;
rowExchangeParityIsOdd = false;
if (null == originalMatrix)
{
return; // No matrix specified.
}
if
(
(originalMatrix.totalRows <= 0)
|| (originalMatrix.totalColumns <= 0)
|| (null == originalMatrix.entries)
)
{
return; // Matrix is empty.
}
if (originalMatrix.totalRows != originalMatrix.totalColumns)
{
return; // Matrix is not square.
}
// Duplicate the original matrix
LUFactorsMatrix = new MatrixF64( originalMatrix );
// (Lines 2-3 of LUP-Decomposition)
// Initialize the row permutation array to the identity
// permutation.
rowPermutationOfOriginalMatrix = new int[LUFactorsMatrix.totalRows];
for (int i = 0; i < LUFactorsMatrix.totalRows; i++)
{
rowPermutationOfOriginalMatrix[i] = i;
}
// For each row, determine the largest absolute value of
// the row elements, and use the reciprocal as the scale for
// that row.
VectorF64 rowScales = VectorF64.Zero( LUFactorsMatrix.totalRows );
for (int i = 0; i < LUFactorsMatrix.totalRows; i++)
{
double largestElementAbsoluteValueInRow = 0.0;
for (int j = 0; j < LUFactorsMatrix.totalColumns; j++)
{
double absoluteElementValue =
Math.Abs( LUFactorsMatrix[i, j] );
if (absoluteElementValue >
largestElementAbsoluteValueInRow)
{
largestElementAbsoluteValueInRow =
absoluteElementValue;
}
}
if (largestElementAbsoluteValueInRow <
singularMatrixIfMaxRowElementIsLessThanThis)
{
return; // Matrix is singular
}
rowScales[i] = (1.0 / largestElementAbsoluteValueInRow);
}
// (Lines 4-18 of LUP-Decomposition)
// Go through the columns of the matrix.
for (int j = 0; j < LUFactorsMatrix.totalColumns; j++)
{
// (Lines 17-18 of LUP-Decomposition)
// (or equation (2.3.12) of NR)
for (int i = 0; i < j; i++)
{
double sum = LUFactorsMatrix[i, j];
for (int k = 0; k < i; k++)
{
sum -= (LUFactorsMatrix[i, k] * LUFactorsMatrix[k, j]);
}
LUFactorsMatrix[i, j] = sum;
}
// Go through all remaining rows and find the best pivot.
double largestScaledSum = 0.0;
int rowIndexOfLargestScaledSum = j;
for (int i = j; i < LUFactorsMatrix.totalRows; i++)
{
// (equation (2.3.13) of NR)
double sum = LUFactorsMatrix[i, j];
for (int k = 0; k < j; k++)
{
sum -= (LUFactorsMatrix[i, k] * LUFactorsMatrix[k, j]);
}
LUFactorsMatrix[i, j] = sum;
double scaledSum = rowScales[i] * Math.Abs( sum );
if (scaledSum >= largestScaledSum)
{
largestScaledSum = scaledSum;
rowIndexOfLargestScaledSum = i;
}
}
// If indeed we found a better pivot, then exchange rows.
if (j != rowIndexOfLargestScaledSum)
{
// (Line 12 of LUP-Decomposition)
// Exchange the row permutation indices
int tempRowIndex = rowPermutationOfOriginalMatrix[j];
rowPermutationOfOriginalMatrix[j] =
rowPermutationOfOriginalMatrix[rowIndexOfLargestScaledSum];
rowPermutationOfOriginalMatrix[rowIndexOfLargestScaledSum] =
tempRowIndex;
// (Lines 13-14 of LUP-Decomposition)
// Exchange the elements of the rows
for (int k = 0; k < LUFactorsMatrix.totalColumns; k++)
{
double temp = LUFactorsMatrix[rowIndexOfLargestScaledSum, k];
LUFactorsMatrix[rowIndexOfLargestScaledSum, k] =
LUFactorsMatrix[j, k];
LUFactorsMatrix[j, k] = temp;
}
// Exchange the row scale factors
double scaleFactor =
rowScales[rowIndexOfLargestScaledSum];
rowScales[rowIndexOfLargestScaledSum] = rowScales[j];
rowScales[j] = scaleFactor;
// Invert the overall row exchange parity
rowExchangeParityIsOdd = (!(rowExchangeParityIsOdd));
}
// Force the pivot element to have at least a certain
// absolute value.
if (Math.Abs( LUFactorsMatrix[j, j] ) <
forcePivotsToHaveAtLeastThisAbsoluteValue)
{
if (LUFactorsMatrix[j, j] < 0.0)
{
LUFactorsMatrix[j, j] =
(-(forcePivotsToHaveAtLeastThisAbsoluteValue));
}
else
{
LUFactorsMatrix[j, j] =
forcePivotsToHaveAtLeastThisAbsoluteValue;
}
}
// If not the final column, then divide all column elements
// below the diagonal by the pivot element, matrixLU[j,j].
// (Lines 15-16 of LUP-Decomposition)
if (j != (LUFactorsMatrix.totalColumns - 1))
{
double reciprocalOfPivot = (1.0 / LUFactorsMatrix[j, j]);
for (int i = (j + 1); i < LUFactorsMatrix.totalRows; i++)
{
LUFactorsMatrix[i, j] *= reciprocalOfPivot;
}
}
}
}
public static void LUBacksubstitution
(
MatrixF64 LUFactorsMatrix,
int[] rowPermutationOfOriginalMatrix,
VectorF64 givenProductVector,
ref VectorF64 solutionVector
)
{
// This implementation is based on pseudo-code appearing in
// ''Introduction to Algorithms'' (Cormen; Lieserson; Rivest;
// 24th printing, 2000; MIT Press; ISBN 0-262-03141-8).
// See the section named ''Overview of LUP decomposition'', in
// chapter 31 (''Matrix Operations''). The implementation here
// follows the ''LUP-Solve'' pseudo-code appearing in a
// section named ''forward and back substitution''.
solutionVector = null;
if (null == LUFactorsMatrix)
{
return; // No matrix specified.
}
if (null == rowPermutationOfOriginalMatrix)
{
return; // No permutation specified.
}
if (null == givenProductVector)
{
return; // No product vector specified.
}
if
(
(null == LUFactorsMatrix.entries)
|| (LUFactorsMatrix.totalRows <= 0)
|| (LUFactorsMatrix.totalColumns <= 0)
)
{
return; // Matrix is empty.
}
if (LUFactorsMatrix.totalRows !=
LUFactorsMatrix.totalColumns)
{
return; // Matrix is not square.
}
if (givenProductVector.Dimensions() !=
LUFactorsMatrix.totalRows)
{
return; // Product vector size not equal to matrix row count.
}
// Copy the product vector in to the result.
solutionVector = new VectorF64( givenProductVector );
// (See LUP-Solve, lines 2-3)
for (int i = 0; i < LUFactorsMatrix.totalRows; i++)
{
double sum =
givenProductVector[rowPermutationOfOriginalMatrix[i]];
for (int j = 0; j < i; j++)
{
sum -= (LUFactorsMatrix[i, j] * solutionVector[j]);
}
solutionVector[i] = sum;
}
// (See LUP-Solve, lines 4-5)
for (int i = (LUFactorsMatrix.totalRows - 1); i >= 0; i--)
{
double sum = solutionVector[i];
for
(
int j = (i + 1);
j < LUFactorsMatrix.totalColumns;
j++
)
{
sum -= (LUFactorsMatrix[i, j] * solutionVector[j]);
}
double diagonalElement = LUFactorsMatrix[i, i];
if (diagonalElement != 0.0)
{
solutionVector[i] = (sum / diagonalElement);
}
}
}
// . . .
public static void Test()
{
// . . .
// Example of LU factoring and back-substitution:
MatrixF64 a3x3 =
new MatrixF64
(
3, 3,
0.0, 1.0, 2.0, // row 0
3.0, 4.0, 5.0, // row 1
6.0, 7.0, 8.0 // row 2
);
int[] rowPermutationOfOriginalMatrix = null;
bool rowInterchangeParityIsOdd = false;
MatrixF64 a3x3LU = null;
MatrixF64.FindLUFactorsOfARowPermutedVersionOfAnOriginalMatrix
(
a3x3,
ref a3x3LU,
ref rowPermutationOfOriginalMatrix,
ref rowInterchangeParityIsOdd
);
a3x3LU.WriteLine();
// [ 6, 7, 8 ]
// [ 0, 1, 2 ]
// [ 0.5, 0.5, 1e-19 ]
// The matrix above is a combination of the L and U
// matrices. The L matrix is 0 above the diagonal,
// 1 on the diagonal, and the shown values below the
// diagonal. The U matrix is 0 below the diagonal,
// and the shown values on and above the diagonal.
// This LU factoring is of a ROW-PERMUTED version
// of the original matrix.
// The following L and U matrices are manually formed
// by inspecting the LU factoring matrix above.
MatrixF64 L3x3 =
new MatrixF64
(
3, 3,
1.0, 0.0, 0.0, // row 0
0.0, 1.0, 0.0, // row 1
0.5, 0.5, 1.0 // row 2
);
MatrixF64 U3x3 =
new MatrixF64
(
3, 3,
6.0, 7.0, 8.0, // row 0
0.0, 1.0, 2.0, // row 1
0.0, 0.0, 0.0 // row 2
);
// The product (L)*(U) produces a ROW-PERMUTED version
// of the original matrix:
result = L3x3 * U3x3;
result.WriteLine( 4 );
// [ 6, 7, 8 ]
// [ 0, 1, 2 ]
// [ 3, 4, 5 ]
VectorF64 b3x1 = new VectorF64( 1.0, 2.0, 3.0 );
VectorF64 x3x1 = null;
MatrixF64.LUBacksubstitution
(
a3x3LU,
rowPermutationOfOriginalMatrix,
b3x1,
ref x3x1
);
// Solution vector
x3x1.WriteLine();
// ( -0.66666667, 1, 0 )
// Both the given b3x1 and this x3x1 are relative to the
// original matrix, a3x3. The a3x3LU for a ROW-PERMUTED
// version of a3x3, but LUBacksubstitution accepts a
// non-permuted product vector (e.g., b3x1) and produces
// a non-permuted solution vector (e.g., x3x1).
// Check solution by multiplying the original matrix a3x3
// by the solution vector x3x1, to get a product vector.
VectorF64 checkb3x1 = a3x3 * x3x1;
checkb3x1.WriteLine();
// ( 1, 2, 3 )
// This matches the b3x1 vector we specified above, so
// the solution is valid.
// . . .
}
}
3.16 Matrix: determinante
“determinante:” Una función que convierte cualquier plaza (n * n) matriz (A) a un número, det(A), tal que:
(1) Si una matriz (B) resultados de intercambio de dos filas, o dos columnas, una matriz de (A), entonces det(B) = (-(det(A)));
(2) Si una matriz (B) multiplicando los resultados de cualquier fila, columna o cualquiera, de una matriz (A), por una serie (c), entonces det(B) = c * det(A);
(3) Si una matriz (B) resultados de la adición de un múltiplo de una fila a otra fila, o la adición de un múltiplo de una columna a otra columna, entonces det(B) = det(A).
(4) Si (A) superior es una matriz triangular (A[i,j]=0 para todos (i>j)) o una menor-matriz triangular (A[i,j]=0 para todos (i<j)), entonces det(A) = (A[0,0] * A[1,1] * ... * A[n-1,n-1]);
(Por ejemplo, el factor determinante de una matriz de identidad es uno; det(I) = 1.)
Reglas (1), (2), y (3), pueden utilizarse, en un proceso llamado “Gaussian eliminación, para convertir cualquier matriz cuadrada a una matriz triangular.
Artículo (4) puede ser utilizado para calcular el determinante de una matriz triangular.
La definición formal de un “determinante:”
Que M[n](K) denotar el conjunto de todos los (n * n) matrices en el campo (K).
El “factor determinante” es la única función (F) tal que F:M[n](K) --> K con dos atributos:
(1) F es multilinear alternando con respecto a las columnas (o filas), y,
(2) F( I ) = 1.
A “multilinear” función (D) de un (n * n) matriz (A) puede escribirse como: D(A) = Sum( A[0,k[0]] * A[1,k[1]] * ... * A[n-1,k[n-1]] * D( e[k[0]], e[k[1]], ..., e[k[n-1]] ) ), donde la suma se toma sobre todas las combinaciones (n^n) de 0 <= k[i] <= (n-1), y donde e[j] representa (j) fila de la matriz identidad.
Esta función puede ser caracterizado por los valores de D( e[k[0]], e[k[1]], ..., e[k[n-1]] ).
Una “alternancia multilinear” función es una función multilinear, (D), de tal manera que D( ..., e[i], ..., e[i], ... ) = 0, y D( ..., e[i], e[j], ... ) = (-( D( ..., e[j], e[i], ... ) )).
De este modo, el canje de las columnas (o filas) cambia el signo de la función.
Por otra parte, el término D( e[k[0]], ..., e[k[n-1]] ) es cero si la combinación de valores k[i] no involucrar a todos los valores únicos, si el conjunto de valores no es una permutación de los números {0,1,...,(n-1)}.
Dejar D( I ) = 1 (D( e[0], e[1], ..., e[n-1] ) = 1), en combinación con la alternancia multilinear atributos descritos anteriormente, significa que todos los valores de la alternancia multilinear función D( e[k[0]], e[k[1]], ..., e[k[n-1]] ) puede determinarse.
La función será distinta de cero sólo si el conjunto de valores k[i] es una permutación de {0,1,...,(n-1)}, y tendrá una magnitud de uno de ellos si esto es distinto de cero, y el signo será el número de intercambios necesarios para convertir la permutación de la permutación identidad.
El valor de un determinante de una (n * n) matriz (A) puede ser calculado sumando los productos de la forma (sgn(p) * A[0,p[0]] * A[1,p[1]] * ... * A[(n-1),p[n-1]]) para todos y cada uno de (p) permutación de los números {0,1,2,...,(n-1)}.
El término denota (sgn(p)) la “firma” (o swap contar) de permutación (p), donde (sgn(p)) es igual a (+1) si (p) es “una permutación,” y es igual a (-1) si (p) es una “permutación impar.”
Debido a que hay (n!) permutaciones de los números {0,1,2,...,(n-1)}, esta fórmula ha (n!) summands.
Otros procedimientos para calcular el factor determinante (ejemplos: las relacionadas con Gaussian eliminación, o LU factoring, o la expansión de los menores, etc) calcular implícitamente una jerarquía intermedia de las cantidades que permiten a los procedimientos para calcular el factor determinante en un orden de (n^3) pasos.
Si el determinante de una matriz es cero, la matriz no tiene una inversa.
Si el determinante de una matriz es distinto de cero, la matriz tiene una inversa.
Si una matriz que representa los coeficientes de un sistema lineal de ecuaciones, y el factor determinante de que la matriz es cero, entonces el sistema de ecuaciones no tiene una solución única.
Si una matriz que representa los coeficientes de un sistema lineal de ecuaciones, y el factor determinante de la matriz que no es cero, entonces el sistema de ecuaciones tiene una solución única.
El código que aparece a continuación incluye los procedimientos para calcular el factor determinante de cualquier matriz cuadrada.
Factores determinantes para 1*1, 2*2, 3*3, y 4*4 matrices son calculadas por las fórmulas explícitas.
Uso de fórmulas explícitas para calcular los factores determinantes para las pequeñas matrices es probable que computar los resultados más rápido que utilizando un procedimiento general para calcular los factores determinantes.
El procedimiento para el caso general LU utiliza el factoring, pero hay muchos otros métodos de cálculo de los factores determinantes.
public class MatrixF64
{
// . . .
private double Determinant1x1()
{
if
(
(this.totalRows != 1)
|| (this.totalColumns != 1)
|| (null == this.entries)
)
{
return (0.0);
}
// The determinant of a 1x1 matrix is simply
// the value of the single element.
return ( this[0,0] );
}
private double Determinant2x2()
{
if
(
(this.totalRows != 2)
|| (this.totalColumns != 2)
|| (null == this.entries)
)
{
return (0.0);
}
return
(
+ this[0,0] * this[1,1]
- this[0,1] * this[1,0]
);
}
private double Determinant3x3()
{
if
(
(this.totalRows != 3)
|| (this.totalColumns != 3)
|| (null == this.entries)
)
{
return (0.0);
}
return
(
+ this[0,0] * this[1,1] * this[2,2]
+ this[0,1] * this[1,2] * this[2,0]
+ this[0,2] * this[1,0] * this[2,1]
- this[0,0] * this[1,2] * this[2,1]
- this[0,1] * this[1,0] * this[2,2]
- this[0,2] * this[1,1] * this[2,0]
);
}
private double Determinant4x4()
{
if
(
(this.totalRows != 4)
|| (this.totalColumns != 4)
|| (null == this.entries)
)
{
return (0.0);
}
return
(
+ this[0, 0] * this[1, 1] * this[2, 2] * this[3, 3]
+ this[0, 0] * this[1, 2] * this[2, 3] * this[3, 1]
+ this[0, 0] * this[1, 3] * this[2, 1] * this[3, 2]
+ this[0, 1] * this[1, 0] * this[2, 3] * this[3, 2]
+ this[0, 1] * this[1, 2] * this[2, 0] * this[3, 3]
+ this[0, 1] * this[1, 3] * this[2, 2] * this[3, 0]
+ this[0, 2] * this[1, 0] * this[2, 1] * this[3, 3]
+ this[0, 2] * this[1, 1] * this[2, 3] * this[3, 0]
+ this[0, 2] * this[1, 3] * this[2, 0] * this[3, 1]
+ this[0, 3] * this[1, 0] * this[2, 2] * this[3, 1]
+ this[0, 3] * this[1, 1] * this[2, 0] * this[3, 2]
+ this[0, 3] * this[1, 2] * this[2, 1] * this[3, 0]
- this[0, 0] * this[1, 1] * this[2, 3] * this[3, 2]
- this[0, 0] * this[1, 2] * this[2, 1] * this[3, 3]
- this[0, 0] * this[1, 3] * this[2, 2] * this[3, 1]
- this[0, 1] * this[1, 0] * this[2, 2] * this[3, 3]
- this[0, 1] * this[1, 2] * this[2, 3] * this[3, 0]
- this[0, 1] * this[1, 3] * this[2, 0] * this[3, 2]
- this[0, 2] * this[1, 0] * this[2, 3] * this[3, 1]
- this[0, 2] * this[1, 1] * this[2, 0] * this[3, 3]
- this[0, 2] * this[1, 3] * this[2, 1] * this[3, 0]
- this[0, 3] * this[1, 0] * this[2, 1] * this[3, 2]
- this[0, 3] * this[1, 1] * this[2, 2] * this[3, 0]
- this[0, 3] * this[1, 2] * this[2, 0] * this[3, 1]
);
}
public double Determinant()
{
if
(
(this.totalRows <= 0)
|| (this.totalColumns <= 0)
|| (null == this.entries)
)
{
return (0.0); // Matrix is empty.
}
if (this.totalRows != this.totalColumns)
{
// Matrix is not square
return (0.0);
}
// Simple cases
if (1 == this.totalRows) { return (this.Determinant1x1()); }
if (2 == this.totalRows) { return (this.Determinant2x2()); }
if (3 == this.totalRows) { return (this.Determinant3x3()); }
if (4 == this.totalRows) { return (this.Determinant4x4()); }
int[] rowPermutationOfOriginalMatrix = null;
bool rowExchangeParityIsOdd = false;
MatrixF64 originalMatrix = new MatrixF64( this );
MatrixF64 LUFactorsMatrix = null;
MatrixF64.FindLUFactorsOfARowPermutedVersionOfAnOriginalMatrix
(
originalMatrix,
ref LUFactorsMatrix,
ref rowPermutationOfOriginalMatrix,
ref rowExchangeParityIsOdd
);
if ((null == LUFactorsMatrix)
|| (null == rowPermutationOfOriginalMatrix))
{
return (0.0);
}
double determinant = 1.0;
if (true == rowExchangeParityIsOdd)
{
determinant = (-1.0);
}
for (int j = 0; j < LUFactorsMatrix.totalColumns; j++)
{
determinant *= LUFactorsMatrix[j, j];
}
return (determinant);
}
// . . .
public static void Test()
{
// . . .
// Examples of matrix determinants:
MatrixF64 d2x2 =
new MatrixF64
(
2, 2,
2.0, 5.0, // row 0
-4.0, -7.0 // row 1
);
double detd2x2 = d2x2.Determinant();
Log.WriteLine( detd2x2 );
// detd2x2 = 6
MatrixF64 d3x3 =
new MatrixF64
(
3, 3,
2.0, 5.0, 7.0, // row 0
-4.0, -1.0, 6.0, // row 1
9.0, 8.0, 3.0 // row 2
);
double detd3x3 = d3x3.Determinant();
Log.WriteLine( detd3x3 );
// detd3x3 = 67
MatrixF64 d4x4 =
new MatrixF64
(
4, 4,
7.0, -5.0, 2.0, 4.0, // row 0
3.0, 2.0, 6.0, 3.0, // row 1
-9.0, 8.0, -3.0, 2.0, // row 2
5.0, 3.0, 2.0, 5.0 // row 3
);
double detd4x4 = d4x4.Determinant();
Log.WriteLine( detd4x4 );
// detd4x4 = 1457
// . . .
}
}
3.17 (d) dimensiones vectores: producto cruzado
“producto cruzado:” Una función asociada con un (d) dimensiones del espacio que convierte una lista ordenada de (d-1) (d) dimensiones vectores, { A(0), A(1), ..., A(d-2) }, a un (d)-dimensional vector, de manera que la función tiene los atributos siguientes:
(1) Si A(i)=0 (vector cero), para cualquier 0 <= i <= (d-2), entonces Cross( A(0), A(1), ..., A(d-2) )=0 (vector cero).
La cruz producto de una lista ordenada de los vectores es cero si todo vector en esa lista es cero;
(2) Si cualquier par de vectores, A(i) y A(j), para cualquier (i!=j) con 0 <= i <= (d-2) y 0 <= j <= (d-2), (Abs(Dot(A(i),A(j))) = Length(A(i)) * Length(A(j))) son paralelos, entonces Cross( A(0), A(1), ..., A(d-2) )=0 (vector cero).
La cruz producto de una lista ordenada de los vectores es cero si dos vectores en esa lista son paralelos;
(3) Si cada vector A(i) es distinto de cero, y si no es A(i) paralelo con A(j) para todos (i!=j) con 0 <= i <= (d-2) y 0 <= j <= (d-2), entonces B = Cross( A(0), A(1), ..., A(d-2) ) es tal que Dot(A(i),B)=0 para todos 0 <= i <= (d-2).
Si el producto cruzado de una lista ordenada de los vectores no es cero, entonces el producto cruzado resultado es perpendicular a cada vector en la lista ordenada de los vectores;
(4) Si (c) es un valor numérico, entonces Cross( (c * A(0)), A(1), ..., A(d-2) ) = Cross( A(0), (c * A(1)), ..., A(d-2) ) = ... = Cross( A(0), A(1), ..., (c * A(d-2)) ) = c * Cross( A(0), A(1), ..., A(d-2) ).
Dada una lista ordenada de los vectores, el producto cruzado de la lista ordenada de vectores con un vector multiplicado por un determinado valor numérico es equivalente a multiplicar el producto cruzado de la lista ordenada de los vectores por el valor numérico;
(5) Habida cuenta de dos números (i) y (j), donde (i!=j) y 0 <= i <= (d-2) y 0 <= j <= (d-2), y una lista ordenada de (d-1) (d) dimensiones vectores, { A(0), A(1), ..., A(d-2) }, y otra lista ordenada de (d-1) (d) dimensiones vectores, { B(0), B(1), ..., B(d-2) }, de tal manera que para todos los B(k) = A(k) con ((k != i) && (k != j)) 0 <= k <= (d-2), y de tal manera que B(i) = A(j) y B(j) = A(i), entonces Cross( A(0), A(1), ..., A(d-2) ) = (-1) * Cross( B(0), B(1), ..., B(d-2) ).
Si el producto cruzado de una lista ordenada de los vectores no es cero, a continuación, que cruzan producto tiene el signo de la cruz producto de esa misma lista ordenada de los vectores, con la excepción de exactamente dos vectores que se intercambian (cambiados).
En vista de los vectores base e(k), para 0 <= k <= (d-1), (d) de dimensiones espaciales, la cruz producto de una lista ordenada de (d-1) vectores pueden ser calculadas por la informática el factor determinante de la matriz (d*d) a continuación:
e(0), e(1), ..., e(d-1),
A( 0 )[0], A( 0 )[1], ..., A( 0 )[d-1],
A( 1 )[0], A( 1 )[1], ..., A( 1 )[d-1],
...,
A(d-2)[0], A(d-2)[1], ..., A(d-2)[d-1]
El factor determinante de la matriz que no será un número, sino que ser un vector, donde se coeficientes numéricos a la base (d) vectores.
Durante dos dimensiones espaciales, y el vector A = ( ax, ay ):
Cross( A )
= Determinant
(
e(x), e(y),
ax, ay
)
= ay * e(x)
- ax * e(y)
= ( ay, -ax )
El resultado es el vector perpendicular al vector de entrada única.
Si el avión x-y se considera que tales aumentos x la derecha en una dirección y y aumenta en una dirección hacia arriba, entonces el producto cruzado resultado es un cuarto de vuelta en “sentido horario” en relación con la entrada de vectores.
Ejemplos de productos en cruz de dos dimensiones del espacio:
Cross( e(x) ) = (-1) * e(y)
Cross( e(y) ) = e(x)
Por espacio tridimensional, los vectores y A = ( ax, ay, az ) y B = ( bx, by, bz ):
Cross( A, B )
= Determinant
(
e(x), e(y), e(z),
ax, ay, az,
bx, by, bz
)
= (ay*bz - az*by) * e(x)
+ (az*bx - ax*bz) * e(y)
+ (ax*by - ay*bx) * e(z)
= ( (ay*bz - az*by), (az*bx - ax*bz), (ax*by - ay*bx) )
Ejemplos de productos cruz en un espacio tridimensional:
Cross( e(x), e(y) ) = e(z)
Cross( e(y), e(x) ) = (-1) * e(z)
Cross( e(y), e(z) ) = e(x)
Cross( e(z), e(y) ) = (-1) * e(x)
Cross( e(z), e(x) ) = e(y)
Cross( e(x), e(z) ) = (-1) * e(y)
Durante cuatro dimensiones espaciales, y los vectores A = ( ax, ay, az, aw ), B = ( bx, by, bz, bw ), y C = ( cx, cy, cz, cw ):
Cross( A, B, C )
= Determinant
(
e(x), e(y), e(z), e(w),
ax, ay, az, aw,
bx, by, bz, bw,
cx, cy, cz, cw
)
= ( ay*bz*cw - ay*bw*cz + az*bw*cy - az*by*cw + aw*by*cz - aw*bz*cy) * e(x)
+ (-ax*bz*cw + ax*bw*cz - az*bw*cx + az*bx*cw - aw*bx*cz + aw*bz*cx) * e(y)
+ ( ax*by*cw - ax*bw*cy + ay*bw*cx - ay*bx*cw + aw*bx*cy - aw*by*cx) * e(z)
+ (-ax*by*cz + ax*bz*cy - ay*bz*cx + ay*bx*cz - az*bx*cy + az*by*cx) * e(w)
= ( ( ay*bz*cw - ay*bw*cz + az*bw*cy - az*by*cw + aw*by*cz - aw*bz*cy),
(-ax*bz*cw + ax*bw*cz - az*bw*cx + az*bx*cw - aw*bx*cz + aw*bz*cx),
( ax*by*cw - ax*bw*cy + ay*bw*cx - ay*bx*cw + aw*bx*cy - aw*by*cx),
(-ax*by*cz + ax*bz*cy - ay*bz*cx + ay*bx*cz - az*bx*cy + az*by*cx) )
Ejemplos de productos en cruz de cuatro dimensiones del espacio:
Cross( e(x), e(y), e(z) ) = (-1) * e(w)
Cross( e(x), e(z), e(y) ) = e(w)
Cross( e(y), e(x), e(z) ) = e(w)
Cross( e(y), e(z), e(x) ) = (-1) * e(w)
Cross( e(z), e(x), e(y) ) = (-1) * e(w)
Cross( e(z), e(y), e(x) ) = e(w)
Cross( e(x), e(y), e(w) ) = e(z)
Cross( e(x), e(w), e(y) ) = (-1) * e(z)
Cross( e(y), e(x), e(w) ) = (-1) * e(z)
Cross( e(y), e(w), e(x) ) = e(z)
Cross( e(w), e(x), e(y) ) = e(z)
Cross( e(w), e(y), e(x) ) = (-1) * e(z)
Cross( e(x), e(z), e(w) ) = (-1) * e(y)
Cross( e(x), e(w), e(z) ) = e(y)
Cross( e(z), e(x), e(w) ) = e(y)
Cross( e(z), e(w), e(x) ) = (-1) * e(y)
Cross( e(w), e(x), e(z) ) = (-1) * e(y)
Cross( e(w), e(z), e(x) ) = e(y)
Cross( e(y), e(z), e(w) ) = e(x)
Cross( e(y), e(w), e(z) ) = (-1) * e(x)
Cross( e(z), e(y), e(w) ) = (-1) * e(x)
Cross( e(z), e(w), e(y) ) = e(x)
Cross( e(w), e(y), e(z) ) = e(x)
Cross( e(w), e(z), e(y) ) = (-1) * e(x)
El código que aparece a continuación calcula la cruz-producto de dos o más vectores.
La prueba ejemplos muestran que la falta de resultado cero vectores son perpendiculares a todos los vectores de entrada.
public class MatrixF64
{
// . . .
public static VectorF64 Cross( params VectorF64[] parameterArray )
{
if (null == parameterArray)
{
return(new VectorF64());
}
int totalVectors = parameterArray.Length;
// The number of dimensions of the space must be one larger
// than the number of specified vectors.
int dimensions = (1 + totalVectors);
if (dimensions < 2)
{
// The cross product does not exist for space with
// fewer than two dimensions.
return (new VectorF64());
}
// All specified vectors must exist and must have a number
// of components equal to the number of dimensions of the space.
for (int k = 0; k < totalVectors; k++)
{
if (null == parameterArray[k])
{
// Vector missing.
return (new VectorF64());
}
else if (null == parameterArray[k].components)
{
// Vector is empty.
return (new VectorF64());
}
else if (parameterArray[k].components.Length != dimensions)
{
// Vector has the wrong number of components.
return (new VectorF64());
}
}
// Form a (d*d) matrix with the (d-1) supplied vectors
// forming the final (d-1) rows of the matrix.
MatrixF64 matrix = MatrixF64.Zero( dimensions, dimensions );
for (int i = 1; i < dimensions; i++)
{
VectorF64 rowVector = parameterArray[ i - 1 ];
for (int j = 0; j < dimensions; j++)
{
matrix[i, j] = rowVector[j];
}
}
VectorF64 result = VectorF64.Zero( dimensions );
// For each component of the result vector, set the first
// row of the matrix to a basis vector, compute the matrix
// determinant, and use the resulting value as the
// component of the result vector. This is inefficient,
// due to the determinant being computed (d) times.
// Clearly it would be better to have explicit formulas
// for 2-dimensional, 3-dimensional, and 4-dimensional
// cross-products.
for (int i = 0; i < dimensions; i++)
{
// Set the first row of the matrix to be the basis
// vector for the coordinate axis with index (i).
for (int j = 0; j < dimensions; j++)
{
matrix[0, j] = 0.0;
}
matrix[0, i] = 1.0;
// Compute the determinant of the modified matrix.
double determinant = matrix.Determinant();
// Set component (i) of the result to the determinant
// that was computed.
result[i] = determinant;
}
return (result);
}
// . . .
public static void Test()
{
// . . .
// Examples of cross products:
// Three-dimensional example:
VectorF64 vcp1x3a = new VectorF64( 1.0, 2.0, 3.0 );
VectorF64 vcp1x3b = new VectorF64( 3.0, 1.0, 2.0 );
VectorF64 vcp1x3 = VectorF64.Cross( vcp1x3a, vcp1x3b );
vcp1x3.WriteLine(); // ( 1, 7, -5 )
// Verify that the result is perpendicular to the original
// vectors:
Log.WriteLine( VectorF64.Perpendicular( vcp1x3, vcp1x3a ) );
Log.WriteLine( VectorF64.Perpendicular( vcp1x3, vcp1x3b ) );
// True, True
// Four-dimensional example:
VectorF64 vcp1x4a = new VectorF64( 1.0, 2.0, 3.0, 4.0 );
VectorF64 vcp1x4b = new VectorF64( 4.0, 1.0, 2.0, 3.0 );
VectorF64 vcp1x4c = new VectorF64( 3.0, 4.0, 1.0, 2.0 );
VectorF64 vcp1x4 = VectorF64.Cross( vcp1x4a, vcp1x4b, vcp1x4c );
vcp1x4.WriteLine(); // ( 4, 4, 44, -36 )
// Verify that the result is perpendicular to the original
// vectors:
Log.WriteLine( VectorF64.Perpendicular( vcp1x4, vcp1x4a ) );
Log.WriteLine( VectorF64.Perpendicular( vcp1x4, vcp1x4b ) );
Log.WriteLine( VectorF64.Perpendicular( vcp1x4, vcp1x4c ) );
// True, True, True
// . . .
}
}
3.18 Matrix: inversa
“matriz inversa:” Una matriz cuadrada en relación con otra matriz cuadrada del mismo tamaño tal que el producto de las dos matrices es la “matriz identidad.”
Si una matriz tiene una matriz asociada inversa, sólo hay uno de esos matriz inversa.
Si el determinante de una matriz es cero, entonces la matriz es el nombre “singular,” y la matriz no tiene asociada una matriz inversa.
El código que aparece a continuación calcula la inversa de cualquier plaza, no singular matriz.
Inversas para 1*1, 2*2, 3*3, y 4*4 matrices son calculadas por las fórmulas explícitas.
Uso de fórmulas explícitas podría ser más rápido que utilizando un procedimiento general.
Las fórmulas explícitas también son fiables.
El procedimiento para el caso general LU utiliza el factoring, pero hay muchos otros métodos de computación de matrices inversas.
public class MatrixF64
{
// . . .
private MatrixF64 Inverse1x1()
{
if
(
(this.totalRows != 1)
|| (this.totalColumns != 1)
|| (null == this.entries)
)
{
return (MatrixF64.Zero( 1, 1 )); // Matrix is empty.
}
if (0.0 == this[0,0])
{
return ( MatrixF64.Zero( 1, 1 ) ); // Matrix has no inverse.
}
MatrixF64 inverse = MatrixF64.Zero( 1, 1 );
// The inverse of a 1x1 matrix is simply the
// reciprocal of the single entry.
inverse[0,0] = (1.0 / this[0,0]);
return (inverse);
}
private MatrixF64 Inverse2x2()
{
if
(
(this.totalRows != 2)
|| (this.totalColumns != 2)
|| (null == this.entries)
)
{
return (MatrixF64.Zero( 2, 2 )); // Matrix is empty.
}
double determinant = this.Determinant2x2();
if (0.0 == determinant)
{
return (MatrixF64.Zero( 2, 2 )); // Matrix has no inverse.
}
MatrixF64 inverse = MatrixF64.Zero( 2, 2 );
inverse[0,0] = this[1,1];
inverse[0,1] = (-(this[0,1]));
inverse[1,0] = (-(this[1,0]));
inverse[1,1] = this[0,0];
double factor = (1.0 / determinant);
for ( int i = 0; i < inverse.totalRows; i++ )
{
for ( int j = 0; j < inverse.totalColumns; j++ )
{
inverse[i,j] *= factor;
}
}
return (inverse);
}
private MatrixF64 Inverse3x3()
{
if
(
(this.totalRows != 3)
|| (this.totalColumns != 3)
|| (null == this.entries)
)
{
return (MatrixF64.Zero( 3, 3 )); // Matrix is empty.
}
double determinant = this.Determinant3x3();
if (0.0 == determinant)
{
return (MatrixF64.Zero( 3, 3 )); // Matrix has no inverse.
}
MatrixF64 inverse = MatrixF64.Zero( 3, 3 );
inverse[0, 0] =
this[1, 1] * this[2, 2]
- this[1, 2] * this[2, 1];
inverse[0, 1] =
this[0, 2] * this[2, 1]
- this[0, 1] * this[2, 2];
inverse[0, 2] =
this[0, 1] * this[1, 2]
- this[0, 2] * this[1, 1];
inverse[1, 0] =
this[1, 2] * this[2, 0]
- this[1, 0] * this[2, 2];
inverse[1, 1] =
this[0, 0] * this[2, 2]
- this[0, 2] * this[2, 0];
inverse[1, 2] =
this[0, 2] * this[1, 0]
- this[0, 0] * this[1, 2];
inverse[2, 0] =
this[1, 0] * this[2, 1]
- this[1, 1] * this[2, 0];
inverse[2, 1] =
this[0, 1] * this[2, 0]
- this[0, 0] * this[2, 1];
inverse[2, 2] =
this[0, 0] * this[1, 1]
- this[0, 1] * this[1, 0];
double factor = (1.0 / determinant);
for (int i = 0; i < inverse.totalRows; i++)
{
for (int j = 0; j < inverse.totalColumns; j++)
{
inverse[i, j] *= factor;
}
}
return (inverse);
}
private MatrixF64 Inverse4x4()
{
if
(
(this.totalRows != 4)
|| (this.totalColumns != 4)
|| (null == this.entries)
)
{
return (MatrixF64.Zero( 4, 4 )); // Matrix is empty.
}
double determinant = this.Determinant4x4();
if (0.0 == determinant)
{
return (MatrixF64.Zero( 4, 4 )); // Matrix has no inverse.
}
MatrixF64 inverse = MatrixF64.Zero( 4, 4 );
inverse[0,0] =
+ this[1,1] * this[2,2] * this[3,3]
+ this[1,2] * this[2,3] * this[3,1]
+ this[1,3] * this[2,1] * this[3,2]
- this[1,1] * this[2,3] * this[3,2]
- this[1,2] * this[2,1] * this[3,3]
- this[1,3] * this[2,2] * this[3,1];
inverse[0,1] =
+ this[0,1] * this[2,3] * this[3,2]
+ this[0,2] * this[2,1] * this[3,3]
+ this[0,3] * this[2,2] * this[3,1]
- this[0,1] * this[2,2] * this[3,3]
- this[0,2] * this[2,3] * this[3,1]
- this[0,3] * this[2,1] * this[3,2];
inverse[0,2] =
+ this[0,1] * this[1,2] * this[3,3]
+ this[0,2] * this[1,3] * this[3,1]
+ this[0,3] * this[1,1] * this[3,2]
- this[0,1] * this[1,3] * this[3,2]
- this[0,2] * this[1,1] * this[3,3]
- this[0,3] * this[1,2] * this[3,1];
inverse[0,3] =
+ this[0,1] * this[1,3] * this[2,2]
+ this[0,2] * this[1,1] * this[2,3]
+ this[0,3] * this[1,2] * this[2,1]
- this[0,1] * this[1,2] * this[2,3]
- this[0,2] * this[1,3] * this[2,1]
- this[0,3] * this[1,1] * this[2,2];
inverse[1,0] =
+ this[1,0] * this[2,3] * this[3,2]
+ this[1,2] * this[2,0] * this[3,3]
+ this[1,3] * this[2,2] * this[3,0]
- this[1,0] * this[2,2] * this[3,3]
- this[1,2] * this[2,3] * this[3,0]
- this[1,3] * this[2,0] * this[3,2];
inverse[1,1] =
+ this[0,0] * this[2,2] * this[3,3]
+ this[0,2] * this[2,3] * this[3,0]
+ this[0,3] * this[2,0] * this[3,2]
- this[0,0] * this[2,3] * this[3,2]
- this[0,2] * this[2,0] * this[3,3]
- this[0,3] * this[2,2] * this[3,0];
inverse[1,2] =
+ this[0,0] * this[1,3] * this[3,2]
+ this[0,2] * this[1,0] * this[3,3]
+ this[0,3] * this[1,2] * this[3,0]
- this[0,0] * this[1,2] * this[3,3]
- this[0,2] * this[1,3] * this[3,0]
- this[0,3] * this[1,0] * this[3,2];
inverse[1,3] =
+ this[0,0] * this[1,2] * this[2,3]
+ this[0,2] * this[1,3] * this[2,0]
+ this[0,3] * this[1,0] * this[2,2]
- this[0,0] * this[1,3] * this[2,2]
- this[0,2] * this[1,0] * this[2,3]
- this[0,3] * this[1,2] * this[2,0];
inverse[2,0] =
+ this[1,0] * this[2,1] * this[3,3]
+ this[1,1] * this[2,3] * this[3,0]
+ this[1,3] * this[2,0] * this[3,1]
- this[1,0] * this[2,3] * this[3,1]
- this[1,1] * this[2,0] * this[3,3]
- this[1,3] * this[2,1] * this[3,0];
inverse[2,1] =
+ this[0,0] * this[2,3] * this[3,1]
+ this[0,1] * this[2,0] * this[3,3]
+ this[0,3] * this[2,1] * this[3,0]
- this[0,0] * this[2,1] * this[3,3]
- this[0,1] * this[2,3] * this[3,0]
- this[0,3] * this[2,0] * this[3,1];
inverse[2,2] =
+ this[0,0] * this[1,1] * this[3,3]
+ this[0,1] * this[1,3] * this[3,0]
+ this[0,3] * this[1,0] * this[3,1]
- this[0,0] * this[1,3] * this[3,1]
- this[0,1] * this[1,0] * this[3,3]
- this[0,3] * this[1,1] * this[3,0];
inverse[2,3] =
+ this[0,0] * this[1,3] * this[2,1]
+ this[0,1] * this[1,0] * this[2,3]
+ this[0,3] * this[1,1] * this[2,0]
- this[0,0] * this[1,1] * this[2,3]
- this[0,1] * this[1,3] * this[2,0]
- this[0,3] * this[1,0] * this[2,1];
inverse[3,0] =
+ this[1,0] * this[2,2] * this[3,1]
+ this[1,1] * this[2,0] * this[3,2]
+ this[1,2] * this[2,1] * this[3,0]
- this[1,0] * this[2,1] * this[3,2]
- this[1,1] * this[2,2] * this[3,0]
- this[1,2] * this[2,0] * this[3,1];
inverse[3,1] =
+ this[0,0] * this[2,1] * this[3,2]
+ this[0,1] * this[2,2] * this[3,0]
+ this[0,2] * this[2,0] * this[3,1]
- this[0,0] * this[2,2] * this[3,1]
- this[0,1] * this[2,0] * this[3,2]
- this[0,2] * this[2,1] * this[3,0];
inverse[3,2] =
+ this[0,0] * this[1,2] * this[3,1]
+ this[0,1] * this[1,0] * this[3,2]
+ this[0,2] * this[1,1] * this[3,0]
- this[0,0] * this[1,1] * this[3,2]
- this[0,1] * this[1,2] * this[3,0]
- this[0,2] * this[1,0] * this[3,1];
inverse[3,3] =
+ this[0,0] * this[1,1] * this[2,2]
+ this[0,1] * this[1,2] * this[2,0]
+ this[0,2] * this[1,0] * this[2,1]
- this[0,0] * this[1,2] * this[2,1]
- this[0,1] * this[1,0] * this[2,2]
- this[0,2] * this[1,1] * this[2,0];
double factor = (1.0 / determinant);
for (int i = 0; i < inverse.totalRows; i++)
{
for (int j = 0; j < inverse.totalColumns; j++)
{
inverse[i, j] *= factor;
}
}
return (inverse);
}
public MatrixF64 Inverse()
{
if
(
(this.totalRows <= 0)
|| (this.totalColumns <= 0)
|| (null == this.entries)
)
{
return (MatrixF64.Zero( 1, 1 )); // Matrix is empty.
}
if (this.totalRows != this.totalColumns)
{
// Matrix is not square
return (MatrixF64.Zero( 1, 1 ));
}
// Simple cases
if (1 == this.totalRows) { return (this.Inverse1x1()); }
if (2 == this.totalRows) { return (this.Inverse2x2()); }
if (3 == this.totalRows) { return (this.Inverse3x3()); }
if (4 == this.totalRows) { return (this.Inverse4x4()); }
int[] rowPermutationOfOriginalMatrix = null;
bool rowExchangeParityIsOdd = false;
MatrixF64 originalMatrix = new MatrixF64( this );
MatrixF64 LUFactorsMatrix = null;
MatrixF64.FindLUFactorsOfARowPermutedVersionOfAnOriginalMatrix
(
originalMatrix,
ref LUFactorsMatrix,
ref rowPermutationOfOriginalMatrix,
ref rowExchangeParityIsOdd
);
if ((null == LUFactorsMatrix)
|| (null == rowPermutationOfOriginalMatrix))
{
return (MatrixF64.Zero( this.totalRows, this.totalColumns ));
}
VectorF64 columnVector = VectorF64.Zero( totalRows );
VectorF64 solutionVector = VectorF64.Zero( totalRows );
MatrixF64 inverse =
MatrixF64.Zero( this.totalRows, this.totalColumns );
for (int j = 0; j < inverse.totalColumns; j++)
{
for (int i = 0; i < inverse.totalRows; i++)
{
columnVector[i] = 0.0;
}
columnVector[j] = 1.0;
MatrixF64.LUBacksubstitution
(
LUFactorsMatrix,
rowPermutationOfOriginalMatrix,
columnVector,
ref solutionVector
);
if (null != solutionVector)
{
for (int i = 0; i < inverse.totalRows; i++)
{
inverse[i, j] = solutionVector[i];
}
}
}
return (inverse);
}
// . . .
public static void Test()
{
// . . .
// Examples of matrix inverses:
MatrixF64 m2x2 =
new MatrixF64
(
2, 2,
2.0, 5.0, // row 0
-4.0, -7.0 // row 1
);
MatrixF64 inversem2x2 = m2x2.Inverse();
inversem2x2.WriteLine();
// [ -1.1666667, -0.83333333 ]
// [ 0.66666667, 0.33333333 ]
MatrixF64 prodinvm2x2andm2x2 = inversem2x2 * m2x2;
prodinvm2x2andm2x2.WriteLine();
// [ 1, 8.8817842e-16 ]
// [ 0, 1 ]
// This product is very close to a
// 2x2 identity matrix, as it should be.
MatrixF64 prodm2x2andinvm2x2 = m2x2 * inversem2x2;
prodm2x2andinvm2x2.WriteLine();
// [ 1, 0 ]
// [ 0, 1 ]
// This product is the 2x2 identity matrix,
// as it should be.
MatrixF64 m3x3 =
new MatrixF64
(
3, 3,
2.0, 5.0, 7.0, // row 0
-4.0, -1.0, 6.0, // row 1
9.0, 8.0, 3.0 // row 2
);
MatrixF64 inversem3x3 = m3x3.Inverse();
inversem3x3.WriteLine();
// [ -0.76119403, 0.6119403, 0.55223881 ]
// [ 0.98507463, -0.85074627, -0.59701493 ]
// [ -0.34328358, 0.43283582, 0.26865672 ]
MatrixF64 prodinvm3x3andm3x3 = inversem3x3 * m3x3;
prodinvm3x3andm3x3.WriteLine();
// [ 1, 0, 4.4408921e-16 ]
// [ 0, 1, 0 ]
// [ 0, 0, 1 ]
// This product is very close to a 3x3 identity matrix,
// as it should be.
MatrixF64 prodm3x3andinvm3x3 = m3x3 * inversem3x3;
prodm3x3andinvm3x3.WriteLine();
// [ 1, -4.4408921e-16, 4.4408921e-16 ]
// [ 0, 1, -2.220446e-16 ]
// [ -4.4408921e-16, 0, 1 ]
// This product is very close to a 3x3 identity matrix,
// as it should be.
// . . .
}
}
3.19 Multiplicativo producto de una matriz y un vector
El producto multiplicativo de una matriz (M) y (a) un vector es un vector (b); (M)*(a) = (b).
Matrix (M) “transforma” vector (a) a (b) vector.
El código siguiente calcula el producto multiplicativo de una matriz y un vector.
public class MatrixF64
{
// . . .
public static VectorF64 operator *( MatrixF64 m, VectorF64 v )
{
if (null == m)
{
return (null); // Matrix is not specified.
}
if
(
(null == m.entries)
|| (m.totalRows <= 0)
|| (m.totalColumns <= 0)
)
{
return (null); // Matrix is empty.
}
if (null == v)
{
// Vector is not specified.
return (VectorF64.Zero( m.totalRows ));
}
if (v.Dimensions() != m.totalColumns)
{
// The number of columns of the matrix must be equal
// to the number of rows of the vector.
return (VectorF64.Zero( m.totalRows ));
}
// The result vector has the same number of rows as the matrix.
VectorF64 result = VectorF64.Zero( m.totalRows );
// Component (i) of the result vector is equal to the
// dot product of row (i) of the matrix with the vector.
for (int i = 0; i < m.totalRows; i++)
{
double dotProduct = 0.0;
for (int j = 0; j < m.totalColumns; j++)
{
dotProduct += (m[i, j] * v[j]);
}
result[i] = dotProduct;
}
return (result);
}
// . . .
public static void Test()
{
// . . .
// Examples of a product of a matrix and a vector:
MatrixF64 t3x3 =
new MatrixF64
(
3, 3,
2.0, 5.0, 7.0, // row 0
-4.0, -1.0, 6.0, // row 1
9.0, 8.0, 3.0 // row 2
);
VectorF64 p3x1 = new VectorF64( 1.0, 2.0, 3.0 );
VectorF64 q3x1 = t3x3 * p3x1;
q3x1.WriteLine();
// ( 33, 12, 34 )
// The multiplicative product of the inverse
// of the matrix and the result vector computed
// above should reproduce the original vector.
MatrixF64 inverset3x3 = t3x3.Inverse();
VectorF64 prodinvt3x3andq3x1 = inverset3x3 * q3x1;
prodinvt3x3andq3x1.WriteLine();
// ( 1, 2, 3 )
// This vector is equal to the original vector,
// as expected.
// . . .
}
}
3.20 El uso de matrices homogéneas para traducir homogénea vectores
“traducción:” Una operación T(s), especificado por un vector (s), que convierte cualquier vector (r) correspondiente a un vector (r + s); T(s) * (r) = (r + s).
T(-s) es la inversa de T(s); Inverse(T(s)) = T(-s); T(-s) * (r) = (r - s).
Una operación de traducción pueden ser representados por una matriz, de tal manera que la multiplicación de la traducción matriz por un vector tiene el efecto de que la traducción de vector.
Las contribuciones de una operación de traducción para los componentes del resultado de vectores no dependen, directa o indirectamente, a los componentes del vector original.
De este modo, una matriz de traducción debe ser estructurado de manera que le permite contribuir con una traducción al resultado sin ser afectados, directa o indirectamente, por los componentes del vector original.
Un “vector homogéneo” es un vector con un último elemento igual a 1.
Una “matriz homogénea” es una matriz con una última fila con todas las entradas 0 igual a excepción de una última entrada de 1.
El producto multiplicativo de una matriz homogénea y otra matriz homogénea es una matriz homogénea.
El multiplicador de un producto homogéneo y de una matriz homogénea vector es un vector homogéneo.
Un atributo interesante de este producto es que las entradas en la última columna de la matriz homogénea, con la excepción de la última entrada, se añade directamente a los componentes del vector resultado, sin ser afectados, directa o indirectamente, por los componentes de el vector original.
De este modo, una operación de traducción pueden ser representados por una matriz homogénea, y esa matriz se puede utilizar para traducir homogénea vectores.
Si queremos representar las posiciones y direcciones en (d)-dimensional del espacio, y queremos usar matrices para traducir las posiciones, luego homogéneo matrices y vectores homogéneo se pueden utilizar.
Sin embargo, los vectores homogéneos debe tener (d+1) componentes, y las matrices homogéneas deben tener ((d+1)*(d+1)) entradas.
El componente final de un vector homogéneo y la última fila de una matriz homogénea, no corresponde a ninguna de las (d) coordenadas que queremos representar en (d) dimensiones espaciales, sino que sólo sirven para que sea posible obtener matrices para hacer algunas operaciones independientemente de la (d) coordinar valores.
El uso de matrices homogéneas y homogénea vectores hace que sea más fácil de compuestos transformaciones y aplicar transformaciones a los vectores, sino otras operaciones cada vez más incómodo.
En representación de posiciones o direcciones en (d) dimensiones espaciales utilizando vectores homogéneos significa que los vectores se han (d+1) componentes, con el último elemento igual a 1.
Cálculos el producto escalar de dos vectores homogéneos utilizando un punto-producto función de combinar implícitamente el producto escalar de la primera (d) componentes con un valor irrelevante 1 de final debido a los componentes de 1.
Por lo tanto, para obtener el producto escalar de los (d) componentes de la homogénea vector (con (d+1) componentes), el producto escalar de vectores homogéneos debe llevarse a cabo mediante un “vector homogéneo producto escalar” función que simplemente se salta el componente final de los vectores.
Otra posibilidad es que el producto escalar puede ser calculada con el producto escalar ordinario de funcionar con el entendimiento de que la suma de 1 debe restarse del resultado antes de interpretar el resultado en lo que se refiere a la (d) las coordenadas pertinentes.
El código siguiente muestra cómo una matriz se pueden formar que puede utilizarse para traducir un vector.
public class MatrixF64
{
// . . .
public static MatrixF64 FormHomogeneousTranslationMatrix( VectorF64 v )
{
// This method assumes the specified vector is a homogeneous
// vector, where the final component is not relevant to the
// translation. Thus, we form a square homogeneous matrix
// with a total number of rows equal to the number of
// components of the specified vector. We form an identity
// matrix of the required size, and copy all but the final
// component of the specified vector to the final column of
// the matrix.
if (null == v)
{
return (MatrixF64.Zero( 1, 1 )); // Vector is not specified.
}
if (v.Dimensions() <= 0)
{
return (MatrixF64.Zero( 1, 1 )); // Vector is empty.
}
MatrixF64 result =
MatrixF64.Identity( v.Dimensions() );
// Add all but the final component of the specified vector
// to the final column of the result matrix. The final
// entry of the final column of the matrix will remain one (1).
for (int i = 0; i < (result.totalRows - 1); i++)
{
result[i, (result.totalColumns - 1)] = v[i];
}
return (result);
}
// . . .
public static void Test()
{
// . . .
// Example of using a homogeneous matrix to
// translate a homogeneous vector:
// Forming a homogeneous matrix that represents
// a translation:
VectorF64 homogeneousTranslationVector4x1 =
new VectorF64( 10.0, 20.0, 30.0, 1.0 );
MatrixF64 translation4x4 =
MatrixF64.FormHomogeneousTranslationMatrix
( homogeneousTranslationVector4x1 );
translation4x4.WriteLine();
// [ 1, 0, 0, 10 ]
// [ 0, 1, 0, 20 ]
// [ 0, 0, 1, 30 ]
// [ 0, 0, 0, 1 ]
// Using the homogeneous matrix to translate
// a homogeneous vector:
VectorF64 homogeneousVector =
new VectorF64( 5.0, 6.0, 7.0, 1.0 );
VectorF64 translatedVector =
translation4x4 * homogeneousVector;
translatedVector.WriteLine();
// ( 15, 26, 37, 1 )
// The relevant vector components (5,6,7)
// were translated by (10,20,30) by the
// homogeneous translation matrix.
// . . .
}
}
3.21 (d) dimensiones espaciales: rotaciones
“(d) rotación en el espacio-dimensional:” Una operación (R) especificado por un período de dos dimensiones (S) plano de rotación, un ángulo de rotación (t), y un centro de rotación punto (c), que convierte cualquier punto (p) a un punto correspondiente (q), tal que:
(1) Length( q - c ) = Length( p - c ), la distancia entre el nuevo punto y el punto centro de rotación es la misma que la distancia entre el punto inicial y el punto central de rotación.
De este modo, el nuevo punto se encuentra sobre la superficie de un particular (d) dimensiones esfera;
(2) Que se (P) un plano que es paralelo al plano de rotación y (S) contiene punto (p).
Punto (q) se (P) en avión.
De este modo, el nuevo punto se encuentra en una de dos dimensiones en el plano (d) dimensiones espaciales;
(3) Dot( (q - c), (p - c) ) = Length( q - c ) * Length( p - c ) * Cos( t ); el producto escalar del vector de la rotación punto central para el punto original y el vector de la rotación punto central para el nuevo punto es igual al producto de las longitudes de los dos vectores y el coseno de la rotación ángulo.
De este modo, el nuevo punto se encuentra sobre la superficie de un particular (d) cono-dimensional;
Constraint (1) limita el nuevo punto de la superficie de un particular (d) dimensiones esfera.
Condición (2) aún más el nuevo punto a una de dos dimensiones en el plano (d) dimensiones espaciales.
Por lo tanto, el nuevo punto se debe a una de dos dimensiones en círculo (d) dimensiones espaciales.
Condición (3) aún más el nuevo punto a una determinada (d) cono-dimensional.
Por lo tanto, el nuevo punto debe estar en un punto en particular (por Cos( t ) = 1 o Cos( t ) = (-1)) o en uno de los dos puntos posibles (para (-1) < Cos( t ) < 1).
Una última limitación está obligado a elegir uno de los dos puntos posibles que se derivan de las limitaciones mencionadas anteriormente.
Que T(s) representa una operación de “traducción,” donde (s) es un vector arbitrario, de tal manera que T(s) * (r) = (r + s), y Inverse(T(s)) * (r) = (r - s).
Una rotación R(S,t,c), por avión S, t ángulo, centro y punto c, puede expresarse en términos de operaciones de traducción y con una rotación de un punto central en el origen: R(S,t,c) = T(c) * R(S,t,0) * T(-c) = T(c) * R(S,t,0) * Inverse(T(c)).
Que R(S,t) representa una operación de rotación con un centro de rotación en el punto de origen (d) dimensiones espaciales.
Por lo tanto, cualquier rotación arbitraria se puede expresar como: R(S,t,c) = T(c) * R(S,t) * T(-c) = T(c) * R(S,t) * Inverse(T(c)).
Vamos a considerar dos dimensiones (S) rotación de aviones que incluyen el origen de (d) dimensiones espaciales y son paralelas a dos ejes de coordenadas.
Que (a) y (b) ser enteros que representan ejes de coordenadas en (d)-dimensional del espacio, de tal manera que 0 <= a <= (d-1), y 0 <= b <= (d-1), y (a != b).
Se define R(a,b,t) a ser una rotación sobre el origen de (d)-dimensional del espacio, con una rotación plano paralelo a ejes de coordenadas indicadas por coordinar los índices (a) y (b), y con un ángulo (t), de tal manera que R(a,b,(pi/2)) sería convertir un valor positivo para coordinar la coordenada eje indicado por (a) a un valor positivo para coordinar coordinar el eje indicado por (b); ejemplo: R(a,b,(pi/2)) * e(a) = e(b), donde e(k) indica una unidad de base de vectores correspondientes a fin de coordinar eje (k).
Esta definición proporciona la última disambiguation necesarios para determinar el nuevo punto producido por una rotación.
Una rotación que se mueve a lo largo del eje hacia a eje b de acuerdo con un ángulo t es equivalente a una rotación que se mueve a lo largo de eje b hacia a eje de un ángulo (-t), es decir, R(a,b,t) = R(b,a,(-t)).
Si (d >= 2), entonces (d) dimensiones espaciales ha ((d*(d-1))/2) pares de ejes de coordenadas.
Así, por (d >= 2), (d) dimensiones espaciales ha ((d*(d-1))/2) distinto eje de rotación de aviones alineados.
(Ejemplos: d=2: () x-y avión; d=3: (x-y avión, y-z plano, plano z-x); d=4: (x-y avión, y-z avión, z-x avión, x-w avión, z-w plano, plano w-y).)
El código siguiente muestra cómo una matriz puede ser utilizado para representar la rotación R(a,b,t).
Un ejemplo se muestra que utiliza una matriz de rotación de manera homogénea a girar una homogénea vector a un nuevo vector homogéneo.
public class MatrixF64
{
// . . .
public static MatrixF64 Rab( int rows, int a, int b, double angle )
{
if (rows <= 0)
{
// Invalid row count specified.
return (MatrixF64.Zero( 1, 1 ));
}
MatrixF64 result = MatrixF64.Identity( rows );
// If (a==b) (rotate coordinate (a) to (a)), then
// there is no rotation to do.
if ((a != b) && (a >= 0) && (a < rows) && (b >= 0) && (b < rows))
{
double sine = Math.Sin( angle );
double cosine = Math.Cos( angle );
// Override m[a,a] = cosine; m[a,b] = -sine;
// m[b,a] = sine; m[b,b] = cosine;
result[a, a] = cosine;
result[b, b] = cosine;
result[a, b] = (-(sine));
result[b, a] = sine;
}
return (result);
}
public static MatrixF64 Rabk( int rows, int a, int b, int k )
{
double angle = (Math.PI / 2.0) * (double)k;
return (MatrixF64.Rab( rows, a, b, angle ));
}
// . . .
public static void Test()
{
// . . .
// Example of using a homogeneous matrix
// to rotate a homogeneous vector:
// Forming a homogeneous rotation matrix that
// rotates a positive component along coordinate
// axis 0 to a positive component along
// cooridnate axis 1, in 3-dimensional space.
// (The homogeneous matrix must be (3+1)=4
// rows by (3+1)=4 columns.)
MatrixF64 rotation4x4 =
MatrixF64.Rab( 4, 0, 1, (0.5 * Math.PI) );
rotation4x4.WriteLine(3);
// [ 6.12e-17, -1, 0, 0 ]
// [ 1, 6.12e-17, 0, 0 ]
// [ 0, 0, 1, 0 ]
// [ 0, 0, 0, 1 ]
// (The entries '6.12e-17' are negligible.)
// The upper-left 3x3 part of the matrix represents
// a counterclockwise rotation of a quarter-turn
// such that a positive component along axis 0
// would be rotated to a positive component along
// axis 1, in 3-dimensional space.
// Using the homogeneous rotation matrix to
// rotate a homogeneous vector:
VectorF64 homogeneousVector2 =
new VectorF64( 5.0, 6.0, 7.0, 1.0 );
VectorF64 rotatedVector =
rotation4x4 * homogeneousVector2;
rotatedVector.WriteLine();
// ( -6, 5, 7, 1 )
// The relevant vector components (5,6,7) were
// rotated to (-6,5,7), which is consistent with
// a counterclockwise quarter-turn rotation
// with a rotation plane parallel to axes 0 and 1.
// . . .
}
}
3.22 Propuesta de definición del “contrario”
En dos dimensiones espaciales, los cuartos de turno de rotación que puede convertir el vector e(0) a e(1) se considera contrario.
La rotación es: R(0,1,(pi/2)) * e(0) = e(1).
Por lo tanto, la rotación es R(0,1,t) contrario.
En el espacio tridimensional, los cuartos de turno de rotación que puede convertir el vector e(0) a e(1) se considera contrario.
El cuarto de rotación que a su vez puede convertir el vector e(1) a e(2) también se considera contrario.
El cuarto de rotación que a su vez puede convertir el vector e(2) a e(0) también se considera contrario.
Por lo tanto, las rotaciones son { R(0,1,t), R(1,2,t), R(2,0,t) } contrario.
El operador booleano fórmula se ajusta a continuación la definición del contrario para bidimensional y tridimensional del espacio, y extrapola el modelo para el caso general de la rotación R(a,b,t):
bool counterClockwise =
(
( (a < b) && (((a + b) % 2) == 1) )
|| ( (a > b) && (((a + b) % 2) == 0) )
);
Esta norma se define por debajo de las rotaciones como “en sentido antihorario:”
R(0,1,t),
R(2,0,t), R(1,2,t),
R(0,3,t), R(3,1,t), R(2,3,t),
R(4,0,t), R(1,4,t), R(4,2,t), R(3,4,t),
R(0,5,t), R(5,1,t), R(2,5,t), R(5,3,t), R(4,5,t),
...
3.23 (d) dimensiones espaciales: rotaciones en sentido horario
Utilizando la definición propuesta de “contrario” en la sección anterior, las rotaciones están por debajo de “las agujas del reloj:”
R(1,0,t),
R(0,2,t), R(2,1,t),
R(3,0,t), R(1,3,t), R(3,2,t),
R(0,4,t), R(4,1,t), R(2,4,t), R(4,3,t),
R(5,0,t), R(1,5,t), R(5,2,t), R(3,5,t), R(5,4,t),
...
Por ejemplo, en cuatro dimensiones espaciales, los seis matrices de rotación en sentido horario a continuación representan las rotaciones:
R(1,0,t) = Cos(t), Sin(t), 0, 0,
-Sin(t), Cos(t), 0, 0,
0, 0, 1, 0,
0, 0, 0, 1
R(0,2,t) = Cos(t), 0, -Sin(t), 0,
0, 1, 0, 0,
Sin(t), 0, Cos(t), 0,
0, 0, 0, 1
R(2,1,t) = 1, 0, 0, 0,
0, Cos(t), Sin(t), 0,
0, -Sin(t), Cos(t), 0,
0, 0, 0, 1
R(3,0,t) = Cos(t), 0, 0, Sin(t),
0, 1, 0, 0,
0, 0, 1, 0,
-Sin(t), 0, 0, Cos(t)
R(1,3,t) = 1, 0, 0, 0,
0, Cos(t), 0, -Sin(t),
0, 0, 1, 0,
0, Sin(t), 0, Cos(t)
R(3,2,t) = 1, 0, 0, 0,
0, 1, 0, 0,
0, 0, Cos(t), Sin(t),
0, 0, -Sin(t), Cos(t)
3.24 (d) dimensiones del espacio: número total de ejes alineados orientaciones
El número de distintos ejes alineados orientaciones en (d) dimensiones espaciales es:
distinctOrientations = Factorial( d ) * IntPow( 2, (d-1) );
Un argumento de apoyo a esta fórmula figura en los párrafos siguientes.
Evidencia experimental verificar la fórmula también se presentan a continuación.
Que { g(0), g(1), ..., g(d-1) } ser un conjunto de (d) mutuamente perpendiculares vectores unitarios adjunta a un objeto en (d) dimensiones espaciales.
Estos vectores unitarios indicar la orientación del objeto en relación con los ejes de coordenadas de la (d) dimensiones espaciales.
Deje que la “orientación predeterminada” del objeto ser tal que g(i) = e(i), para todos los 0 <= i <= (d-1); cada vector de orientación del objeto está asociado a una unidad distinta de vectores base de un sistema de coordenadas en la (d) dimensiones espaciales.
El primer vector de orientación del objeto puede ser paralelo a cualquiera de los ejes de (d) la (d) dimensiones espaciales, y puede ser en cualquiera de dos direcciones en relación con ese eje.
Por lo tanto, el primer vector de orientación ha (2 * d) opciones.
El segundo vector de orientación del objeto puede ser paralela a cualquiera de los restantes (d-1) ejes, con dos posibles direcciones en relación con ese eje.
Por lo tanto, el segundo vector ha orientación (2 * (d-1)) opciones.