Algebra liniara
Colin Fahey
1. Software
LinearAlgebra.zip
Algebra liniara codul sursă (C#)
19910 bytes
MD5: 11d8c8035cac30ba543e5e0b72ee9767
2. Introducere
Acest articol descrie vectori şi matrici (d)-dimensionale în spaţiu.
3. (d)-dimensional space: Atribute
3.1 Array
"matrice:" O colectie de variabile, astfel încât fiecare variabilă are un nume unic, şi de asemenea că numele poate fi atribuit un ordin.
Valori pot fi utilizate ca nume de variabile într-o matrice.
De exemplu, dacă un array (d) conţine variabile, atunci valori { 0, 1, 2, ..., (d-1) } poate fi numele atribuit de variabile în matrice.
3.2 (d)-dimensional vector
"(d)-dimensional vector:" un array de (d) variabile.
"vector pe componente:" o variabilă într-un vector.
3.3 (d)-dimensional vector spaţiu
"o-dimensional space:" set complet de valori care pot fi stocate de către o variabilă.
"(d)-dimensional space:" set complet de combinaţii de valori care pot fi stocate de către o serie de variabile (d).
Formale, definirea unui "spatiu vectorial:"
Să (T) fi o bază de tip (exemple: numărul real, întreg, numărul de complexe, număr raţional, etc.)
Orice variabila de un tip de bază este numit un "scalar."
Un "tip (T)-(d)-dimensional vector spaţiu" este un set (S) de (d)-dimensional vectori având două operaţii, plus (+) vector, scalar şi multiplicare (*), care îndeplinesc condiţiile de mai jos:
(1) Dacă (v) şi (w) sunt doi vectori în orice (S), apoi (v + w) este, de asemenea, un vector în (S);
(2) Dacă (u), (v), şi orice (w) sunt trei vectori în (S), apoi (u + v) + w = u + (v + w);
[aditiv commutativity]
(3) Dacă (v) şi (w) sunt doi vectori în orice (S), apoi (v + w) = (w + v);
[aditiv associativity]
(4) Există un "vector de zero," (0), în (S), astfel că pentru orice vector (v) în (S), (v + (0)) = v;
[aditiv de identitate]
(5) Dacă este (c) orice tip de scalare (T), şi este (v) orice vector în (S), apoi produsul (c * v) este un vector în (S);
(6) Dacă (a), (b), şi (c) sunt scalars orice tip de (T), şi (v) şi (w) sunt vectori în orice (S), apoi (a + b) * v = a*v + b*v, şi c*(v + w) = c*v + c*w;
[multiplicativ distributivity]
(7) Dacă (a) şi (b) sunt scalars orice tip de (T), şi este (v) orice vector în (S), apoi (a*b)*v = a*(b*v);
(8) Dacă "1" este un tip de scalare (T) astfel că (1*1)=1, şi este (v) orice vector în (S), apoi (1*v) = v;
(9) Pentru fiecare vector (v) în (S), vectorul (-1)*v = -v îndeplineşte v + (-v) = (0);
[aditiv invers]
3.4 (d)-dimensional vector de cod
Codul de mai jos arată modul în care o (d)-dimensional vector, cu 64 de biţi în virgulă flotantă punct de componente, pot fi puse în aplicare.
O matrice este suficientă pentru a reprezenta un vector.
Codul de mai jos are un array conţinute într-o clasă, doar pentru comoditate.
Acest cod nu este destinat a fi eficiente.
O structură (de exemplu: "struct", un tip de valoare) reprezintă vectori de un număr fix de dimensiuni (exemple: 3 sau 4) este probabil să fie mult mai eficient decât generale (d)-dimensional vector clasa arătat aici.
Desi codul de mai jos defineşte un vector cu plutitor-punct de componente, acest document, de asemenea, face uz de vectori cu componente număr întreg.
Codul de mai jos pot fi uşor modificate pentru punerea în aplicare a vectorilor cu număr întreg de componente.
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:" Un vector cu toate componentele (d) egal cu 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 spatiu: puncte
"(d)-dimensional spaţiu punct:" un array de (d) specifice, cu valori variabile "(valori de coordonate)."
"(d)-dimensional spaţiu origine:" un array de (d) toate variabilele cu valori egale cu zero.
3.6 (d)-dimensional vectori: non-relative şi relativă
Un (d)-dimensional vector care "nu" este "relativă" este o (d)-dimensional vector care, direct sau reprezinta o stare de configurare.
Un (d)-dimensional vector care este "relativă" este o (d)-dimensional, care reprezinta vectorul de modificări la un set de componente.
Un vector relativă poate reprezenta diferenţa dintre doi vectori nu-relativă.
Având în vedere o ruda vectorială, determinarea unei non-relativ de statutul sau de configurare, folosind ca vector relativă, care necesită combinarea cu un vector relativ non-vectorială relativă.
Non-relativă şi relativă vectorilor vectorii sunt atât vectori.
Dacă un anumit vector nu este relativă sau relativă trebuie să fie specificat, atunci când vectorul este definit.
Dacă un (d)-dimensional vector este interpretat ca fiind non-relativă, atunci (d)-dimensional vector poate reprezenta un punct în spaţiu (d)-dimensional.
3.7 (d)-dimensional vector plus, Navigare, şi scalarea
Vector plus, Navigare, şi scalarea:
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 baza vectorilor
Formale, definirea unui "baza" unui spatiu vectorial:
Să (T) fi o bază de tip (exemple: numărul real, întreg, numărul de complexe, număr raţional, etc.)
Orice variabila de un tip de bază este numit un "scalar."
Să (V) fi un "(T) de tip (d)-dimensional vectorul spatiu."
În cazul în care non-zero vectori { u1, u2, ..., ud } în (V) sunt de aşa natură încât în fiecare vector (v) în (V) poate fi scrisă ca o "combinaţie liniară" a acestor vectori, v = c1*u1 + c2*u2 + ... + cd*ud, în cazul în care { c1, c2, ..., cd } sunt scalars de tip (T), apoi (V) este "spanned" prin vectori { u1, u2, ..., ud }.
Orice set de non-zero vectori { u1, u2, ..., ud } că "span" spatiu vectorial (V) este numit pe "baza" de (V).
Un simplu "baza" unui (d)-dimensional vector spaţiu este un set de (d) distincte (d)-dimensional vectori, fiecare cu o componentă egală cu unul şi toate celelalte componente egale cu zero.
Astfel de baza vectorilor sunt "orthonormal," ceea ce înseamnă că acestea sunt reciproc "perpendiculare-(orthogonal)" şi că fiecare unitate a vectorului lungime.
Fiecare astfel de unitate vectorul este un vector paralel cu una dintre axele de coordonate (d).
Exprimându-un vector arbitrar ca o "combinaţie liniară" a acestor vectori este de baza directe; fiecare componentă a vectorului este arbitrară, înmulţită cu vectorul baza corespunzătoare, şi aceste produse sunt adăugate împreună pentru a forma de vector arbitrar.
Codul de mai jos arată cum un vector poate fi exprimat ca o "combinaţie liniară de" vectorii de baza.
Codul de mai jos în mod arbitrar orthonormal defineşte un set de baza vectorilor.
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 )
// . . .
}
}
Orice vector în spaţiul (d)-dimensionale pot fi exprimate ca o sumă de produse de baza si numere de vectori:
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 spaţiu: distanţa dintre punctele
Să (P) fi un (d)-dimensional vector care reprezintă un punct în spaţiu (d)-dimensional.
Să (Q) fi un (d)-dimensional vector care reprezintă un punct în spaţiu (d)-dimensional.
Să (R) fi un (d)-dimensional, care reprezinta vectorul de schimbare în (d) coordonate pentru a obţine de la punctul de a (P) punct (Q); R = (Q - P).
Într-un spaţiu-dimensional, P = ( p0 ), Q = ( q0 ), şi R = (Q - P) = ( q0 - p0 ).
Distanţa dintre cele două puncte este: Abs( q0 - p0 ).
În spaţiul bidimensional, P = ( p0, p1 ), Q = ( q0, q1 ), şi R = (Q - P) = ( q0-p0, q1-p1 ).
Interpretarea celor doua perpendiculare deplasarilor ca perpendicular laturile unui triunghi din dreapta, distanţa între puncte corespunde cu lungimea de hypotenuse din acest triunghi.
Pythagorean de formulă, (a*a) + (b*b) = (c*c), în cazul în care (a) şi (b) sunt perpendiculare de lungimi de laturi a unui triunghi din dreapta, şi (c) este lungimea de hypotenuse (skew parte), pot fi utilizate pentru a determina distanţa dintre cele două puncte: Sqrt( Sq(q0-p0) + Sq(q1-p1) ).
În tri-dimensională de spatiu, P = ( p0, p1, p2 ), Q = ( q0, q1, q2 ), şi R = (Q - P) = ( q0-p0, q1-p1, q2-p2 ).
Interpretarea deplasament ( q0-p0, q1-p1, 0 ) ca perpendicular laturile unui triunghi din dreapta, şi folosind formula de Pythagorean, dă distanţa dintre punctul (P) şi punctul ( q0, q1, p2 ): d01 = Sqrt( Sq(q0-p0) + Sq(q1-p1) ).
Deplasarea ( q0-p0, q1-p1, 0 ) este perpendicular pe ( 0, 0, q2-p2 ) deplasare, precum şi a unui alt drept-triunghi se pot forma, precum şi Pythagorean formulă poate fi utilizat din nou.
Astfel, distanta de la punctul de a (P) punct (Q) este data de: 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) ).
Metoda de extindere a formulei de la distanta bidimensional spaţiu pentru a tri-dimensională de spaţiu să poată fi aplicate în mod repetat pentru a determina distanţa formulă pentru (d)-dimensional spaţiu: Sqrt( (Sq(q0-p0) + Sq(q1-p1)) + Sq(q2-p2) + ... + Sq(qd-pd) ).
Codul de mai jos defineşte o funcţie numită "Lungime" care computes lungimea unei (d)-dimensional vector.
Atunci când un vector reprezintă deplasarea între două puncte în spaţiu (d)-dimensional, lungimea vectorului reprezintă distanţa dintre cele două puncte.
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 vectori: dot produs
"Punct" de "produs" converteşte doi (d)-dimensională la un număr de vectori.
Codul de mai jos computes punct un produs de doi vectori:
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);
}
// . . .
}
Pentru orice vector, (A), Length(A) = Sqrt(Dot(A,A)).
3.11 (d)-dimensional vectori: definiţie a "paralel"
Vectorilor (A) şi (B) sunt "paralele," dacă toate condiţiile următoare sunt adevărate:
(1) A.Length() > 0;
(2) B.Length() > 0;
(3) Abs(Dot(A,B)) = A.Length()*B.Length().
Codul de mai jos determină dacă o pereche de vectori sunt paralele (eventual anti-aliniate).
Plutitor-punct de numere pot acumula erori in parte din cauza fragmente limitate de precizie, şi, prin urmare, computerul ar trebui să includă codul de non-zero, atunci când toleranţele compararea plutitoare-punct de numere.
Codul include valorile exemplu de toleranţă, ci toleranţă de exemplu, valorile ar putea să nu fie potrivite pentru anumite sarcini.
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 vectori: definiţie a "perpendicular"
Vectorilor (A) şi (B) sunt "perpendiculare" dacă toate condiţiile următoare sunt adevărate:
(1) A.Length() > 0;
(2) B.Length() > 0;
(3) Abs(Dot(A,B)) = 0.
Codul de mai jos determină dacă o pereche de vectori sunt perpendiculare. Plutitor-punct de numere pot acumula erori in parte din cauza fragmente limitate de precizie, şi, prin urmare, computerul ar trebui să includă codul de non-zero, atunci când toleranţele compararea plutitoare-punct de numere.
Codul include valorile exemplu de toleranţă, ci toleranţă de exemplu, valorile ar putea să nu fie potrivite pentru anumite sarcini.
public class VectorF64
{
// . . .
public static bool Perpendicular( VectorF64 a, VectorF64 b )
{
// Smallest normalized float : 1.1754943e-38
// Smallest normalized double: 2.2250738585072020e-308
double nonZeroThreshold = 1.0e-38; // conservative for double
// double: (52+1)-bit mantissa; log10(2^53)=15.95 decimal digits
double fractionalDifferenceThreshold = 1.0e-14; // conservative
if ((null == a) || (null == b))
{
return (false); // Vector is not specified.
}
if ((null == a.components) || (null == b.components))
{
return (false); // Vector is empty.
}
if (a.Dimensions() != b.Dimensions())
{
return (false); // Vectors not the same size.
}
double lengthA = a.Length();
if (lengthA <= nonZeroThreshold)
{
return (false);
}
double lengthB = b.Length();
if (lengthB <= nonZeroThreshold)
{
return (false);
}
double zeroImpliesPerpendicular =
Math.Abs( VectorF64.Dot( a, b ) ) / (lengthA * lengthB);
double absoluteDifferenceFromZero =
Math.Abs( zeroImpliesPerpendicular );
if (absoluteDifferenceFromZero <= fractionalDifferenceThreshold)
{
return (true);
}
return (false);
}
// . . .
public static void Test()
{
// . . .
// Example of testing for perpendicular vectors:
// A 6-dimensional vector:
VectorF64 vf2 = new VectorF64( 0.0, 1.0, 2.0, 0.0, 4.0, 5.0 );
vf2.WriteLine(); // ( 0, 1, 2, 0, 4, 5 )
// A 6-dimensional vector:
VectorF64 vg2 = new VectorF64( 10.0, 0.0, 0.0, -5.0, 0.0, 0.0 );
vg2.WriteLine(); // ( 10, 0, 0, -5, 0, 0 )
// Determine if the specified vectors are perpendicular
bool perpendicular = VectorF64.Perpendicular( vf2, vg2 );
Log.WriteLine( perpendicular ); // True
// Add a non-negligible displacement to a component of a vector:
vf2[0] += 1.0e-13;
vf2.WriteLine(); // ( 1E-13, 1, 2, 0, 4, 5 )
// Determine if the specified vectors are perpendicular
perpendicular = VectorF64.Perpendicular( vf2, vg2 );
Log.WriteLine( perpendicular ); // False
// . . .
}
}
3.13 Matrices
"matrice:" O colectie de variabile, astfel încât fiecare variabilă are o combinaţie unică a unui nume de "rând" şi o "coloană" nume.
"intrare:" o variabilă într-o matrice.
Valori pot fi utilizate ca nume de "rând" şi de "coloană" pentru numele de variabile într-o matrice.
De exemplu, dacă o matrice a (totalRows) rânduri şi de coloane (totalColumns), atunci valori { 0, 1, ..., (totalRows-1) } poate fi numele atribuit de rânduri, şi de valori { 0, 1, ..., (totalColumns-1) } poate fi numele atribuit de coloane.
Astfel, o variabile într-o matrice pot fi specificate prin specificarea unei perechi de întregi, ( row, column ), indicând combinaţie a rândul şi coloana corespunzătoare pentru a indicilor specifice de variabilă.
Dimensiunea unei matrice este specificat ca "(totalRows) * (totalColumns)" (sau "(totalRows) de (totalColumns))."
Acest ordin de dimensiuni este aceeaşi cu ordinea de dimensiuni utilizate pentru a specifica intrări în matrice "(( row, column ))."
[Această convenţie este oarecum unfortunate, deoarece pentru multe utilizări bidimensional (exemple: imagini, grafice, etc) ca de obicei de convenţie este de a specifica dimensiunile "width * height" ca şi coordonate ca "( horizontal, vertical )" (sau "( x, y ))."
Aceasta este opusul ordinea de dimensiuni şi coordonate folosit pentru a descrie matrici şi intrările lor. ]
O matrice cu (totalRows) egală cu (totalColumns) este numit "pătrat;" altfel, este numit matrice "dreptunghiular."
O matrice poate fi privită ca conţinând un set de "vectori rând," în cazul în care variabilele în fiecare rând sunt interpretate ca aparţinând unui vector.
O matrice poate fi, de asemenea, considerate ca fiind conţine un set de "vectori coloană," în cazul în care variabilele în fiecare coloană sunt interpretate ca aparţinând unui vector.
Matrices poate reprezenta o mare varietate de relaţii matematice.
În sensul unei matrice, precum şi operaţiunile care ar putea fi potrivite pentru prelucrare, inregistrarile de o matrice, depinde de context.
Cu toate acestea, există reguli de baza de matrice aritmetic care sunt relevante pentru mai multe contexte, iar aceste norme de bază vor fi definite într-o secţiune ulterioară.
Un array si o valoare (totalColumns) sunt suficiente pentru a reprezenta o matrice.
De matrice poate avea (totalRows * totalColumns) variabile, de la intrarea şi la ( row, column ) poate să corespundă cu matrice variabilă de la indexul ((totalColumns * row) + column).
Codul de mai jos defineşte o matrice, cu 64 de biţi în virgulă flotantă punct de intrări.
Un array si o valoare (totalColumns) sunt suficiente pentru a reprezenta o matrice.
Codul de mai jos are un array si o valoare (totalColumns) conţinute într-o clasă, doar pentru comoditate.
Codul de mai jos nu este destinat a fi eficiente.
O structură (de exemplu: "struct", un tip de valoare) reprezintă matrici de dimensiuni special fix (exemple: 2*2, 3*3, sau 4*4) este probabil să fie mult mai eficient decât generale (totalRows * totalColumns) clasa arătat aici.
Desi codul de mai jos defineşte o matrice cu plutitor-punct de intrări, de asemenea, acest articol face uz de matrici cu intrări număr întreg.
Codul de mai jos pot fi usor modificate pentru punerea în aplicare a matrici cu intrări număr întreg.
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 matrice:" O matrice cu toate intrările egal cu 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 plus, Navigare, multiplicare şi
Matrix plus, Navigare, şi multiplicare.
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 ]
// . . .
}
}
"identitate matrice:" O matrice pătrată (totală este egală total de rânduri de coloane) cu intrări egal cu (1) pe diagonală (intrări, cu un rând index egală cu coloana indexul său), şi cu toate celelalte intrări egal cu (0).
Dacă o matrice patrata, (M), este înmulţită cu o "matrice de identitate," (I), de acelaşi număr de rânduri (şi coloane), produsul este egal cu (M).
Multiplying (I) de (M) de asemenea, produce un produs egal cu (M).
Astfel, "matricea de identitate" este similară cu numărul "1" pentru multiplicare de numere (scalars).
Codul de mai jos creează o identitate matrice cu un anumit număr de rânduri.
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 de factoring; inapoi-de substituţie
"Matrix LU de factoring:" o procedura care sa converteasca un pătrat cu matrice, (M), pentru a doua noi matrici pătrate, (L) şi (U), având aceeaşi dimensiune ca (M), astfel că (L) * (U) = (M), şi că o astfel de matrice (U) este "superior triunghiular" (toate intrările de mai jos sunt diagonalei zero), şi este "mai mic" (L) matrice "triunghiular" (toate intrările de deasupra diagonalei sunt zero).
Matrix LU de factoring poate fi utilizat ca parte a unui mai mare, procedura de a rezolva un sistem de ecuaţii, sau invers, pentru a găsi o matrice, sau pentru a găsi o matrice de determinant.
Matrix LU de factoring este o alternativă la "Gauss-Jordan procedura de eliminare.
Gauss-Jordan eliminarea necesită un sistem de ecuaţii (A)*(x)=(b), întrucât doar de factoring LU necesită o matrice (A).
De asemenea, dupa stabilirea LU de factoring de o matrice (A), este foarte uşor să se determine (x) dat orice (b).
Procedura pentru a rezolva vectorului (x) în (A)*(x)=(b), dat fiind (b) şi LU de factoring (L)*(U)=(A), implică o procedură care va fi numit "de back-substituţie" aici.
Codul de mai jos include o procedură de calcul a LU factori de orice patrat, non-to matrice.
Codul de mai jos include o procedură de a face back-LU de substituţie.
Atentie: Computerul codul de mai jos computes de L şi U un sir de factori permuted-o anumită versiune de matrice.
De asemenea, L şi U matrice rezultatele sunt combinate într-o singură matrice de ieşire.
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: determinant
"determinant:" o funcţie care transformă orice pătrat (n * n) matrice (A) la un număr, det(A), astfel că:
(1) Dacă o matrice (B) rezultate din schimbul de două rânduri, coloane sau două, de o matrice (A), apoi det(B) = (-(det(A)));
(2) Dacă o matrice (B) rezultate din înmulţirea orice rand, orice coloana sau, de o matrice (A), de către un număr (c), apoi det(B) = c * det(A);
(3) Dacă o matrice (B) rezultate de la un multiplu de a adăuga un rând la un alt rând, sau adăugarea la un multiplu de o coloană într-o altă coloană, apoi det(B) = det(A).
(4) Dacă (A) este o matrice superior triunghiular-(A[i,j]=0 pentru toate (i>j)) sau de o mică matrice triunghiular (A[i,j]=0 pentru toate (i<j)), apoi det(A) = (A[0,0] * A[1,1] * ... * A[n-1,n-1]);
(De exemplu, un determinant de identitate este o matrice; det(I) = 1.)
Reguli (1), (2), şi (3), pot fi utilizate, într-un proces numit "Gaussian de eliminare, pentru a converti orice matrice pătrate într-o matrice triunghiular.
Articolul (4) pot fi utilizate pentru a calcula a determinant a unei matrice triunghiular.
Formale, definirea unui "determinant:"
Să M[n](K) desemna set de toate (n * n) matrici peste domeniul (K).
De "determinant" este unicul funcţia (F) astfel că F:M[n](K) --> K cu două atribute:
(1) F este multilinear alternand cu privire la coloane (sau de rânduri); şi,
(2) F( I ) = 1.
Un "multilinear" funcţia (D) de (n * n) o matrice (A) pot fi scrise ca: 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]] ) ), în cazul în care suma este preluat de toate (n^n) combinaţii de 0 <= k[i] <= (n-1), şi în cazul în care e[j] reprezinta rând (j) de identitate matrice.
Această funcţie poate fi caracterizată de valorile de D( e[k[0]], e[k[1]], ..., e[k[n-1]] ).
O "alternativ multilinear" funcţia este o funcţie multilinear, (D), astfel că D( ..., e[i], ..., e[i], ... ) = 0, şi D( ..., e[i], e[j], ... ) = (-( D( ..., e[j], e[i], ... ) )).
Astfel, schimbarea de coloane (sau de rânduri) în semn de modificări din funcţie.
De asemenea, termenul,, D( e[k[0]], ..., e[k[n-1]] ) este zero dacă k[i] combinaţie de valori nu implică toate valori unice; dacă setul de valori nu este un permutation din numerele de {0,1,...,(n-1)}.
Lăsându-D( I ) = 1 (D( e[0], e[1], ..., e[n-1] ) = 1), în combinaţie cu alternative multilinear atribute descrise mai sus, înseamnă că toate valorile de alternativ multilinear funcţia D( e[k[0]], e[k[1]], ..., e[k[n-1]] ) poate fi determinată.
Funcţia va fi non-zero doar dacă setul de valori este o k[i] permutation de {0,1,...,(n-1)}, şi va avea o magnitudine de una, dacă aceasta este non-zero, şi semn va fi numărul de swap necesare pentru a converti la permutation permutation de identitate.
Valoarea unui determinant de (n * n) o matrice (A) pot fi calculate prin asocierea de produse de forma (sgn(p) * A[0,p[0]] * A[1,p[1]] * ... * A[(n-1),p[n-1]]) pentru fiecare permutation (p) din numerele de {0,1,2,...,(n-1)}.
Termenul desemnează (sgn(p)) "semnătura" (sau swap count) din permutation (p), în cazul în care (sgn(p)) este egală cu (+1) dacă (p) este "chiar" un "permutation," şi este egală cu (-1) dacă (p) este o "permutation nui adevărat."
Pentru că există (n!) permutations din numerele de {0,1,2,...,(n-1)}, această formulă a (n!) summands.
Alte proceduri pentru calculul determinant (exemple: cele care implică eliminarea Gaussian, sau LU de factoring, sau extinderea de către minori, etc), implicit o ierarhie a calcula cantităţile intermediare, care permit ca aceste proceduri de a calcula un determinant în scopul de a (n^3) paşi.
În cazul în care o matrice de determinant este zero, matrice nu are un invers.
În cazul în care determinant este o matrice de non-zero, matricea are o invers.
În cazul în care o reprezinta matricea coeficienţilor de un sistem de ecuaţii liniare, şi faptul că determinant al matricei este zero, atunci sistemul de ecuatii nu au o soluţie unică.
În cazul în care o reprezinta matricea coeficienţilor de un sistem de ecuaţii liniare, şi determinant al matricei care este non-zero, atunci sistemul de ecuatii are o soluţie unică.
Codul de mai jos include proceduri pentru calculul determinant de orice matrice pătrate.
Determinanţi pentru 1*1, 2*2, 3*3, şi 4*4 matrici sunt calculate prin formule explicite.
Folosind formule explicite pentru a calcula de determinanţi pentru matrici mici, este de natură să calculez rezultate mai repede decât utilizând o procedură generală de a calcula aceste determinanţi.
Procedura pentru utilizări generale caz LU de factoring, dar sunt multe alte metode de calcul determinanţi.
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 vectori: Maşină de produs
"Maşină de produs:" o funcţie asociată cu un (d)-dimensional care sa converteasca un spaţiu ordonat lista de (d-1) (d)-dimensional vectori, { A(0), A(1), ..., A(d-2) }, la o (d)-dimensional vector, astfel că funcţia are atributele de mai jos:
(1) Dacă A(i)=0 (vectorul zero), pentru orice 0 <= i <= (d-2), apoi Cross( A(0), A(1), ..., A(d-2) )=0 (zero vector).
Crucea de produs comandat de o listă a vectorilor este zero dacă orice vector în care lista este zero;
(2) În cazul în care orice pereche de vectori, A(i) şi A(j), pentru orice (i!=j) cu 0 <= i <= (d-2) şi 0 <= j <= (d-2), sunt paralele (Abs(Dot(A(i),A(j))) = Length(A(i)) * Length(A(j))), apoi Cross( A(0), A(1), ..., A(d-2) )=0 (zero vector).
Crucea de produs comandat de o listă a vectorilor este zero dacă există doi vectori în această listă sunt paralele;
(3) Dacă fiecare vector A(i) este non-zero, şi, dacă A(i) nu este paralel cu A(j) pentru toate (i!=j) cu 0 <= i <= (d-2) şi 0 <= j <= (d-2), apoi B = Cross( A(0), A(1), ..., A(d-2) ) este de aşa natură încât Dot(A(i),B)=0 pentru toate 0 <= i <= (d-2).
În cazul în care eco produs comandat de o listă a vectorilor este non-zero, atunci crucea produsul rezultat este perpendiculară, pentru fiecare vector de ordonate în lista de vectori;
(4) Dacă (c) este o valoare numerică, apoi 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) ).
Dat fiind comandat o lista de vectori, cu crucea de produse de comandat lista de vectori cu nici un vector înmulţită cu o anumită valoare numerică este echivalentă cu înmulţirea cruce produs din lista ordonată a vectorilor de valoare numerică;
(5) Având în vedere două numere (i) şi (j), în cazul în care (i!=j) şi 0 <= i <= (d-2) şi 0 <= j <= (d-2), şi a ordonat o listă de (d-1) (d)-dimensional vectori, { A(0), A(1), ..., A(d-2) }, şi a ordonat-o altă listă de (d-1) (d)-dimensional vectori, { B(0), B(1), ..., B(d-2) }, astfel că B(k) = A(k) pentru toate ((k != i) && (k != j)) cu 0 <= k <= (d-2), şi că astfel de B(i) = A(j) şi B(j) = A(i), apoi Cross( A(0), A(1), ..., A(d-2) ) = (-1) * Cross( B(0), B(1), ..., B(d-2) ).
În cazul în care eco produs comandat de o listă a vectorilor este non-zero, atunci crucea pe care produsul a opus semn de cruce din acelaşi produs comandat lista de vectori, cu excepţia exact doi vectori în curs de schimb (swapped).
Având în vedere baza vectorilor e(k), pentru 0 <= k <= (d-1), de (d)-dimensional spaţiu, cu crucea de un produs comandat (d-1) lista de vectori pot fi calculate de calcul a determinant al (d*d) matricea de mai jos:
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]
De determinant al matricei care nu va fi un număr, dar în schimb va fi un vector, unde va fi numerică a coeficienţilor de (d) baza vectorilor.
Pentru bidimensional spaţiu, şi vectorul A = ( ax, ay ):
Cross( A )
= Determinant
(
e(x), e(y),
ax, ay
)
= ay * e(x)
- ax * e(y)
= ( ay, -ax )
Rezultatul este vectorul perpendicular pe vectorul unic de intrare.
În cazul în care planul este x-y vizitate, astfel încât x creşte într-o direcţie rightward şi y creşte într-o directie in sus, apoi cruce produsul rezultat este de un sfert de-rândul său, "în sens orar" relativ la vectorul de intrare.
Exemple de produse eco bidimensional în spaţiu:
Cross( e(x) ) = (-1) * e(y)
Cross( e(y) ) = e(x)
Pentru tri-dimensională de spaţiu, şi a vectorilor A = ( ax, ay, az ) şi 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) )
Exemple de produse în cruce tri-dimensională de spaţiu:
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)
Timp de patru-dimensionale spaţiu, şi a vectorilor A = ( ax, ay, az, aw ), B = ( bx, by, bz, bw ), şi 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) )
Exemple de produse cruce în patru-dimensionale spaţiu:
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)
Codul de mai jos computes de eco-produsul a două sau mai multe vectori.
Testul exemple arată că non-zero rezultat vectorilor sunt perpendiculare la spre tot de la vectori de intrare.
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: invers
"matricea invers:" un pătrat cu matrice referitoare la o altă matrice pătrate de aceeaşi dimensiune, astfel încât produsul a două matrici este "matricea identitate."
Dacă o matrice are o matrice asociat invers, nu există decât o singură astfel de matrice inversă.
În cazul în care o matrice de determinant este zero, atunci matricea se numeste "la singular," şi matricea nu au o matrice asociat invers.
Codul de mai jos computes de invers de orice patrat, non-to matrice.
Inverses pentru 1*1, 2*2, 3*3, şi 4*4 matrici sunt calculate prin formule explicite.
Utilizarea formulelor explicite ar putea fi mai repede decât utilizând o procedură generală.
În mod expres de formule sunt, de asemenea, de incredere.
Procedura pentru utilizări generale caz LU de factoring, dar sunt multe alte metode de calcul inverses de matrici.
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 Multiplicativ produs de o matrice şi un vector
Multiplicativ de produs dintr-o matrice (M) şi un vector (a) este un vector (b); (M)*(a) = (b).
Matrix (M) "transforma" vectorul (a) pentru a vectorului (b).
Codul de mai jos computes de multiplicativ produs de o matrice şi un vector.
public class MatrixF64
{
// . . .
public static VectorF64 operator *( MatrixF64 m, VectorF64 v )
{
if (null == m)
{
return (null); // Matrix is not specified.
}
if
(
(null == m.entries)
|| (m.totalRows <= 0)
|| (m.totalColumns <= 0)
)
{
return (null); // Matrix is empty.
}
if (null == v)
{
// Vector is not specified.
return (VectorF64.Zero( m.totalRows ));
}
if (v.Dimensions() != m.totalColumns)
{
// The number of columns of the matrix must be equal
// to the number of rows of the vector.
return (VectorF64.Zero( m.totalRows ));
}
// The result vector has the same number of rows as the matrix.
VectorF64 result = VectorF64.Zero( m.totalRows );
// Component (i) of the result vector is equal to the
// dot product of row (i) of the matrix with the vector.
for (int i = 0; i < m.totalRows; i++)
{
double dotProduct = 0.0;
for (int j = 0; j < m.totalColumns; j++)
{
dotProduct += (m[i, j] * v[j]);
}
result[i] = dotProduct;
}
return (result);
}
// . . .
public static void Test()
{
// . . .
// Examples of a product of a matrix and a vector:
MatrixF64 t3x3 =
new MatrixF64
(
3, 3,
2.0, 5.0, 7.0, // row 0
-4.0, -1.0, 6.0, // row 1
9.0, 8.0, 3.0 // row 2
);
VectorF64 p3x1 = new VectorF64( 1.0, 2.0, 3.0 );
VectorF64 q3x1 = t3x3 * p3x1;
q3x1.WriteLine();
// ( 33, 12, 34 )
// The multiplicative product of the inverse
// of the matrix and the result vector computed
// above should reproduce the original vector.
MatrixF64 inverset3x3 = t3x3.Inverse();
VectorF64 prodinvt3x3andq3x1 = inverset3x3 * q3x1;
prodinvt3x3andq3x1.WriteLine();
// ( 1, 2, 3 )
// This vector is equal to the original vector,
// as expected.
// . . .
}
}
3.20 Utilizarea matrici omogene de a traduce omogene de vectori
"traducere:" o operaţie T(s), specificată de către un anumit vector (s), care converteşte orice vector (r) la un vector (r + s) corespunzătoare; T(s) * (r) = (r + s).
T(-s) este invers de T(s); Inverse(T(s)) = T(-s); T(-s) * (r) = (r - s).
O operaţiune de traducere poate fi reprezentat de o matrice, astfel că înmulţirea de traducere matrice de un vector are ca efect faptul că traducerea vector.
Contribuţiile de o traducere funcţionare a componentelor vectorului de rezultatul nu depinde, direct sau indirect, pe de componente ale vectorului original.
Astfel, o matrice de traducere trebuie să fie structurat într-un mod care îi permite să contribuie la o traducere a rezultat, fără a fi afectate, direct sau indirect, de către componentele vectorului original.
Un "vector omogen" este un vector cu un final componentă egală cu 1.
O "matrice omogenă" este o matrice cu un final rând cu toate intrările egală cu excepţia 0 pentru un final de intrarea 1.
Multiplicativ de produs dintr-o matrice de omogene şi de un alt omogene este o matrice matrice omogenă.
Multiplicativ de produs dintr-o matrice omogena omogene şi un vector este un vector de omogene.
Un interesant atribut al acestui produs este că intrările în final coloană din matricea omogenă, cu excep ț ia finală de intrare, se adaugă direct la componentele vectorului de rezultat, fără a fi afectate, direct sau indirect, de către componentele vectorul original.
Astfel, o traducere operatia poate fi reprezentat de o matrice omogenă, şi astfel o matrice poate fi folosit pentru a traduce omogene de vectori.
Dacă dorim să reprezinte poziţiile şi direcţiile în spaţiu (d)-dimensional, şi am dori să utilizaţi matrici de a traduce poziţii, apoi matrici omogene şi omogene de vectori pot fi folosite.
Cu toate acestea, omogene vectorilor trebuie să aibă (d+1) componente, matrici omogene şi trebuie să aibă ((d+1)*(d+1)) intrări.
Rezultatul final, o componentă a vectorului de omogene şi ultimul rând de o matrice omogenă nu corespund nici uneia dintre (d) coordonatele pe care dorim să o reprezinte în (d)-dimensional spaţiu, dar în loc să servească numai la spre a face it posibil la spre a lua matrici de a face unele operaţiuni de independent de a coordona (d) valori.
Utilizarea matrici omogene şi omogene de vectori face mai uşor de compus transformări şi de a aplica transformarile de vectori, dar şi alte operaţiuni devin mai inconvenient.
Representing poziţii sau direcţii în (d)-dimensional spaţiu omogen, folosind vectori de vectori înseamnă că va avea (d+1) componente, cu final componentă egală cu 1.
Calculul dot produsul a două omogene de vectori, folosind un produs regulate dot-ar combina implicit funcţia de punct de produse de prima (d) componente, cu o valoare de 1 irelevante din cauza finală a componentelor de 1.
Prin urmare, pentru a ajunge la punctul de produs în cauză (d) componente ale vectorului de omogene (cu (d+1) componente), punct de produse omogene de vectori ar trebui să se facă utilizând un "vector de omogene punct de produse" pur şi simplu că funcţia skips final componentă a vectorilor.
Alternativ, punct de produs ar putea fi calculate cu regularitate punct de produs cu funcţia, în sensul că o sumă de 1 ar trebui să se scade din rezultatul înainte de interpretare a rezultat ca legate de (d) relevante coordonate.
Codul de mai jos arată cum o matrice poate fi format care poate fi apoi folosit pentru a traduce un vector.
public class MatrixF64
{
// . . .
public static MatrixF64 FormHomogeneousTranslationMatrix( VectorF64 v )
{
// This method assumes the specified vector is a homogeneous
// vector, where the final component is not relevant to the
// translation. Thus, we form a square homogeneous matrix
// with a total number of rows equal to the number of
// components of the specified vector. We form an identity
// matrix of the required size, and copy all but the final
// component of the specified vector to the final column of
// the matrix.
if (null == v)
{
return (MatrixF64.Zero( 1, 1 )); // Vector is not specified.
}
if (v.Dimensions() <= 0)
{
return (MatrixF64.Zero( 1, 1 )); // Vector is empty.
}
MatrixF64 result =
MatrixF64.Identity( v.Dimensions() );
// Add all but the final component of the specified vector
// to the final column of the result matrix. The final
// entry of the final column of the matrix will remain one (1).
for (int i = 0; i < (result.totalRows - 1); i++)
{
result[i, (result.totalColumns - 1)] = v[i];
}
return (result);
}
// . . .
public static void Test()
{
// . . .
// Example of using a homogeneous matrix to
// translate a homogeneous vector:
// Forming a homogeneous matrix that represents
// a translation:
VectorF64 homogeneousTranslationVector4x1 =
new VectorF64( 10.0, 20.0, 30.0, 1.0 );
MatrixF64 translation4x4 =
MatrixF64.FormHomogeneousTranslationMatrix
( homogeneousTranslationVector4x1 );
translation4x4.WriteLine();
// [ 1, 0, 0, 10 ]
// [ 0, 1, 0, 20 ]
// [ 0, 0, 1, 30 ]
// [ 0, 0, 0, 1 ]
// Using the homogeneous matrix to translate
// a homogeneous vector:
VectorF64 homogeneousVector =
new VectorF64( 5.0, 6.0, 7.0, 1.0 );
VectorF64 translatedVector =
translation4x4 * homogeneousVector;
translatedVector.WriteLine();
// ( 15, 26, 37, 1 )
// The relevant vector components (5,6,7)
// were translated by (10,20,30) by the
// homogeneous translation matrix.
// . . .
}
}
3.21 (d)-dimensional spaţiu: rotaţii
"de rotaţie în spaţiu (d)-dimensionala:" (R) specificat o operaţie de către o bidimensional (S) planul de rotaţie, un unghi de rotaţie (t), de rotaţie şi un punct central (c), care converteşte orice punct (p) la un punct (q) corespunzătoare, astfel că:
(1) Length( q - c ) = Length( p - c ); distanţa dintre noi şi punctul de punct central de rotaţie este aceeaşi cu distanţa dintre punctul iniţial şi punctul de centrul de rotaţie.
Astfel, noul punct de pe suprafaţa de un anumit (d)-dimensional sferă;
(2) Sa (P) fi un plan care este paralel cu planul de rotaţie (S) şi conţine punctul (p).
Punctul (q) va fi în avionul (P).
Astfel, noul punct este pe un anumit plan bidimensional în spaţiul (d)-dimensională;
(3) Dot( (q - c), (p - c) ) = Length( q - c ) * Length( p - c ) * Cos( t ); produs din punct de vectorul de la centrul de rotaţie la punctul iniţial şi punctul de vector de la centrul de rotaţie punct la punct de noi, este egal cu produsul dintre cele două lungimi de vectori şi cosine de rotaţie unghi.
Astfel, noul punct de pe suprafaţa de un anumit (d)-dimensional cone;
Constrângeri de (1) limitele noul punct de la suprafata de un anumit domeniu (d)-dimensional.
Stare (2) suplimentare constrânge pe noi la un anumit punct bidimensional plane în spaţiu (d)-dimensional.
De aceea, noi trebuie să fie punct de pe un anumit cerc bidimensional în spaţiul (d)-dimensional.
Stare (3) suplimentare constrânge pe noi la un anumit punct (d) cone-dimensional.
Prin urmare, noul punct trebuie să fie la un anumit punct (de Cos( t ) = 1 sau Cos( t ) = (-1)) sau la una din cele două puncte posibile (pentru (-1) < Cos( t ) < 1).
Un final constrângere este necesară pentru a alege una dintre cele două puncte care, posibil rezultat de la constrângerile de mai sus.
Să T(s) reprezintă o operaţiune "de traducere," în cazul în care (s) este un vector arbitrar, astfel că T(s) * (r) = (r + s), şi Inverse(T(s)) * (r) = (r - s).
O rotatie R(S,t,c), pentru S plane, unghiul t, şi punct central c, pot fi exprimate în termeni de traducere şi a unei operaţiuni de rotaţie cu un punct central de la origine: R(S,t,c) = T(c) * R(S,t,0) * T(-c) = T(c) * R(S,t,0) * Inverse(T(c)).
Să R(S,t) reprezintă o operaţie de rotaţie cu un punct central de rotatie la originea (d)-dimensional space.
Astfel, orice arbitrar de rotaţie poate fi exprimată ca: R(S,t,c) = T(c) * R(S,t) * T(-c) = T(c) * R(S,t) * Inverse(T(c)).
Să ne considerăm bidimensional de rotaţie (S) care includ planuri de origine a (d)-dimensionale şi spaţiul sunt paralel la două axe de coordonate.
Să (a) şi (b) fie întregi care reprezintă axele de coordonate în spaţiu (d)-dimensional, astfel că 0 <= a <= (d-1), şi 0 <= b <= (d-1), şi (a != b).
Noi definim R(a,b,t) a fi o rotaţie despre originea (d)-dimensional spaţiu, cu un avion de rotaţie paralele cu axele de coordonate indicate de a coordona Indicii (a) şi (b), şi cu un unghi (t), astfel că R(a,b,(pi/2)) ar converti o valoare pozitivă coordoneze pentru a coordona axa indicată de (a) de a coordona o valoare pozitivă pentru a coordona axa indicată de (b); exemplu: R(a,b,(pi/2)) * e(a) = e(b), în cazul în care e(k) indică o unitate de baza vectorului corespunzătoare pentru a coordona (k) axă.
Această definiţie oferă final de dezambiguizare necesare pentru a determina pe deplin noul punct de produse de o rotaţie.
O rotatie care se mută-a lungul axei a faţă de axa b în conformitate cu un unghi t este echivalentă cu o rotaţie, care se mută-a lungul axei b faţă de axa a de către un unghi (-t); adică, R(a,b,t) = R(b,a,(-t)).
Dacă (d >= 2), apoi (d)-dimensionala a spatiu ((d*(d-1))/2) perechi de axe de coordonate.
Astfel, pentru (d >= 2), (d)-dimensionala a spatiu-((d*(d-1))/2) distincte axa de rotaţie avioane aliniate.
(Exemple: d=2: () x-y plane; d=3: (x-y avionul, y-z avionul, z-x avionul); d=4: (x-y avionul, y-z avionul, z-x avionul, x-w avionul, z-w avionul, w-y avionul).)
Codul de mai jos arată cum o matrice poate fi utilizată pentru a reprezenta R(a,b,t) de rotaţie.
Un exemplu este arătat că foloseşte o matrice de rotaţie într-o manieră omogenă la spre rotate un vector de omogene la un nou vector de omogene.
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 Propuneri pentru definirea "counterclockwise"
În spaţiul bidimensional, la rândul său, un sfert de rotaţie, care poate converti vectorului e(0) la e(1) este considerat counterclockwise.
De rotaţie este: R(0,1,(pi/2)) * e(0) = e(1).
Prin urmare, este R(0,1,t) rotaţie counterclockwise.
În tri-dimensională de spaţiu, şi-a rândul său, un sfert de rotaţie, care poate converti vectorului e(0) la e(1) este considerat counterclockwise.
La rândul său, un sfert de rotaţie, care poate converti vectorului e(1) la e(2) este, de asemenea, considerate counterclockwise.
La rândul său, un sfert de rotaţie, care poate converti vectorului e(2) la e(0) este, de asemenea, considerate counterclockwise.
Prin urmare, rotaţii { R(0,1,t), R(1,2,t), R(2,0,t) } sunt counterclockwise.
Logic formulă de mai jos se încadrează în definiţia de counterclockwise pentru bidimensional şi tri-dimensională de spaţiu, şi extrapolates modelul general pentru cazul de rotaţie R(a,b,t):
bool counterClockwise =
(
( (a < b) && (((a + b) % 2) == 1) )
|| ( (a > b) && (((a + b) % 2) == 0) )
);
Această regulă ar defini ca "counterclockwise" rotaţii de mai jos:
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 spaţiu: în sens orar rotaţii
Folosind definiţia propusă de "counterclockwise" în secţiunea anterioară, de rotaţii de mai jos sunt "sens orar:"
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),
...
De exemplu, în patru-dimensionale spaţiu, cele şase matrici de rotaţie de mai jos reprezintă rotaţii în sens orar:
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 spaţiu: numărul total de axă-aliniat orientări
Numărul de distincte axa-aliniat orientări în (d)-dimensional spaţiu este:
distinctOrientations = Factorial( d ) * IntPow( 2, (d-1) );
Un argument în sprijinul acestei formule apare în paragrafele de mai jos.
Dovezi experimentale verificarea formulă este, de asemenea, prezentate mai jos.
Să { g(0), g(1), ..., g(d-1) } fi un set de (d) reciproc perpendiculare-unitate de vectori ataşate la un obiect în spaţiu (d)-dimensional.
Aceste indica vectorii unitate de orientare a obiectului în raport cu axele de coordonate ale (d)-dimensional space.
Să "implicit" de "orientare" cu privire la obiectul g(i) = e(i) fie astfel încât, pentru toate 0 <= i <= (d-1); fiecare orientare a vectorului de obiect este asociat cu o unitate distinctă de baza vectorului de un sistem de coordonate în spaţiu (d)-dimensional.
Prima orientare a vectorului de obiect pot fi paralel cu oricare dintre (d) axe de (d)-dimensional spaţiu, şi poate fi în oricare din doua directii relativ la care axa.
Astfel, prima de orientare a vectorului (2 * d) opţiuni.
Cea de-a doua orientare a vectorului de obiect pot fi paralel cu oricare dintre (d-1) restul de axe, cu doua directii posibile în raport cu care axa.
Astfel, cea de-a doua orientare a vectorului (2 * (d-1)) opţiuni.
În acelaşi mod, în următorii orientare a vectorului (2 * (d-2)) opţiuni.