Algebra lineare
Colin Fahey
1. Software
LinearAlgebra.zip
Algebra lineare codice sorgente (C#)
19910 bytes
MD5: 11d8c8035cac30ba543e5e0b72ee9767
2. Introduzione
In questo articolo viene descritto vettori e matrici in (d)-dimensionale spazio.
3. (d)-dimensionale spazio: attributi
3.1 Array
“array:„ Una raccolta di tali variabili che ogni variabile ha un nome univoco, e tale che il nome può essere assegnato un ordine.
Valori interi possono essere utilizzati come nomi per le variabili in un array.
Ad esempio, se un array contiene (d) variabili, quindi la { 0, 1, 2, ..., (d-1) } valori interi possono essere i nomi assegnati alle variabili nella matrice.
3.2 (d)-dimensionale Vector
“(d)-dimensionale vettore:„ una serie di variabili (d).
“vettore componente:„ una variabile all'interno di un vettore.
3.3 (d)-dimensionale spazio vettoriale
“uno spazio-dimensionale:„ Il set completo di valori che possono essere archiviati da una variabile.
“(d)-dimensionale spazio:„ La serie completa di combinazioni di valori che possono essere archiviati da una serie di variabili (d).
La definizione formale di uno “spazio vettoriale:„
Lasciate (T) essere un tipo di base (esempi: numero reale, intero, numero complesso, numero razionale, ecc.)
Qualsiasi variabile di un tipo di base è chiamato un “scalare.„
Un “(T) di tipo (d)-dimensionale spazio vettoriale„ è l'insieme (S) di (d)-dimensionale vettori che hanno due operazioni, oltre (+) vettore, e scalare (*) moltiplicazione, che soddisfano le condizioni seguenti:
(1) Se (v) e (w) sono due vettori in (S), quindi (v + w) è anche un vettore in (S);
(2) Se (u), (v), e (w) sono di tre vettori in (S), quindi (u + v) + w = u + (v + w);
[additivo commutativity]
(3) Se (v) e (w) sono due vettori in (S), quindi (v + w) = (w + v);
[additivo Associatività]
(4) Vi è un “vettore zero,„ (0), in (S), tale che per ogni vettore (v) in (S), (v + (0)) = v;
[additivo identità]
(5) Se è (c) qualsiasi tipo di scalare (T), e (v) è un qualsiasi vettore in (S), quindi il prodotto (c * v) è un vettore in (S);
(6) Se (a), (b), e (c) sono scalari di qualsiasi tipo (T), e (v) e (w) sono vettori in qualsiasi (S), quindi (a + b) * v = a*v + b*v, e c*(v + w) = c*v + c*w;
[moltiplicativo Distributività]
(7) Se (a) e (b) sono scalari di qualsiasi tipo (T), e (v) è un qualsiasi vettore in (S), quindi (a*b)*v = a*(b*v);
(8) Se “1„ è un tipo di scalare (T) tale che (1*1)=1, e (v) è un qualsiasi vettore in (S), quindi (1*v) = v;
(9) Per ogni vettore (v) in (S), il vettore (-1)*v = -v soddisfa v + (-v) = (0);
[additivo inverso]
3.4 (d)-dimensionale Vector codice
Il codice qui sotto mostra come una (d) vettore-dimensionale, con 64-bit in virgola mobile componenti, può essere implementato.
Un array è sufficiente a rappresentare un vettore.
Il codice qui sotto ha una matrice contenuta in una classe, solo per comodità.
Il codice non è destinato ad essere efficace.
Una struttura (ad esempio: “struct„, un tipo di valore) che rappresentano vettori di un numero fisso di dimensioni (esempi: 3 o 4) rischia di essere molto più efficienti di quelli generali (d)-vettore classe dimensionale mostrato qui.
Sebbene il codice qui sotto definisce un vettore con virgola mobile componenti, questo documento fa anche uso di vettori con componenti intero.
Il codice qui sotto può essere facilmente modificato per attuare con vettori intero componenti.
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)-dimensionale zero vettore:„ un vettore con tutti i componenti (d) pari a zero.
public class VectorF64
{
// . . .
public static VectorF64 Zero( int dimensions )
{
VectorF64 zero = new VectorF64();
zero.components = new double[dimensions];
for (int i = 0; i < dimensions; i++)
{
zero[i] = 0.0;
}
return (zero);
}
// . . .
public static void Test()
{
// . . .
// An 8-dimensional vector with all 8 64-bit floating-point
// components set to zero:
VectorF64 z = VectorF64.Zero( 8 );
z.WriteLine(); // ( 0, 0, 0, 0, 0, 0, 0, 0 )
// . . .
}
}
3.5 (d)-dimensionale spazio: punti
“(d)-dimensionale spazio punto:„ una serie di (d) variabili con valori specifici “(valori delle coordinate).„
“(d)-dimensionale spazio origine:„ una serie di variabili (d) con tutti i valori pari a zero.
3.6 (d)-dimensionale vettori: non relativa e relativa
Un (d)-dimensionale vettore che “non„ è “parente„ è un (d)-dimensionale vettore che direttamente rappresenta una situazione o configurazione.
Un (d) vettore-dimensionale che è “relativo„ è un (d)-dimensionale vettore che rappresenta le modifiche a un insieme di componenti.
Un parente vettore può rappresentare la differenza tra i due non relativi vettori.
Dato un vettore relativa, determinando un non-stato e relativa configurazione utilizzando relativa vettore che richiede che unisce la relativa vettore con una percentuale di non rispetto vettore.
Non relativi vettori e relativi vettori sono entrambi vettori.
Se un particolare vettore non è parente o un parente deve essere specificato quando il vettore è definito.
(d) se un vettore-dimensionale viene interpretato come non relativa, quindi la (d)-dimensionale vettore può rappresentare un punto in (d)-dimensionale spazio.
3.7 (d)-dimensionale vettore addizione, sottrazione, e scaling
Vettore di addizione, sottrazione, e scaling:
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)-dimensionale base vettori
La definizione formale di una “base„ di uno spazio vettoriale:
Lasciate (T) essere un tipo di base (esempi: numero reale, intero, numero complesso, numero razionale, ecc.)
Qualsiasi variabile di un tipo di base è chiamato un “scalare.„
Lasciate (V) essere un “(T) di tipo (d)-dimensionale spazio vettoriale.„
Se il non-zero vettori { u1, u2, ..., ud } in (V) sono tali che ogni vettore (v) in (V) può essere scritto come “combinazione lineare„ di tali vettori, v = c1*u1 + c2*u2 + ... + cd*ud, dove sono { c1, c2, ..., cd } scalari di tipo (T), allora è (V) “spanning„ di vettori { u1, u2, ..., ud }.
Un gruppo di non-zero vettori che { u1, u2, ..., ud } “span„ spazio vettoriale (V) è chiamato a “base„ di (V).
Una semplice “base„ di una (d)-dimensionale spazio vettoriale è l'insieme di (d) distinti (d)-dimensionale vettori, ciascuno con un elemento pari a uno e tutti gli altri componenti pari a zero.
Base di tali vettori sono “Orthonormal,„ il che significa che esse sono reciprocamente-perpendicolare “(ortogonale)„ e che ogni vettore ha unità di lunghezza.
Ciascuno di tali vettori è un vettore unitario parallelo ad uno degli coordinerà (d) assi.
Esprimendo un arbitrario vettore come “combinazione lineare„ di questi vettori di base è diretta; componente di ogni arbitrario il vettore è moltiplicato per il corrispondente vettore di base, e questi prodotti sono aggiunti insieme per formare il vettore arbitrario.
Il codice qui sotto mostra come un vettore può essere espresso come “combinazione lineare„ di vettori di base.
Il codice qui sotto arbitrariamente Orthonormal definisce un insieme di vettori di base.
public class VectorF64
{
// . . .
public static VectorF64 BasisVector( int dimensions, int componentIndex )
{
if (dimensions < 0)
{
// Invalid number of dimensions specified.
return (new VectorF64());
}
VectorF64 basisVector = VectorF64.Zero( dimensions );
if ((componentIndex >= 0) && (componentIndex < dimensions))
{
basisVector[ componentIndex ] = 1.0;
}
return (basisVector);
}
// . . .
public static void Test()
{
// . . .
// The 4 basis vectors of 4-dimensional space,
// each with 4 64-bit floating-point components:
VectorF64 b0 = VectorF64.BasisVector( 4, 0 );
b0.WriteLine(); // ( 1, 0, 0, 0 )
VectorF64 b1 = VectorF64.BasisVector( 4, 1 );
b1.WriteLine(); // ( 0, 1, 0, 0 )
VectorF64 b2 = VectorF64.BasisVector( 4, 2 );
b2.WriteLine(); // ( 0, 0, 1, 0 )
VectorF64 b3 = VectorF64.BasisVector( 4, 3 );
b3.WriteLine(); // ( 0, 0, 0, 1 )
// . . .
}
}
Qualsiasi vettore in (d)-dimensionale spazio può essere espresso come somma dei prodotti di numeri e vettori di base:
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)-dimensionale spazio: la distanza tra i punti
Lasciate (P) essere un (d)-dimensionale vettore che rappresenta un punto in (d)-dimensionale spazio.
Lasciate (Q) essere un (d)-dimensionale vettore che rappresenta un punto in (d)-dimensionale spazio.
Lasciate (R) essere un (d)-dimensionale vettore che rappresenta la variazione del (d) coordinate per andare dal punto (P) a punto (Q); R = (Q - P).
In uno spazio-dimensionale, P = ( p0 ), Q = ( q0 ), e R = (Q - P) = ( q0 - p0 ).
La distanza tra i due punti è: Abs( q0 - p0 ).
In due dimensioni spazio, P = ( p0, p1 ), Q = ( q0, q1 ), e R = (Q - P) = ( q0-p0, q1-p1 ).
Interpretare le due perpendicolari spostamenti come la perpendicolare lati di un triangolo-destra, la distanza fra i punti corrisponde alla lunghezza del hypotenuse di quel triangolo.
Pythagorean la formula, (a*a) + (b*b) = (c*c), dove (a) e (b) sono le lunghezze dei lati perpendicolari di un diritto-triangolo, e (c) è la lunghezza della hypotenuse (skew lato), può essere usato per determinare la distanza tra i due punti: Sqrt( Sq(q0-p0) + Sq(q1-p1) ).
In spazio tridimensionale, P = ( p0, p1, p2 ), Q = ( q0, q1, q2 ), e R = (Q - P) = ( q0-p0, q1-p1, q2-p2 ).
Interpretando lo spostamento ( q0-p0, q1-p1, 0 ) come la perpendicolare lati di un triangolo-destra, e utilizzando la formula Pythagorean, ha pronunciato la distanza tra il punto (P) e il punto ( q0, q1, p2 ): d01 = Sqrt( Sq(q0-p0) + Sq(q1-p1) ).
Lo spostamento ( q0-p0, q1-p1, 0 ) è perpendicolare allo sfollamento ( 0, 0, q2-p2 ), e un altro destro del triangolo può essere formato, e la Pythagorean formula può essere usato di nuovo.
Così, la distanza dal punto (P) a punto (Q) è dato da: 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) ).
Il metodo di estendere la formula di distanza da due dimensioni spazio a spazio tridimensionale può essere applicato più volte per poi determinare la distanza formula per (d)-dimensionale spazio: Sqrt( (Sq(q0-p0) + Sq(q1-p1)) + Sq(q2-p2) + ... + Sq(qd-pd) ).
Il codice qui sotto definisce una funzione denominata “Lunghezza„ che calcola la lunghezza di una (d) vettore-dimensionale.
Quando un vettore rappresenta lo spostamento tra due punti in (d)-dimensionale spazio, la lunghezza del vettore rappresenta la distanza tra questi due punti.
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)-dimensionale vettori: dot prodotto
Il “punto prodotto„ converte due (d)-dimensionale vettori a un numero.
Il codice qui sotto calcola un punto prodotto di due vettori:
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);
}
// . . .
}
Per ogni vettore, (A), Length(A) = Sqrt(Dot(A,A)).
3.11 (d)-dimensionale vettori: definizione di “parallelo„
(A) vettori e (B) sono “paralleli„ se tutte le seguenti condizioni:
(1) A.Length() > 0;
(2) B.Length() > 0;
(3) Abs(Dot(A,B)) = A.Length()*B.Length().
Il codice qui sotto determina se una coppia di vettori sono paralleli (eventualmente anti-allineati).
Numeri floating-point possono accumulare errori nella parte frazionaria a causa della limitata precisione, e, quindi, di codice dovrebbe includere non-zero tolleranze, se si confrontano i numeri floating-point.
Il codice di esempio include i valori di tolleranza, ma l'esempio di tolleranza valori potrebbero non essere appropriato per alcuni compiti.
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)-dimensionale vettori: definizione di “perpendicolare„
(A) vettori e (B) sono “perpendicolari„ se tutte le seguenti condizioni:
(1) A.Length() > 0;
(2) B.Length() > 0;
(3) Abs(Dot(A,B)) = 0.
Il codice qui sotto determina se una coppia di vettori sono perpendicolari. Numeri floating-point possono accumulare errori nella parte frazionaria a causa della limitata precisione, e, quindi, di codice dovrebbe includere non-zero tolleranze, se si confrontano i numeri floating-point.
Il codice di esempio include i valori di tolleranza, ma l'esempio di tolleranza valori potrebbero non essere appropriato per alcuni compiti.
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 Matrici
“matrice:„ Una raccolta di tali variabili che ogni variabile è una combinazione unica di una “riga„ nome e un nome di “colonna.„
“ingresso:„ una variabile all'interno di una matrice.
Valori interi possono essere utilizzati come nomi di “riga„ e “colonna di„ nomi per le variabili in una matrice.
Ad esempio, se una matrice ha (totalRows) righe e colonne (totalColumns), allora la valori interi { 0, 1, ..., (totalRows-1) } può essere il nome assegnato al file, e la { 0, 1, ..., (totalColumns-1) } valori interi possono essere i nomi assegnati alle colonne.
Pertanto, una variabili in un modello che può essere specificato indicando un paio di numeri interi, ( row, column ), con l'indicazione della combinazione di riga e colonna indici corrispondenti alle specifiche variabile.
Le dimensioni di una matrice è specificato come “(totalRows) * (totalColumns)„ (o “(totalRows) di (totalColumns)).„
Questo ordine di dimensioni è lo stesso che l'ordine di dimensioni utilizzato per specificare le voci nella matrice “(( row, column )).„
[Questa convenzione è un po 'un peccato, perché per molti bidimensionale usi (esempi: immagini, grafici, ecc) la solita convenzione è quello di specificare le dimensioni come “width * height„ e coordina come “( horizontal, vertical )„ (o “( x, y )).„
Questo è l'opposto dell 'ordine di dimensioni e le coordinate usato per descrivere le matrici e le loro voci. ]
Una matrice con (totalRows) pari al (totalColumns) è chiamato “quadrato,„ altrimenti la matrice è chiamato “rettangolare.„
Una matrice può essere considerata come contenente una serie di “vettori riga,„ dove le variabili in ciascuna riga vengono interpretati come appartenenti a un vettore.
Una matrice può anche essere considerata come contenente una serie di “vettori colonna,„ in cui le variabili in ogni colonna vengono interpretati come appartenenti a un vettore.
Matrici possono rappresentare una grande varietà di rapporti matematici.
Il significato di una matrice, e le operazioni che potrebbero essere appropriati per la trasformazione di voci di una matrice, dipende dal contesto.
Tuttavia, vi sono regole di base della matrice aritmetica che sono rilevanti per molti contesti, e queste regole di base verranno definite in una successiva sezione.
Un array e un valore (totalColumns) sono sufficienti a rappresentare una matrice.
L'array può avere (totalRows * totalColumns) variabili, e la data di entrata in ( row, column ) può corrispondere a matrice variabile a Indice ((totalColumns * row) + column).
Il codice qui sotto definisce una matrice, con 64-bit in virgola mobile voci.
Un array e un valore (totalColumns) sono sufficienti a rappresentare una matrice.
Il codice qui sotto è un array e un valore (totalColumns) contenute in una classe, solo per comodità.
Il codice qui sotto non è destinato ad essere efficace.
Una struttura (ad esempio: “struct„, un tipo di valore), che rappresentano particolare matrici di dimensioni fisse (esempi: 2*2, 3*3, o 4*4) rischia di essere molto più efficienti di quelli generali di classe (totalRows * totalColumns) mostrato qui.
Sebbene il codice qui sotto definisce una matrice con virgola mobile voci, questo articolo fa anche uso di matrici con voci intero.
Il codice qui sotto può essere facilmente modificata in modo da attuare con matrici intero voci.
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:„ Una matrice con tutte le voci pari a zero.
public class MatrixF64
{
// . . .
public static MatrixF64 Zero( int rows, int columns )
{
MatrixF64 zero = new MatrixF64( rows, columns );
// The constructor above sets all entries to zero.
return (zero);
}
// . . .
public static void Test()
{
// . . .
// An 8x2 matrix (8 rows x 2 columns) with all 16
// 64-bit floating-point entries set to zero:
MatrixF64 z = MatrixF64.Zero( 8, 2 );
z.WriteLine();
// [ 0, 0 ] // row 0
// [ 0, 0 ] // row 1
// [ 0, 0 ] // row 2
// [ 0, 0 ] // row 3
// [ 0, 0 ] // row 4
// [ 0, 0 ] // row 5
// [ 0, 0 ] // row 6
// [ 0, 0 ] // row 7
// . . .
}
}
3.14 Matrice addizione, sottrazione, moltiplicazione e
Matrice addizione, sottrazione, moltiplicazione e.
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 ]
// . . .
}
}
“matrice identità:„ una matrice quadrata (in totale righe colonne totale è pari a) a voci pari a (1) sulla diagonale (voci con un indice di riga uguale al suo indice di colonna), e con tutte le altre voci pari a (0).
Se una matrice quadrata, (M), è moltiplicata da una “matrice identità,„ (I), dello stesso numero di righe (e colonne), il prodotto è pari a (M).
Moltiplicando (I) di (M) produce anche un prodotto pari al (M).
In tal modo, “l'identità matrice„ è analogo al numero “1„ per la moltiplicazione dei numeri (scalari).
Il codice qui sotto crea una matrice identità con un determinato numero di righe.
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 Matrice: LU di factoring; sostituzione di back -
“LU matrice di factoring:„ una procedura che converte una matrice quadrata, (M), a due nuove matrici quadrate, (L) e (U), avendo la stessa dimensione di (M), in modo tale che (L) * (U) = (M), e tale che (U) matrice è “triangolare superiore„ (tutte le voci al di sotto della diagonale sono zero), e matrice (L) è “inferiore triangolare„ (tutte le voci di sopra della diagonale sono pari a zero).
LU matrice di factoring può essere utilizzato come parte di una più ampia procedura per risolvere un sistema di equazioni, o per trovare l'inversa di una matrice, o per trovare il determinante di una matrice.
LU matrice di factoring è un'alternativa al “Gauss-Jordan la procedura di eliminazione.
Gauss-Jordan eliminazione richiede un sistema di equazioni (A)*(x)=(b), LU di factoring che richiede solo una matrice (A).
Inoltre, dopo la determinazione del LU di factoring di una matrice (A), è molto facile da determinare (x) avere qualsiasi (b).
La procedura per risolvere il vettore (x) in (A)*(x)=(b), dato (b) e la LU factoring (L)*(U)=(A), comporta una procedura che verrà chiamato “sostituzione di back-qui.„
Il codice qui sotto include una procedura per il calcolo LU fattori di qualsiasi piazza, non singolare matrice.
Il codice qui sotto include una procedura per fare LU di back-sostituzione.
Attenzione: il computer calcola codice qui sotto il L e U fattori di una riga-versione permeata di una determinata matrice.
Inoltre, i L e U matrice risultati sono combinati in un unico output matrice.
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 Matrice: determinante
“determinante:„ Una funzione che converte qualsiasi quadrato (n * n) matrice (A) a un numero, det(A), tale che:
(1) Se una matrice (B) risultati da scambiare due righe o due colonne, di una matrice (A), quindi det(B) = (-(det(A)));
(2) Se una matrice (B) risultati da moltiplicando ogni riga, colonna o qualsiasi, di una matrice (A), da un numero (c), quindi det(B) = c * det(A);
(3) Se una matrice (B) risultati di aggiungere un multiplo di una riga per riga un altro, o di aggiungere un multiplo di una colonna a un'altra colonna, poi det(B) = det(A).
(4) Se (A) è un superiore-matrice triangolare (A[i,j]=0 per tutti i (i>j)) o un basso matrice triangolare (A[i,j]=0 per tutti i (i<j)), poi det(A) = (A[0,0] * A[1,1] * ... * A[n-1,n-1]);
(Per esempio, il determinante di una identità è una matrice; det(I) = 1.)
Norme (1), (2), e (3), può essere utilizzato in un processo chiamato “Gaussian eliminazione, per convertire qualsiasi matrice quadrata a una matrice triangolare.
(4) regola può essere usato per calcolare il determinante di una matrice triangolare.
La definizione formale di un “determinante:„
Lasciate M[n](K) indicare l'insieme di tutti i (n * n) matrici oltre il campo (K).
“Determinante„ è l'unica funzione (F) tale che F:M[n](K) --> K con due attributi:
(1) F è alternata multilinear per quanto riguarda le colonne (o righe) e,
(2) F( I ) = 1.
Un “multilinear„ funzione (D) di un (n * n) matrice (A) può essere scritto come: 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]] ) ), se la somma è ripreso tutti i (n^n) combinazioni di 0 <= k[i] <= (n-1) e, se del e[j] riga (j) rappresenta l'identità di matrice.
Tale funzione può essere caratterizzato da valori di D( e[k[0]], e[k[1]], ..., e[k[n-1]] ).
“Un'alternarsi multilinear„ funzione è una funzione multilinear, (D), in modo tale che D( ..., e[i], ..., e[i], ... ) = 0, e D( ..., e[i], e[j], ... ) = (-( D( ..., e[j], e[i], ... ) )).
Così, a scambiare le colonne (o righe) cambia il segno della funzione.
Inoltre, il termine D( e[k[0]], ..., e[k[n-1]] ) è pari a zero se la combinazione di valori k[i] non coinvolgere tutti i valori univoci, se l'insieme di valori non è una permutazione dei numeri {0,1,...,(n-1)}.
D( I ) = 1 (D( e[0], e[1], ..., e[n-1] ) = 1) locazione, in combinazione con l'alternarsi multilinear attributi di cui sopra, significa che tutti i valori dei alternata multilinear funzione D( e[k[0]], e[k[1]], ..., e[k[n-1]] ) può essere determinata.
La funzione sarà diverso da zero solo se l'insieme di valori k[i] è una permutazione di {0,1,...,(n-1)}, e avrà una magnitudine di uno se non è non-zero, e il segno sarà il numero di operazioni di swap necessario per convertire la permutazione a permutazione identità.
Il valore di un determinante di un (n * n) matrice (A) può essere calcolata sommando i prodotti sotto forma di (sgn(p) * A[0,p[0]] * A[1,p[1]] * ... * A[(n-1),p[n-1]]) per ogni permutazione (p) dei numeri {0,1,2,...,(n-1)}.
Il termine denota (sgn(p)) la “firma„ (o di swap count) di permutazione (p), dove (sgn(p)) è pari a (+1) se (p) è “ancora permutazione,„ ed è pari a (-1) se (p) è una “permutazione dispari.„
Perché ci sono (n!) permutazioni dei numeri {0,1,2,...,(n-1)}, questa formula ha (n!) summands.
Altre procedure per il calcolo il determinante (esempi: quelle che comportano Gaussian eliminazione, o LU factoring, o di espansione da parte dei minori, ecc) implicitamente calcolare una gerarchia di intermedi quantitativi che consentono di tali procedure per calcolare il determinante in un ordinamento di (n^3) passi.
Se il determinante di una matrice è pari a zero, la matrice non ha un inverso.
Se il determinante di una matrice è non-zero, la matrice è un inverso.
Se una matrice rappresenta i coefficienti di un sistema lineare di equazioni, e il determinante di tale matrice è pari a zero, allora il sistema di equazioni non dispone di una soluzione unica.
Se una matrice rappresenta i coefficienti di un sistema lineare di equazioni, e il determinante di tale matrice è non-zero, allora il sistema di equazioni ha una soluzione unica.
Il codice qui sotto comprende le procedure per il calcolo il determinante di qualsiasi matrice quadrata.
Determinanti per 1*1, 2*2, 3*3, e 4*4 matrici sono calcolate esplicito di formule.
Esplicito utilizzando formule per calcolare la determinanti per le piccole matrici è probabile per calcolare risultati più velocemente che utilizzando una procedura generale per calcolare quelli determinanti.
La procedura per il caso generale utilizza LU factoring, ma ci sono molti altri metodi di calcolo determinanti.
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)-dimensionale vettori: cross prodotto
“croce prodotto:„ Una funzione associata a un (d) spazio-dimensionale che converte un elenco ordinato di (d-1) (d)-dimensionale vettori, { A(0), A(1), ..., A(d-2) }, ad un (d) vettore-dimensionale, in modo tale che la funzione ha gli attributi di seguito:
(1) Se A(i)=0 (zero vector), per qualsiasi 0 <= i <= (d-2), quindi Cross( A(0), A(1), ..., A(d-2) )=0 (zero vettoriale).
La croce prodotto di un elenco ordinato di vettori è pari a zero se ogni vettore in tale elenco è pari a zero;
(2) Se una qualsiasi coppia di vettori, A(i) e A(j), per qualsiasi (i!=j) con 0 <= i <= (d-2) e 0 <= j <= (d-2), sono paralleli (Abs(Dot(A(i),A(j))) = Length(A(i)) * Length(A(j))), quindi Cross( A(0), A(1), ..., A(d-2) )=0 (vettore zero).
La croce prodotto di un elenco ordinato di vettori è pari a zero se due vettori in tale elenco sono paralleli;
(3) Se ogni vettore A(i) è non-zero, e se non è A(i) parallelo con A(j) per tutti i (i!=j) con 0 <= i <= (d-2) e 0 <= j <= (d-2), quindi B = Cross( A(0), A(1), ..., A(d-2) ) è tale che Dot(A(i),B)=0 per tutti i 0 <= i <= (d-2).
Se la croce prodotto di un elenco ordinato di vettori è non-zero, allora la croce prodotto risultato è perpendicolare ad ogni vettore in ordine la lista dei vettori;
(4) Se (c) è un valore numerico, quindi 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) ).
Dato un elenco ordinato di vettori, la croce del prodotto elenco ordinato di vettori con una qualsiasi vettore moltiplicato per un determinato valore numerico è equivalente a croce moltiplicando il prodotto del elenco ordinato di vettori dal valore numerico;
(5) Dato due numeri (i) e (j), dove (i!=j) e 0 <= i <= (d-2) e 0 <= j <= (d-2), e un elenco ordinato di (d-1) (d)-dimensionale vettori, { A(0), A(1), ..., A(d-2) }, e un altro elenco ordinato di (d-1) (d)-dimensionale vettori, { B(0), B(1), ..., B(d-2) }, tali che per tutti i B(k) = A(k) con ((k != i) && (k != j)) 0 <= k <= (d-2), e tale che B(i) = A(j) e B(j) = A(i), quindi Cross( A(0), A(1), ..., A(d-2) ) = (-1) * Cross( B(0), B(1), ..., B(d-2) ).
Se la croce prodotto di un elenco ordinato di vettori è non-zero, allora che attraversano il prodotto ha di fronte il segno della croce prodotto dello stesso elenco ordinato di vettori con l'eccezione di due vettori esattamente essere scambiate (scambiati).
Data base vettori e(k), per 0 <= k <= (d-1), di (d) spazio-dimensionale, la croce prodotto di un elenco ordinato di vettori (d-1) può essere calcolata calcolando il determinante della matrice (d*d) qui di seguito:
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]
Il determinante della matrice che non sarà un numero, ma sarà invece un vettore, dove ci saranno i coefficienti numerici per il (d) base vettori.
Per due dimensioni spazio, e il vettore A = ( ax, ay ):
Cross( A )
= Determinant
(
e(x), e(y),
ax, ay
)
= ay * e(x)
- ax * e(y)
= ( ay, -ax )
Il risultato vettore è perpendicolare al vettore di ingresso unico.
x-y se l'aereo è vista in modo che x aumenta in una direzione destra e y aumenti in una direzione verso l'alto, quindi attraversare il prodotto risultato è un quarto di giro in “senso orario„ rispetto al vettore di ingresso.
Esempi di croce prodotti in due dimensioni spazio:
Cross( e(x) ) = (-1) * e(y)
Cross( e(y) ) = e(x)
Per spazio tridimensionale, i vettori e A = ( ax, ay, az ) e B = ( bx, by, bz ):
Cross( A, B )
= Determinant
(
e(x), e(y), e(z),
ax, ay, az,
bx, by, bz
)
= (ay*bz - az*by) * e(x)
+ (az*bx - ax*bz) * e(y)
+ (ax*by - ay*bx) * e(z)
= ( (ay*bz - az*by), (az*bx - ax*bz), (ax*by - ay*bx) )
Esempi di prodotti in croce spazio tridimensionale:
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)
Per quattro-dimensionale spazio, i vettori e A = ( ax, ay, az, aw ), B = ( bx, by, bz, bw ), e C = ( cx, cy, cz, cw ):
Cross( A, B, C )
= Determinant
(
e(x), e(y), e(z), e(w),
ax, ay, az, aw,
bx, by, bz, bw,
cx, cy, cz, cw
)
= ( ay*bz*cw - ay*bw*cz + az*bw*cy - az*by*cw + aw*by*cz - aw*bz*cy) * e(x)
+ (-ax*bz*cw + ax*bw*cz - az*bw*cx + az*bx*cw - aw*bx*cz + aw*bz*cx) * e(y)
+ ( ax*by*cw - ax*bw*cy + ay*bw*cx - ay*bx*cw + aw*bx*cy - aw*by*cx) * e(z)
+ (-ax*by*cz + ax*bz*cy - ay*bz*cx + ay*bx*cz - az*bx*cy + az*by*cx) * e(w)
= ( ( ay*bz*cw - ay*bw*cz + az*bw*cy - az*by*cw + aw*by*cz - aw*bz*cy),
(-ax*bz*cw + ax*bw*cz - az*bw*cx + az*bx*cw - aw*bx*cz + aw*bz*cx),
( ax*by*cw - ax*bw*cy + ay*bw*cx - ay*bx*cw + aw*bx*cy - aw*by*cx),
(-ax*by*cz + ax*bz*cy - ay*bz*cx + ay*bx*cz - az*bx*cy + az*by*cx) )
Esempi di croce prodotti in quattro dimensioni spazio:
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)
Il codice qui sotto calcola la cross-prodotto di due o più vettori.
La prova esempi dimostrano che il non-zero risultato vettori sono perpendicolari a tutti i vettori di input.
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 Matrice: inversa
“matrice inversa:„ Una matrice quadrata collegata ad un altro quadrato matrice delle stesse dimensioni tali che il prodotto di due matrici è la “matrice identità.„
Se una matrice è associato un inversa della matrice, non vi è solo una di queste inversa della matrice.
Se il determinante di una matrice è pari a zero, allora la matrice è chiamato “singolare,„ e la matrice non dispone di un associato inversa della matrice.
Il codice qui sotto calcola l'inverso di ogni piazza, non singolare matrice.
Inverses per 1*1, 2*2, 3*3, e 4*4 matrici sono calcolate esplicito di formule.
Esplicito utilizzando formule potrebbe essere più veloce di utilizzare una procedura generale.
L'esplicito formule sono anche affidabili.
La procedura per il caso generale utilizza LU factoring, ma ci sono molti altri metodi di calcolo inverses di 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 Moltiplicativo prodotto di una matrice e un vettore
Moltiplicativo il prodotto di una matrice (M) e un vettore (a) è un vettore (b); (M)*(a) = (b).
Matrice (M) “trasforma„ (a) vettore a vettore (b).
Il codice qui sotto moltiplicativo calcola il prodotto di una matrice e un vettore.
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 Omogeneo utilizzando matrici tradurre omogeneo vettori
“traduzione:„ Un'operazione T(s), specificato da un particolare vettore (s), che converte qualsiasi vettore (r) corrispondente ad un vettore (r + s); T(s) * (r) = (r + s).
T(-s) è l'inverso del T(s); Inverse(T(s)) = T(-s); T(-s) * (r) = (r - s).
Una traduzione operazione può essere rappresentato da una matrice, in modo tale che moltiplicando la matrice di traduzione di un vettore ha l'effetto di traduzione che vettoriali.
I contributi di una traduzione per il funzionamento dei componenti il risultato vettore non dipendono, direttamente o indirettamente, sui componenti del vettore originale.
Pertanto, una matrice di traduzione deve essere strutturato in un modo che consente di contribuire in una traduzione per il risultato, senza essere colpiti, direttamente o indirettamente, da componenti del vettore originale.
Un “omogeneo vettore„ è un vettore con un ultimo elemento pari al 1.
Un “omogenea matrice„ è una matrice con una riga finale con tutte le voci pari a 0 fatta eccezione per un finale di entrata 1.
Moltiplicativo il prodotto di una matrice omogenea e un altro omogenea matrice è una matrice omogenea.
La moltiplicazione di un prodotto omogeneo e di una matrice omogenea vettore è un vettore omogeneo.
Un attributo interessante di questo prodotto è che le voci in ultima colonna della matrice omogenea, con l'eccezione della finale entrata, sono aggiunti direttamente ai componenti del vettore il risultato, senza essere colpiti, direttamente o indirettamente, dai componenti di il vettore originale.
Così, una traduzione operazione può essere rappresentata da una matrice omogenea, e tale matrice può essere utilizzato per tradurre omogeneo vettori.
Se vogliamo rappresentare le posizioni e le direzioni in (d)-dimensionale spazio, e vogliamo usare matrici tradurre posizioni, quindi omogenea matrici ed omogenea vettori possono essere utilizzati.
Tuttavia, la omogeneo vettori devono avere (d+1) componenti, e la omogeneo matrici devono avere ((d+1)*(d+1)) voci.
L'ultima componente di un vettore omogeneo e la finale di una riga omogenea matrice non corrispondono ad uno qualsiasi degli (d) coordinate che vogliamo rappresentare in (d)-dimensionale spazio, ma invece servono solo a rendere possibile per ottenere matrici di fare alcune operazioni indipendentemente dal (d) valori delle coordinate.
Utilizzando matrici omogenee ed omogenea vettori rende più facile composti trasformazioni e applicare trasformazioni a vettori, ma anche altre operazioni di diventare più scomode.
Che rappresentano posizioni o le direzioni in (d)-dimensionale spazio omogeneo utilizzando vettori significa che i vettori dovranno (d+1) componenti, con l'ultimo elemento pari al 1.
Il punto di calcolo prodotto di due vettori omogenei utilizzando un normale prodotto dot-funzione permetterebbe di far implicitamente il punto di prodotto il primo (d) componenti con un valore di irrilevante 1 finale a causa di componenti di 1.
Pertanto, per ottenere il punto del prodotto (d) le componenti del vettore omogenea (con (d+1) componenti), il punto di prodotto omogeneo vettori dovrebbe essere fatto utilizzando una speciale “omogeneo vettore dot prodotto„ funzione che semplicemente ignora la componente finale dei vettori.
In alternativa, il punto prodotto potrebbe essere calcolata con la regolare funzione di punto di prodotto con l'intesa che un importo di 1 dovrebbe essere sottratto dal risultato prima di interpretare il risultato in riguardanti il (d) pertinenti coordinate.
Il codice qui sotto mostra come una matrice può essere formato che può essere quindi utilizzato per tradurre un vettore.
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)-dimensionale spazio: rotazioni
“rotazione in (d)-dimensionale spazio:„ Un'operazione (R) specificato da una a due dimensioni (S) piano di rotazione, un angolo di rotazione (t), e di un centro di rotazione (c) punto, che converte qualsiasi punto (p) ad un punto corrispondente (q), tale che:
(1) Length( q - c ) = Length( p - c ); la distanza tra il nuovo punto di rotazione e il punto centrale è la stessa della distanza tra il punto di partenza e il punto centrale di rotazione.
Di conseguenza, il nuovo punto è sulla superficie di un particolare (d)-dimensionale sfera;
(2) Che (P) essere un aereo che è parallelo al piano di rotazione e (S) contiene (p) punto.
(q) punto sarà in aereo (P).
Di conseguenza, il nuovo punto è su un particolare bidimensionale in aereo (d)-dimensionale spazio;
(3) Dot( (q - c), (p - c) ) = Length( q - c ) * Length( p - c ) * Cos( t ); prodotto il punto del vettore dal centro di rotazione punto a punto l'originale e il vettore dal centro di rotazione a punto il nuovo punto è uguale al prodotto delle lunghezze di questi due vettori e il coseno di rotazione angolo.
Di conseguenza, il nuovo punto è sulla superficie di un particolare (d) cono-dimensionale;
(1) vincolo limita il nuovo punto sulla superficie di un particolare (d)-dimensionale sfera.
(2) ulteriore condizione vincolante il nuovo punto a un particolare bidimensionale in aereo (d)-dimensionale spazio.
Pertanto, il nuovo punto deve essere su un particolare bidimensionale cerchio in (d)-dimensionale spazio.
(3) ulteriore condizione vincolante il nuovo punto a un particolare (d) cono-dimensionale.
Pertanto, il nuovo punto deve essere a un punto particolare (per Cos( t ) = 1 o Cos( t ) = (-1)) o in uno dei due possibili punti (per (-1) < Cos( t ) < 1).
Un ultimo vincolo è richiesto di scegliere uno dei due possibili punti che derivano da vincoli di cui sopra.
Lasciate T(s) rappresenta una “traduzione„ operazione, se (s) è un vettore arbitrario, tale che T(s) * (r) = (r + s), e Inverse(T(s)) * (r) = (r - s).
Una rotazione R(S,t,c), per S aereo, angolo t, e punto centrale c, può essere espresso in termini di operazioni di traduzione e di una rotazione con un punto centrale all'origine: R(S,t,c) = T(c) * R(S,t,0) * T(-c) = T(c) * R(S,t,0) * Inverse(T(c)).
Lasciate R(S,t) rappresenta una operazione di rotazione con un punto di rotazione al centro l'origine dei (d)-dimensionale spazio.
Così, ogni arbitraria rotazione può essere espresso come: R(S,t,c) = T(c) * R(S,t) * T(-c) = T(c) * R(S,t) * Inverse(T(c)).
Si consideri bidimensionale rotazione (S) aerei che includono l'origine del (d)-dimensionale spazio e sono in parallelo a due assi di coordinate.
Lasciate (a) e (b) essere interi che rappresentano assi di coordinare in (d)-dimensionale spazio, tale che 0 <= a <= (d-1), e 0 <= b <= (d-1), e (a != b).
Definiamo R(a,b,t) ad essere una rotazione circa l'origine del (d) spazio-dimensionale, con una rotazione piano parallelo a coordinare assi indicato da coordinare indici (a) e (b), e con un angolo (t), in modo tale che R(a,b,(pi/2)) sarebbe convertire un valore positivo per coordinare la coordinare asse indicato da (a) ad un valore positivo per coordinare la coordinare asse indicato da (b); esempio: R(a,b,(pi/2)) * e(a) = e(b), dove e(k) indica una unità di base vettore corrispondente a coordinare asse (k).
Questo prevede la definizione finale disambigua necessaria onde accertare il nuovo punto prodotto da una rotazione.
Una rotazione che si muove lungo l'asse a verso l'asse b secondo un angolo t è equivalente a una rotazione che si muove lungo l'asse b verso a asse di un angolo (-t); vale a dire, R(a,b,t) = R(b,a,(-t)).
Se (d >= 2), quindi (d)-dimensionale spazio ha ((d*(d-1))/2) coordinare le coppie di assi.
Così, per (d >= 2), (d)-dimensionale spazio ha ((d*(d-1))/2) distinti asse di rotazione allineati aerei.
(Esempi: d=2: () x-y aereo; d=3: (x-y aereo, y-z aereo, z-x piano); d=4: (x-y aereo, y-z aereo, z-x aereo, x-w aereo, z-w piano, piano w-y).)
Il codice qui sotto mostra come una matrice può essere utilizzata per rappresentare la rotazione R(a,b,t).
Un esempio è dimostrato che utilizza una matrice di rotazione in un modo omogeneo per ruotare un vettore omogeneo a un nuovo vettore omogeneo.
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 Definizione proposta di “antiorario„
In due dimensioni spazio, il quarto di girare a rotazione che può convertire il vettore e(0) a e(1) è considerato antiorario.
La rotazione è: R(0,1,(pi/2)) * e(0) = e(1).
Pertanto, la rotazione R(0,1,t) è antiorario.
In spazio tridimensionale, i quarti di girare a rotazione che può convertire il vettore e(0) a e(1) è considerato antiorario.
I quarti di girare a rotazione che può convertire il vettore e(1) a e(2) viene considerata anche antiorario.
I quarti di girare a rotazione che può convertire il vettore e(2) a e(0) viene considerata anche antiorario.
Pertanto, la { R(0,1,t), R(1,2,t), R(2,0,t) } rotazioni sono antiorario.
L'operatore booleano formula qui di seguito si inserisce la definizione di antiorario per bidimensionale e tridimensionale spazio, e extrapolates il modello per il caso generale di rotazione R(a,b,t):
bool counterClockwise =
(
( (a < b) && (((a + b) % 2) == 1) )
|| ( (a > b) && (((a + b) % 2) == 0) )
);
Questa norma definisce le rotazioni “antiorario„ come qui di seguito:
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)-dimensionale spazio: rotazioni in senso orario
Utilizzando la proposta di definizione di “antiorario„ nella sezione precedente, le rotazioni “in senso orario„ sono qui di seguito:
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),
...
Per esempio, in quattro dimensioni spazio, il termine di sei matrici di rotazione al di sotto rappresentano rotazioni in senso orario:
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)-dimensionale spazio: numero totale di asse-allineati orientamenti
Il numero di distinte asse allineati orientamenti in (d)-dimensionale spazio è:
distinctOrientations = Factorial( d ) * IntPow( 2, (d-1) );
Un argomento sostenendo questa formula appare nei punti seguenti.
Prove sperimentali verificare la formula è anche presentato qui di seguito.
Lasciate { g(0), g(1), ..., g(d-1) } essere un insieme di (d) reciprocamente-perpendicolare unità vettori allegato a un oggetto in (d)-dimensionale spazio.
Queste unità di vettori di indicare l'orientamento dell'oggetto rispetto agli assi di coordinare le (d)-dimensionale spazio.
Lasciate il valore di “default orientamento„ dell'oggetto essere tale che g(i) = e(i), per tutti i 0 <= i <= (d-1); ogni vettore di orientamento l'oggetto è associato ad una distinta base vettoriale di un sistema di coordinate in (d)-dimensionale spazio.
Il primo vettore di orientamento l'oggetto può essere parallelo con uno qualsiasi degli assi (d) dei (d)-dimensionale spazio, e può essere in una delle due direzioni rispetto a quello asse.
Così, il primo vettore di orientamento ha (2 * d) opzioni.
Il secondo vettore di orientamento l'oggetto può essere parallelo ad uno qualsiasi degli (d-1) restanti assi, con due possibili direzioni rispetto a quello asse.
Così, il secondo orientamento vettore ha (2 * (d-1)) opzioni.
Allo stesso modo, la prossima orientamento vettore ha (2 * (d-2)) opzioni.
Il vettore di orientamento finale è costretto a finale restanti asse, ed è costretto a uno dei due possibili direzioni lungo tale asse (per conservare la parità dei set di orientamento vettori; determinante di un ordine di vettori l'orientamento deve essere costante ).
Così, l'ultimo orientamento vettore ha solo (1) opzione.