Álgebra Linear
Colin Fahey
1. Software
LinearAlgebra.zip
Álgebra Linear código fonte (C#)
19910 bytes
MD5: 11d8c8035cac30ba543e5e0b72ee9767
2. Introdução
Este artigo descreve vetores e matrizes, em (d)-dimensional espaço.
3. (d)-dimensional espaço: atributos
3.1 Array
“array:” Uma coleção de variáveis de modo a que cada variável tem um nome único, e de tal modo que os nomes podem ser atribuídos um fim.
Valores inteiros podem ser usadas como nomes para as variáveis em um array.
Por exemplo, se um array contém (d) variáveis e, em seguida, os valores inteiros { 0, 1, 2, ..., (d-1) } os nomes podem ser atribuídos às variáveis do array.
3.2 (d)-dimensional vector
“(d)-dimensional vector:” Um conjunto de variáveis (d).
“vetor componente:” Uma variável dentro de um vector.
3.3 (d)-dimensional vector espacial
“one-dimensional espaço:” O conjunto completo de valores que podem ser armazenados em uma variável.
“(d)-dimensional espaço:” O conjunto completo de combinações de valores que podem ser armazenadas em um array de (d) variáveis.
A definição formal de um “espaço vetorial:”
Deixe (T) ser um tipo básico (exemplos: número real, inteiro, número complexo, número racional, etc).
Qualquer uma variável do tipo básico é nomeado um “escalar.”
Um “(T)-tipo (d)-dimensional vector espaço” é o conjunto (S) de (d)-dimensional vetores ter duas operações, além (+) vetorial, escalar e multiplicação (*), satisfazendo as condições abaixo:
(1) Se (v) e (w) são quaisquer dois vectores em (S) e, em seguida, (v + w) é também um vector de (S);
(2) Se (u), (v), e (w) são três vectores, em qualquer (S) e, em seguida, (u + v) + w = u + (v + w);
[aditivo Comutatividade]
(3) Se (v) e (w) são quaisquer dois vectores em (S) e, em seguida, (v + w) = (w + v);
[aditivo associatividade]
(4) Existe um “vetor zero,” (0), em (S), de tal forma que para qualquer vector (v) em (S), (v + (0)) = v;
[aditivo identidade]
(5) Se (c) é qualquer tipo de escalar (T), e (v) é qualquer vetor em (S) e, em seguida, o produto (c * v) é um vetor em (S);
(6) Se (a), (b), e (c) são escalares de qualquer tipo (T), e (v) e (w) são vetores em qualquer (S) e, em seguida, (a + b) * v = a*v + b*v, e c*(v + w) = c*v + c*w;
[multiplicativo distributivity]
(7) Se (a) e (b) são escalares de qualquer tipo (T), e (v) é qualquer vetor em (S) e, em seguida, (a*b)*v = a*(b*v);
(8) Se “1” é um tipo de escalar (T) tal que (1*1)=1, e (v) é qualquer vetor em (S) e, em seguida, (1*v) = v;
(9) Para cada vector (v) em (S), o vetor (-1)*v = -v satisfaz v + (-v) = (0);
[aditivo inversa]
3.4 (d)-dimensional vector código
O código abaixo mostra como um vetor (d)-dimensional, com 64 bits de ponto flutuante de componentes, pode ser implementada.
Um array é suficiente para representar um vetor.
O código abaixo tem um array contidos em uma classe, apenas por conveniência.
O código não se destina a ser eficientes.
Uma estrutura (exemplo: “struct”, um tipo de valor), representando vetores de um número fixo de dimensões (exemplos: 3 ou 4) é provável que seja muito mais eficiente do que a geral (d)-dimensional vector classe mostrado aqui.
Apesar de o código abaixo define um vetor com componentes de ponto flutuante, este documento também faz uso de vetores com componentes inteiro.
O código abaixo pode ser facilmente modificado para implementar vetores com componentes inteiro.
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) zero-dimensional vector:” Um vetor com todos os componentes (d) igual a zero.
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)-dimensional espaço: pontos
“(d)-dimensional espaço ponto:” um arranjo de (d) variáveis com valores específicos “(valores coordenar).”
“(d)-dimensional espaço origem:” um arranjo de (d) todas as variáveis com valores iguais a zero.
3.6 (d)-dimensional vectores: não-parente e parente
Um (d)-dimensional vector que “não” é “parente” é um vetor (d)-dimensional que representa directamente a um estatuto ou configuração.
Um (d)-dimensional vector que é “relativa” (d)-dimensional é um vetor que representa a evolução de um conjunto de componentes.
Um parente vetor pode representar a diferença entre os dois vetores não-parente.
Dado um parente vetor, determinando um estado não-parente ou configuração usando esse parente vetor exige que combinam com um parente vetor não-parente vetor.
Não-parente vetores e relativos são vetores dois vectores.
Se um determinado vetor é não-familiar ou parente deve ser especificado quando o vector está definido.
Se um (d)-dimensional vector é interpretado como sendo não-parente, então o (d)-dimensional vector pode representar um ponto no espaço (d)-dimensional.
3.7 (d)-dimensional vector adição, subtração, e escamação
Vector adição, subtração, e escamação:
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 vetores
A definição formal de uma “base” de um espaço vetorial:
Deixe (T) ser um tipo básico (exemplos: número real, inteiro, número complexo, número racional, etc).
Qualquer uma variável do tipo básico é nomeado um “escalar.”
Deixe (V) ser um “(T)-tipo (d)-dimensional vector espaço.”
Se o não-zero vetores { u1, u2, ..., ud } em (V) são de tal ordem que cada vector (v) em (V) pode ser escrito como uma “combinação linear” desses vetores, v = c1*u1 + c2*u2 + ... + cd*ud, onde estão { c1, c2, ..., cd } escalares do tipo (T) e, em seguida, (V) é “calibrado” por vetores { u1, u2, ..., ud }.
Qualquer conjunto de vetores não-zero { u1, u2, ..., ud } que “abrangem” espaço (V) vetor é chamado de uma “base” de (V).
Uma simples “função” de um espaço vetorial (d)-dimensional é o conjunto de (d) distintas (d)-dimensional vetores, cada um com um elemento igual a um e de todos os outros componentes igual a zero.
Essas são vetores “orthonormal” base, o que significa que eles são mutuamente “perpendiculares-(ortogonais),” e que cada vector tem unidade de comprimento.
Cada um desses vetores é um vetor unitário paralelo a um dos eixos (d) coordenar.
Expressando um vector arbitrário como uma “combinação linear” destes vetores base é direta; cada componente do vetor arbitrário é multiplicado pela correspondente base vetorial, e esses produtos são acrescentados juntos para formar o vetor arbitrário.
O código abaixo mostra como um vetor pode ser expresso como uma “combinação linear” dos vetores base.
O código abaixo arbitrariamente orthonormal define um conjunto de vetores 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 )
// . . .
}
}
Qualquer vetor em (d)-dimensional espaço pode ser expressa como uma soma de produtos a base de números e vetores:
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)-dimensional espaço: distância entre os pontos
Deixe (P) ser um vector (d)-dimensional que representa um ponto no espaço (d)-dimensional.
Deixe (Q) ser um vector (d)-dimensional que representa um ponto no espaço (d)-dimensional.
Deixe (R) ser um (d)-dimensional vector que representa a mudança na (d) coordenadas para chegar ao ponto de (P) ponto (Q); R = (Q - P).
Em um espaço de dimensão, P = ( p0 ), Q = ( q0 ), e R = (Q - P) = ( q0 - p0 ).
A distância entre os dois pontos é a seguinte: Abs( q0 - p0 ).
Em duas dimensões espaciais, P = ( p0, p1 ), Q = ( q0, q1 ), e R = (Q - P) = ( q0-p0, q1-p1 ).
Interpretando as duas perpendicular deslocamentos como a perpendicular lados de um triângulo de direita, a distância entre os pontos corresponde ao comprimento da hipotenusa desse triângulo.
A fórmula Pythagorean, (a*a) + (b*b) = (c*c), onde (a) e (b) são os comprimentos dos lados de uma perpendicular à direita-triângulo, e (c) é o comprimento da hipotenusa (deforme lado), podem ser utilizados para determinar a distância entre os dois pontos: Sqrt( Sq(q0-p0) + Sq(q1-p1) ).
No espaço tridimensional, P = ( p0, p1, p2 ), Q = ( q0, q1, q2 ), e R = (Q - P) = ( q0-p0, q1-p1, q2-p2 ).
Interpretando o deslocamento ( q0-p0, q1-p1, 0 ) como a perpendicular lados de um triângulo de direita, e usando a fórmula Pythagorean, dá a distância entre o ponto a ponto (P) e ( q0, q1, p2 ): d01 = Sqrt( Sq(q0-p0) + Sq(q1-p1) ).
O deslocamento ( q0-p0, q1-p1, 0 ) é perpendicular ao deslocamento ( 0, 0, q2-p2 ), e outro de direita triângulo pode ser formada, e os Pythagorean fórmula pode ser usado novamente.
Assim, a distância do ponto (P) ao ponto (Q) é 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) ).
O método de extensão a distância fórmula de duas dimensões do espaço a espaço tridimensional pode ser aplicado repetidamente para, eventualmente, determinar a distância fórmula para (d)-dimensional espaço: Sqrt( (Sq(q0-p0) + Sq(q1-p1)) + Sq(q2-p2) + ... + Sq(qd-pd) ).
O código abaixo define uma função que calcula o chamado “Comprimento” comprimento de um vetor (d)-dimensional.
Quando um vetor representa o deslocamento entre dois pontos em (d)-dimensional espaço, o comprimento do vetor representa a distância entre esses dois pontos.
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)-dimensional vetores: produto escalar
O “produto escalar” converte-dimensional (d) dois vetores de um número.
O código abaixo calcula um produto escalar de dois 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 qualquer vector, (A), Length(A) = Sqrt(Dot(A,A)).
3.11 (d)-dimensional vectores: definição do “paralelo”
Vetores (A) e (B) são “paralelos” se todas as seguintes são verdadeiras:
(1) A.Length() > 0;
(2) B.Length() > 0;
(3) Abs(Dot(A,B)) = A.Length()*B.Length().
O código abaixo determina se um par de vetores são paralelos (possivelmente anti-alinhados).
Números de ponto flutuante pode acumular erros no fracionário parte devido à precisão limitada, e, portanto, deve incluir código computador diferente de zero quando se comparam as tolerâncias de ponto flutuante de números.
O código inclui exemplo tolerância valores, mas o exemplo tolerância valores poderá não ser apropriado para algumas tarefas.
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)-dimensional vectores: definição de “perpendicular”
Vetores (A) e (B) são “perpendiculares” se todas as seguintes são verdadeiras:
(1) A.Length() > 0;
(2) B.Length() > 0;
(3) Abs(Dot(A,B)) = 0.
O código abaixo determina se um par de vetores são perpendiculares. Números de ponto flutuante pode acumular erros no fracionário parte devido à precisão limitada, e, portanto, deve incluir código computador diferente de zero quando se comparam as tolerâncias de ponto flutuante de números.
O código inclui exemplo tolerância valores, mas o exemplo tolerância valores poderá não ser apropriado para algumas tarefas.
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 Matrizes
“matriz:” Uma coleção de variáveis de modo a que cada variável tem uma combinação única de uma “linha” e uma “coluna” nome nome.
“entrada:” uma variável dentro de uma matriz.
Valores inteiros podem ser usadas como nomes “fila” e “coluna” nomes para as variáveis em uma matriz.
Por exemplo, se uma matriz tem (totalRows) linhas e colunas (totalColumns), então os valores inteiros { 0, 1, ..., (totalRows-1) } pode ser o nome atribuído às fileiras, e os valores inteiros { 0, 1, ..., (totalColumns-1) } podem ser os nomes atribuídos às colunas.
Assim, uma variáveis em uma matriz pode ser especificado pela especificação de um par de inteiros, ( row, column ), indicando a combinação dos índices linha e coluna correspondente à variável específica.
O tamanho de uma matriz é especificada como “(totalRows) * (totalColumns)” (ou “(totalRows) por (totalColumns)).”
Este fim de dimensões é a mesma que a ordem das dimensões usadas para especificar as entradas na matriz “(( row, column )).”
[Esta convenção é um pouco infeliz, porque para muitos usos bidimensional (exemplos: imagens, gráficos, etc) é o habitual convenção para especificar as dimensões “width * height” e coordena as “( horizontal, vertical )” (ou “( x, y )).”
Este é o oposto da ordem de dimensões e coordenadas usado para descrever matrizes e os seus cadastros. ]
Uma matriz com (totalRows) igual a (totalColumns) é chamado “quadrado;” de outra forma, a matriz é chamado “retangular.”
A matriz pode ser considerada como contendo um conjunto de “vetores linha,” onde as variáveis em cada linha são interpretados como pertencentes a um vetor.
Uma matriz também pode ser considerado como contendo um conjunto de “vetores coluna,” onde as variáveis em cada coluna são interpretadas como pertencentes a um vetor.
Matrizes podem representar uma grande variedade de relações matemáticas.
O significado de uma matriz, e as operações que possam ser apropriadas para a transformação entradas de uma matriz, depende do contexto.
No entanto, existem regras básicas de aritmética matriz que sejam relevantes para muitos contextos, e estas regras básicas será definido em uma seção posterior.
Um array e um valor (totalColumns) são suficientes para representar uma matriz.
A matriz pode ter (totalRows * totalColumns) variáveis, e da entrada em ( row, column ) pode corresponder à matriz a variável Índice ((totalColumns * row) + column).
O código abaixo define uma matriz, com 64 bits de ponto flutuante de entradas.
Um array e um valor (totalColumns) são suficientes para representar uma matriz.
O código abaixo tem um array e um (totalColumns) valor contido em uma classe, apenas por conveniência.
O código abaixo não se destina a ser eficientes.
Uma estrutura (exemplo: “struct”, um tipo de valor), representando matrizes de especial dimensões fixas (exemplos: 2*2, 3*3, ou 4*4) é provável que seja muito mais eficiente do que a classe geral (totalRows * totalColumns) mostrado aqui.
Apesar de o código abaixo define uma matriz com entradas de ponto flutuante, este artigo também faz uso de matrizes com entradas inteiras.
O código abaixo pode ser facilmente modificado para implementar matrizes com entradas inteiras.
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
// . . .
}
}
“zero matriz:” uma matriz com todas as entradas iguais a zero.
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 Matrix adição, subtração, multiplicação e
Matrix adição, subtração, e multiplicação.
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 identidade:” Uma matriz quadrada (total de linhas equivale a total colunas) com entradas iguais a (1) sobre a diagonal (entradas com uma fila índice igual ao seu índice coluna), e com todas as outras entradas igual a (0).
Se uma matriz quadrada, (M), é multiplicada por uma “matriz identidade,” (I), do mesmo número de linhas (e colunas), o produto é igual a (M).
(I) multiplicando por (M) também produz um produto igual a (M).
Assim, a “matriz identidade” é análogo ao número “1” para a multiplicação dos números (escalares).
O código abaixo cria uma matriz identidade com um determinado número de linhas.
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; back-substitution
“Matrix LU factoring:” Um processo que converte uma matriz quadrada, (M), a nova praça duas matrizes, (L) e (U), tendo o mesmo tamanho que (M), de tal forma que (L) * (U) = (M), e tal que (U) matriz é “triangular superior” (todas as entradas estão abaixo da diagonal zero), e é “mais baixo” (L) matriz “triangular” (todas as entradas acima da diagonal são zero).
Matrix LU factoring pode ser usado como parte de um processo maior de resolver um sistema de equações, ou para encontrar a inversa de uma matriz, ou para encontrar o determinante de uma matriz.
Matrix LU factoring é uma alternativa para a eliminação “Gauss-Jordan procedimento.
Gauss-Jordan eliminação requer um sistema de equações (A)*(x)=(b), enquanto LU factoring apenas exige uma matriz (A).
Além disso, após a determinação do LU factoring de uma matriz (A), é muito fácil de determinar (x) dada qualquer (b).
O procedimento para resolver para o vetor (x) em (A)*(x)=(b), dada a (b) e LU factoring (L)*(U)=(A), envolve um processo que venha a ser chamado “de volta-substituição” aqui.
O código abaixo inclui um procedimento para a computação LU factores de qualquer quadrado, não-singular matriz.
O código abaixo inclui um procedimento para fazer back-LU substituição.
Atenção: O computador calcula o código abaixo L e U factores de uma linha de permuta versão de uma determinada matriz.
Além disso, a matriz L e U resultados são combinadas em uma única saída a matriz.
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:” Uma função que converte qualquer quadrado (n * n) matriz (A) para um número de telefone, det(A), de tal forma que:
(1) Se uma matriz (B) resultados de trocar duas linhas ou duas colunas, de uma matriz (A) e, em seguida, det(B) = (-(det(A)));
(2) Se uma matriz (B) multiplicando os resultados de qualquer fila, ou de qualquer coluna, de uma matriz (A), por uma série (c) e, em seguida, det(B) = c * det(A);
(3) Se os resultados a partir de uma matriz (B) adicionando um múltiplo de uma linha para outra linha, ou a partir de adicionar um múltiplo de uma coluna para outra coluna e, em seguida, det(B) = det(A).
(4) Se (A) é uma matriz triangular superior (A[i,j]=0 para todos os (i>j)) ou de uma matriz triangular inferior (A[i,j]=0 para todos os (i<j)) e, em seguida, det(A) = (A[0,0] * A[1,1] * ... * A[n-1,n-1]);
(Por exemplo, o determinante de uma matriz identidade é uma; det(I) = 1.)
Regras (1), (2), e (3), podem ser utilizadas, em um processo denominado “Gaussian eliminação, para converter qualquer matriz quadrada de uma matriz triangular.
Regra (4) pode ser usado para calcular o determinante de uma matriz triangular.
A definição formal de um “determinante:”
Deixe M[n](K) denotar o conjunto de todas as matrizes (n * n) ao longo do campo (K).
O “determinante” é a única função (F) tal que F:M[n](K) --> K com dois atributos:
(1) F é multilinear alternada no que diz respeito às colunas (ou linhas), e,
(2) F( I ) = 1.
Um “multilinear” função (D) de um (n * n) matriz (A) pode ser escrito 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]] ) ), onde a soma é tomada sobre todas as (n^n) combinações de 0 <= k[i] <= (n-1), e onde e[j] representa (j) fila da matriz identidade.
Essa função pode ser caracterizada pelos valores de D( e[k[0]], e[k[1]], ..., e[k[n-1]] ).
Uma função é uma “multilinear alternada” multilinear função, (D), de tal forma que D( ..., e[i], ..., e[i], ... ) = 0, e D( ..., e[i], e[j], ... ) = (-( D( ..., e[j], e[i], ... ) )).
Assim, trocando colunas (ou linhas) muda o sinal da função.
Além disso, o termo D( e[k[0]], ..., e[k[n-1]] ) é zero se a combinação de k[i] valores não implicar todos os valores exclusivos; se o conjunto de valores não é uma permutação dos números {0,1,...,(n-1)}.
Solte D( I ) = 1 (D( e[0], e[1], ..., e[n-1] ) = 1), em combinação com a alternância multilinear atributos descritos acima, significa que todos os valores dos multilinear alternada D( e[k[0]], e[k[1]], ..., e[k[n-1]] ) função pode ser determinado.
A função será diferente de zero somente se o conjunto de valores é k[i] uma permutação de {0,1,...,(n-1)}, e terá uma magnitude de um se este for diferente de zero, eo sinal será o número de swaps necessária para converter a permutação para o permutação identidade.
O valor de uma (n * n) determinante de uma matriz (A) pode ser computada a soma dos produtos de forma (sgn(p) * A[0,p[0]] * A[1,p[1]] * ... * A[(n-1),p[n-1]]) para todos e cada permutação (p) dos números {0,1,2,...,(n-1)}.
O termo (sgn(p)) denota a “assinatura” (ou contar swap) de permutação (p), onde (sgn(p)) é igual a (+1) se (p) é “mesmo” uma “permutação,” e é igual a (-1) se (p) é uma “permutação ímpar.”
Porque há (n!) permutações dos números {0,1,2,...,(n-1)}, esta fórmula tem (n!) summands.
Outros procedimentos de computação o determinante (exemplos: as que envolvem Gaussian eliminação, ou LU factoring, expansão ou por menores, etc) computar implicitamente uma hierarquia das quantidades intermediárias que permitem que os procedimentos para calcular o determinante de uma ordem de (n^3) etapas.
Se o determinante de uma matriz é zero, a matriz não tem um inverso.
Se o determinante de uma matriz é diferente de zero, a matriz tem uma inversa.
Se uma matriz representa os coeficientes de um sistema de equações lineares, e que o determinante da matriz é zero, então o sistema de equações não tem uma solução única.
Se uma matriz representa os coeficientes de um sistema de equações lineares, e que o determinante da matriz é diferente de zero, então o sistema de equações tem uma solução única.
O código inclui procedimentos a seguir para calcular o determinante de qualquer matriz quadrada.
Determinantes para 1*1, 2*2, 3*3, e 4*4 matrizes são calculados pela explícita fórmulas.
Usando explícita fórmulas para calcular os determinantes para as pequenas matrizes é provável que a computar os resultados mais rápido do que usando um procedimento geral para calcular os determinantes.
O procedimento para o caso geral usa LU factoring, mas existem muitos outros métodos de computação 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)-dimensional vectores: cruzada produto
“cross produto:” Uma função associada a um espaço (d)-dimensional que converte uma lista ordenada de (d-1) (d)-dimensional vetores, { A(0), A(1), ..., A(d-2) }, para uma (d)-dimensional vector, de tal forma que a função tem os atributos abaixo:
(1) Se A(i)=0 (vetor zero), para qualquer 0 <= i <= (d-2) e, em seguida, Cross( A(0), A(1), ..., A(d-2) )=0 (vetor zero).
A cruz produto de uma lista ordenada de vetores é zero se algum vetor nessa lista é zero;
(2) Se qualquer par de vetores, A(i) e A(j), para qualquer (i!=j) com 0 <= i <= (d-2) e 0 <= j <= (d-2), (Abs(Dot(A(i),A(j))) = Length(A(i)) * Length(A(j))) são paralelas e, em seguida, Cross( A(0), A(1), ..., A(d-2) )=0 (vetor zero).
A cruz produto de uma lista ordenada de vetores é zero se houver dois vectores nessa lista são paralelos;
(3) Se cada vetor A(i) não é zero, e se não for A(i) paralelo com A(j) para todos os (i!=j) com 0 <= i <= (d-2) e 0 <= j <= (d-2) e, em seguida, B = Cross( A(0), A(1), ..., A(d-2) ) é tal que Dot(A(i),B)=0 para todos os 0 <= i <= (d-2).
Se o produto de uma cruz ordenados lista de vetores é diferente de zero, então o resultado é produto transversal perpendicular a cada vetor na lista ordenada de vectores;
(4) Se (c) é um valor numérico e, em seguida, 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 uma lista ordenada de vetores, a cruz produto da lista ordenada de vectores com qualquer um vetor multiplicado por um valor numérico específico é equivalente à multiplicação do produto cruzada da lista de vetores ordenados pelo valor numérico;
(5) Dado dois números (i) e (j), onde (i!=j) e 0 <= i <= (d-2) e 0 <= j <= (d-2), e uma lista ordenada de (d-1) (d)-dimensional vetores, { A(0), A(1), ..., A(d-2) }, e uma outra lista ordenada de (d-1) (d)-dimensional vetores, { B(0), B(1), ..., B(d-2) }, de tal forma que para todos os ((k != i) && (k != j)) com B(k) = A(k) 0 <= k <= (d-2), e tal que B(i) = A(j) e B(j) = A(i) e, em seguida, Cross( A(0), A(1), ..., A(d-2) ) = (-1) * Cross( B(0), B(1), ..., B(d-2) ).
Se o produto de uma cruz ordenados lista de vetores é diferente de zero e, em seguida, o produto tem que atravessar de sinal contrário, a cruz de produto da mesma lista ordenada de vetores, com excepção de exatamente dois vetores são trocadas (trocados).
Dado e(k) vetores base, para 0 <= k <= (d-1), (d)-dimensional do espaço, a cruz de um produto encomendado (d-1) lista de vetores pode ser computada a computação o determinante da matriz (d*d) abaixo:
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]
O determinante da matriz que não será um número, mas em vez irá ser um vetor, onde haverá coeficientes numéricos para a base (d) vetores.
Por duas dimensões espaciais, eo vetor A = ( ax, ay ):
Cross( A )
= Determinant
(
e(x), e(y),
ax, ay
)
= ay * e(x)
- ax * e(y)
= ( ay, -ax )
O resultado vetor é perpendicular à única entrada de vetores.
Se o avião x-y é visto que tais aumentos de x um rightward direção e y aumenta em um sentido ascendente e, em seguida, o resultado é um produto cruzado de quarto de volta “no sentido horário” em relação à entrada vetor.
Exemplos de produtos cruzados, em duas dimensões espaciais:
Cross( e(x) ) = (-1) * e(y)
Cross( e(y) ) = e(x)
Para espaço tridimensional, e os vectores A = ( ax, ay, az ) e 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) )
Exemplos de produtos cruzados no espaço 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)
Para a quatro dimensões espaciais, e os vectores A = ( ax, ay, az, aw ), B = ( bx, by, bz, bw ), e 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) )
Exemplos de produtos cruzados, em quatro dimensões espaciais:
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)
O código abaixo calcula a cross-produto de dois ou mais vetores.
O teste exemplos mostram que o resultado não-zero para todos os vectores são perpendiculares à entrada de vetores.
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:” Uma matriz quadrada relacionadas com a outra praça matriz do mesmo tamanho tal que o produto das duas matrizes é a “matriz identidade.”
Se uma matriz tem uma matriz inversa associados, só existe uma tal matriz inversa.
Se o determinante de uma matriz é zero, então a matriz é chamado “singular,” e de a matriz não tem uma matriz inversa associados.
O código abaixo calcula o inverso de qualquer quadrado, não-singular matriz.
Inverses para 1*1, 2*2, 3*3, e 4*4 matrizes são calculados pela explícita fórmulas.
Usando explícita fórmulas poderia ser mais rápido do que usando um procedimento geral.
As fórmulas também são explícitas fiáveis.
O procedimento para o caso geral usa LU factoring, mas existem muitos outros métodos de computação inverses de matrizes.
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 produto de uma matriz e um vetor
O multiplicativo (M) produto de uma matriz e um vector (a) é um vetor (b); (M)*(a) = (b).
Matrix (M) “transforma” vetor (a) a (b) vetor.
O código abaixo calcula a multiplicativo produto de uma matriz e um 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 Usando matrizes para traduzir homogênea homogênea vetores
“Tradução:” Uma operação T(s), especificados por um vetor (s) particular, que converte qualquer vector (r) correspondente a um vetor (r + s); T(s) * (r) = (r + s).
T(-s) é a inversa da T(s); Inverse(T(s)) = T(-s); T(-s) * (r) = (r - s).
Uma tradução operação pode ser representada por uma matriz, de tal forma que a multiplicação da tradução matriz por um vetor tem o efeito de traduzir esse vetor.
As contribuições de uma tradução para o funcionamento dos componentes do vetor resultado não dependem, directa ou indirectamente, sobre os componentes do vetor original.
Assim, uma tradução matriz deve ser estruturado de uma forma que lhe permita contribuir para o resultado de uma tradução sem ser afectada, directa ou indirectamente, pelos componentes do vetor original.
Um “vetor” é um vetor “homogêneo” com um elemento igual a 1 final.
A “matriz” é uma matriz “homogénea” com uma linha final com todas as entradas iguais a 0 à excepção de uma entrada de 1 final.
A multiplicação de um produto homogêneo e outra matriz é uma matriz homogênea homogênea matriz.
A multiplicação de um produto homogêneo e de uma matriz é um vetor homogêneo vetor homogêneo.
Um atributo interessante deste produto é que as entradas na última coluna da matriz homogênea, com exceção da última entrada, são adicionados diretamente para os componentes do vetor resultado, sem ser afectada, directa ou indirectamente, pelos componentes do o vetor original.
Assim, uma tradução operação pode ser representada por uma matriz homogénea, e de uma tal matriz pode ser utilizado para traduzir homogênea vetores.
Se queremos representar posições e direções em (d)-dimensional espaço, e nós queremos usar matrizes para traduzir posições e, em seguida, homogênea e homogênea vetores matrizes podem ser usados.
No entanto, a homogeneidade vetores devem ter (d+1) componentes, e as matrizes homogéneas devem ter ((d+1)*(d+1)) entradas.
O último componente de um vetor homogêneo e da última fila de uma matriz homogênea que não correspondem a nenhum dos (d) coordenadas que queremos representar em (d)-dimensional espaço, mas sim servir apenas para que se possa obter matrizes de fazer algumas operações independentemente do (d) coordenar valores.
Usando homogênea homogênea vetores e matrizes facilita a compostos transformações e aplicar transformações a vetores, mas outras operações tornar-se mais incómodo.
Representando posições ou nas direcções (d)-dimensional usando o espaço homogéneo vetores que significa que os vetores terá (d+1) componentes, com o último elemento igual a 1.
Computing o produto escalar de dois vetores homogénea usando um produto regular dot-função seria combinar implicitamente o produto escalar do primeiro (d) componentes com um valor de 1 irrelevante devido ao final de componentes 1.
Por isso, para obter o produto escalar da respectiva (d) componentes do vetor homogêneo (com (d+1) componentes), o produto escalar de vetores homogéneos deverá ser feito utilizando um “vector” especial “homogênea produto escalar” função que simplesmente ignora o final da componente de vetores.
Alternativamente, o produto escalar pode ser computada com o produto escalar função regular com o entendimento de que uma quantidade de 1 deve ser deduzida do resultado antes de interpretar o resultado como pertencente ao (d) relevantes coordenadas.
O código abaixo mostra como uma matriz podem ser formadas que podem então ser usados para traduzir um vetor.
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)-dimensional espaço: rotações
“Rotação em (d)-dimensional espaço:” Uma operação (R) especificado por um plano bidimensional rotação (S), uma rotação (t) ângulo, rotação e um ponto central (c), que converte qualquer ponto (p) correspondente a um ponto (q), de tal forma que:
(1) Length( q - c ) = Length( p - c ), a distância entre o novo ponto a rotação eo ponto central é o mesmo que a distância entre o ponto inicial de rotação e de centro.
Assim, o novo ponto é sobre a superfície de uma esfera particular (d)-dimensional;
(2) Let (P) ser um avião que é paralela à avião rotação (S) e contém ponto (p).
Ponto (q) será (P) no avião.
Assim, o novo ponto é sobre um determinado plano bidimensional no espaço (d)-dimensional;
(3) Dot( (q - c), (p - c) ) = Length( q - c ) * Length( p - c ) * Cos( t ); o produto escalar do vector a partir do centro de rotação do ponto inicial e do vetor de rotação do ponto central para o novo ponto é igual ao produto do comprimento dos dois vectores e co-seno da rotação ângulo.
Assim, o novo ponto é sobre a superfície de uma determinada (d)-dimensional cone;
Constrangimento (1) limites do novo ponto na superfície de uma esfera particular (d)-dimensional.
Condição (2) ainda mais o novo ponto para um determinado plano bidimensional no espaço (d)-dimensional.
Por isso, o novo ponto deve estar em um determinado círculo bidimensional no espaço (d)-dimensional.
Condição (3) ainda mais o novo ponto para um particular (d) cone-dimensional.
Por isso, o novo ponto deve estar em um ponto específico (para Cos( t ) = 1 ou Cos( t ) = (-1)) ou em um de dois possíveis pontos (para (-1) < Cos( t ) < 1).
Uma última limitação é necessária para escolher um dos dois possíveis pontos que resultam das restrições acima mencionadas.
Deixe T(s) representa uma “tradução” operação, onde (s) é um vector arbitrário, de tal forma que T(s) * (r) = (r + s), e Inverse(T(s)) * (r) = (r - s).
A rotação R(S,t,c), por avião S, t ângulo, e ponto central c, pode ser expressa em termos de tradução operações e uma rotação com um ponto central na origem: R(S,t,c) = T(c) * R(S,t,0) * T(-c) = T(c) * R(S,t,0) * Inverse(T(c)).
Deixe R(S,t) representa uma operação com uma rotação rotação ponto central na origem do (d)-dimensional espaço.
Assim, qualquer arbitrariedade rotação pode ser expresso como: R(S,t,c) = T(c) * R(S,t) * T(-c) = T(c) * R(S,t) * Inverse(T(c)).
Vamos considerar duas dimensões rotação aviões (S) que incluem a origem de (d)-dimensional espaço e estão em paralelo com dois eixos coordenar.
Deixe (a) e (b) ser inteiros que representam coordenar eixos em (d)-dimensional espaço, de tal forma que 0 <= a <= (d-1), e 0 <= b <= (d-1), e (a != b).
Nós definimos R(a,b,t) a ser uma rotação sobre a origem dos (d)-dimensional espaço, com uma rotação plano paralelo coordenar eixos indicado por coordenar índices (a) e (b), e com um ângulo (t), de tal forma que R(a,b,(pi/2)) iria converter um valor positivo para coordenar a coordenar Eixo indicado por (a) para coordenar um valor positivo para coordenar o eixo indicado por (b); exemplo: R(a,b,(pi/2)) * e(a) = e(b), onde e(k) indica uma unidade base vetor correspondente ao eixo coordenar (k).
Esta definição proporciona o final desambiguação necessário para determinar o novo ponto totalmente produzido por uma rotação.
Uma rotação que se move em direção a pontos ao longo eixo a eixo b de acordo com um ângulo t é equivalente a uma rotação que se move em direção a pontos ao longo eixo b eixo a por um ângulo (-t); ou seja, R(a,b,t) = R(b,a,(-t)).
Se (d >= 2) e, em seguida, (d)-dimensional espaço tem ((d*(d-1))/2) pares de coordenar eixos.
Assim, para (d >= 2), (d)-dimensional espaço tem ((d*(d-1))/2) distintas eixo-alinhados rotação aviões.
(Exemplos: d=2: (x-y avião); d=3: (x-y avião, y-z avião, avião z-x); d=4: (x-y avião, y-z avião, z-x avião, x-w avião, z-w avião, avião w-y).)
O código abaixo mostra como uma matriz pode ser usado para representar a rotação R(a,b,t).
Um exemplo é mostrado que usa uma matriz de rotação de forma homogénea para rodar um vetor homogêneo para um novo vetor 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 Definição proposta de “counterclockwise”
No espaço bidimensional, o quarto de volta rotação que pode converter o vector e(0) a e(1) é considerado counterclockwise.
A rotação é: R(0,1,(pi/2)) * e(0) = e(1).
Por isso, a rotação R(0,1,t) é counterclockwise.
No espaço tridimensional, o quarto de volta rotação que pode converter o vector e(0) a e(1) é considerado counterclockwise.
O quarto de volta rotação que pode converter o vector e(1) a e(2) também é considerado counterclockwise.
O quarto de volta rotação que pode converter o vector e(2) a e(0) também é considerado counterclockwise.
Por conseguinte, as rotações { R(0,1,t), R(1,2,t), R(2,0,t) } são counterclockwise.
A fórmula booleana abaixo se encaixa na definição de counterclockwise para bidimensionais e tridimensionais espaciais, e extrapola o padrão para o caso geral da rotação R(a,b,t):
bool counterClockwise =
(
( (a < b) && (((a + b) % 2) == 1) )
|| ( (a > b) && (((a + b) % 2) == 0) )
);
Esta regra seria definir as rotações “counterclockwise” como abaixo:
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)-dimensional espaço: no sentido horário rotações
Usando a definição proposta de “counterclockwise” na secção anterior, as rotações “no sentido horário” a seguir são:
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 exemplo, em quatro dimensões espaciais, a rotação seis matrizes abaixo representam rotações no sentido horário:
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)-dimensional espaço: número total de eixo-alinhados orientações
O número de eixos distintos-alinhados nas orientações (d)-dimensional espaço é:
distinctOrientations = Factorial( d ) * IntPow( 2, (d-1) );
Um argumento apoiando esta fórmula aparece nos parágrafos seguintes.
Experimental provas verificar a fórmula também é apresentada a seguir.
Deixe { g(0), g(1), ..., g(d-1) } ser um conjunto de (d) mutuamente-vetores unitários perpendiculares anexado a um objeto no espaço (d)-dimensional.
Estes vectores unitários indicar a orientação do objeto em relação a coordenar os eixos da (d)-dimensional espaço.
Deixe o “padrão” de “orientação” a objeto g(i) = e(i) ser tal que, para todos os 0 <= i <= (d-1); cada um vetor de orientação para o objeto está associado a uma unidade distinta de um vector base no sistema de coordenadas (d)-dimensional espaço.
O primeiro vector de orientação para o objeto pode ser paralelo com qualquer um dos eixos (d) do (d)-dimensional espaço, e pode ser em qualquer uma das duas direções em relação a esse eixo.
Assim, a primeira orientação vetor tem (2 * d) opções.
O segundo vector de orientação para o objeto pode ser paralela a qualquer dos restantes (d-1) eixos, com duas direções possíveis em relação a esse eixo.
Assim, a segunda orientação vetor tem (2 * (d-1)) opções.
Do mesmo modo, orientação para o próximo vetor tem (2 * (d-2)) opções.
O último vector é forçado a orientação para o final restante eixo, e é forçado a uma das duas possíveis direções ao longo desse eixo (para conservar a paridade de orientação para o conjunto de vetores; o determinante de um conjunto ordenado de vetores a orientação deve ser constante ).
Assim, o vetor tem apenas uma orientação final (1) opção.
O número total de combinações é dado pelo produto: ((2 * d) * (2 * (d-1)) * ... * (2 * 3) * (2 * 2) * (1)) = (2^(d-1)) * (d*(d-1)*(d-2)*...*2*1) = (2^(d-1)) * Factorial(d).