Lineaire Algebra
Colin Fahey
1. Software
LinearAlgebra.zip
Lineaire algebra broncode (C#)
19910 bytes
MD5: 11d8c8035cac30ba543e5e0b72ee9767
2. Inleiding
Dit artikel beschrijft de vectoren en matrices in (d)-dimensionale ruimte.
3. (d)-dimensionale ruimte: attributen
3.1 Array
„array:“ Een verzameling van variabelen zodanig dat elke variabele heeft een unieke naam, en zodanig te zijn dat de namen kunnen worden toegewezen van een order.
Integer waarden kunnen worden gebruikt als namen voor de variabelen in een array.
Als bijvoorbeeld een array bevat (d) variabelen, dan is de gehele waarden { 0, 1, 2, ..., (d-1) } kunnen de namen die aan de variabelen in de array.
3.2 (d)-dimensionale vector
„(d)-dimensionale vector:“ Een array van (d) variabelen.
„vector component:“ Een variabele binnen een vector.
3.3 (d)-dimensionale vectorruimte
„een-dimensionale ruimte:“ De complete set van waarden die kunnen worden opgeslagen door een variabele.
„(d)-dimensionale ruimte:“ De complete set van combinaties van waarden die kunnen worden opgeslagen door een array van (d) variabelen.
De formele definitie van een „vectorruimte:“
Laten (T) een fundamenteel type (voorbeelden: echte nummer, integer, complex getal, rationeel nummer, etc).
Elke variabele van een type heet een „scalar.“
Een „(T)-type (d)-dimensionale vectorruimte“ is de verzameling (S) van (d)-dimensionale vectoren die twee operaties, vector Bovendien (+), en scalaire vermenigvuldiging (*), voldoet aan de voorwaarden hieronder:
(1) Als (v) en (w) zijn twee vectoren in (S), dan (v + w) is ook een vector in (S);
(2) Als (u), (v), en (w) zijn alle drie vectoren in (S), dan (u + v) + w = u + (v + w);
[toevoegingsmiddel Commutativiteit]
(3) Indien (v) en (w) zijn twee vectoren in (S), dan (v + w) = (w + v);
[toevoegingsmiddel Associativiteit]
(4) Er is een „nul-vector,“ (0), in (S), zodanig dat voor een vector (v) in (S), (v + (0)) = v;
[toevoegingsmiddel identiteit]
(5) Als (c) is een scalaire van het type (T), en (v) is het vector in (S), dan is het product (c * v) is een vector in (S);
(6) Als (a), (b), en (c) zijn alle scalars van het type (T), en (v) en (w) zijn alle vectoren in (S), dan (a + b) * v = a*v + b*v, en c*(v + w) = c*v + c*w;
[vermenigvuldigingsmethode Distributiviteit]
(7) Als (a) en (b) zijn alle scalars van het type (T), en (v) is het vector in (S), dan (a*b)*v = a*(b*v);
(8) Als „1“ is een scalaire van het type (T) zodanig zijn dat (1*1)=1, en (v) is het vector in (S), dan (1*v) = v;
(9) Voor elke vector (v) in (S), de vector (-1)*v = -v voldoet aan v + (-v) = (0);
[additieve inverse]
3.4 (d)-dimensionale vector-code
De code hieronder toont hoe een (d)-dimensionale vector, met 64-bit floating-point componenten, kunnen worden uitgevoerd.
Een array is voldoende om een vector.
De onderstaande code is een array in een klasse, maar voor het gemak.
De code is niet bedoeld als efficiënt.
Een structuur (bijvoorbeeld: „struct“, een waarde type) van vectoren van een vast aantal dimensies (voorbeelden: 3 of 4) wordt waarschijnlijk veel efficiënter dan de algemene (d)-dimensionale vector klasse hier weergegeven.
Hoewel de onderstaande code definieert een vector met floating-point componenten, dit document maakt ook gebruik van vectoren met integer componenten.
De onderstaande code kan gemakkelijk worden aangepast voor de uitvoering van vectoren met integer componenten.
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 nul vector:“ een vector met alle componenten (d) gelijk is aan nul.
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 ruimte: de punten
„(d)-dimensionale ruimte punt:“ Een array van (d) variabelen met specifieke waarden „(coördinatie).“
„(d)-dimensionale ruimte oorsprong:“ Een array van (d) variabelen met alle waarden die gelijk is aan nul.
3.6 (d)-dimensionale vectoren: niet-relatieve en relatieve
Een (d)-dimensionale vector die een „niet-familielid“ is een (d)-dimensionale vector, dat direct is een status of configuratie.
Een (d)-dimensionale vector is dat „relatief“ is een (d)-dimensionale vector is dat wijzigingen in een aantal componenten.
Een relatieve vector kan het verschil tussen twee niet-relatieve vectoren.
Aangezien een relatieve vector, de bepaling van een niet-relatieve status of configuratie met behulp van dat relatieve vector vereist een combinatie van dat relatieve vector met een niet-relatieve vector.
Niet-relatieve vectoren en relatieve vectoren zijn beide vectoren.
Of een bepaalde vector is niet verwant of familielid moet worden bepaald wanneer de vector is gedefinieerd.
Als een (d)-dimensionale vector is geïnterpreteerd als niet-familielid, dan is de (d)-dimensionale vector kan een punt in (d)-dimensionale ruimte.
3.7 (d)-dimensionale vector optellen, aftrekken, en schaal -
Vector optellen, aftrekken, en schaal:
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 basis vectoren
De formele definitie van een „basis“ van een vectorruimte:
Laten (T) een fundamenteel type (voorbeelden: echte nummer, integer, complex getal, rationeel nummer, etc).
Elke variabele van een type heet een „scalar.“
Laten (V) een „(T)-type (d)-dimensionale vectorruimte.“
Als de niet-nul vectoren { u1, u2, ..., ud } in (V) dien aard zijn dat elke vector (v) in (V) kan worden geschreven als een „lineaire combinatie“ van de vectoren, v = c1*u1 + c2*u2 + ... + cd*ud, waar { c1, c2, ..., cd } scalars zijn van het type (T), dan (V) wordt „overspannen“ door vectoren { u1, u2, ..., ud }.
Elke set van niet-nul vectoren { u1, u2, ..., ud } dat „span“ vectorruimte (V) is genoemd op „basis“ van (V).
Een eenvoudige „basis“ van een (d)-dimensionale vectorruimte is de verzameling van (d) onderscheiden (d)-dimensionale vectoren, elk met een element dat gelijk is aan een en alle andere componenten die gelijk is aan nul.
Deze basis-vectoren zijn „orthonormal,“ wat betekent dat zij elkaar over-loodrecht „(loodrecht)“ en dat elke vector heeft lengte-eenheid.
Elk van deze vector is een vector eenheid parallel aan een van de zwaartepunten (d) coördineren.
De uiting van een willekeurige vector als een „lineaire combinatie“ van deze basis vectoren is direct; elke component van de willekeurige vector wordt vermenigvuldigd met de overeenkomstige basis vector, en deze producten worden bij elkaar opgeteld om de willekeurige vector.
De onderstaande code laat zien hoe een vector kan worden uitgedrukt als een „lineaire combinatie“ van basis vectoren.
De onderstaande code willekeurig definieert een set van orthonormal basis vectoren.
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 )
// . . .
}
}
Elke vector in (d)-dimensionale ruimte kan worden uitgedrukt als een som van de producten van nummers en basis vectoren:
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 ruimte: de afstand tussen de punten
Laten (P) een (d)-dimensionale vector is een punt dat in (d)-dimensionale ruimte.
Laten (Q) een (d)-dimensionale vector is een punt dat in (d)-dimensionale ruimte.
Laten (R) een (d)-dimensionale vector, dat komt overeen met de verandering in (d) coördinaten te krijgen (P) van punt tot punt (Q); R = (Q - P).
In een-dimensionale ruimte, P = ( p0 ), Q = ( q0 ), en R = (Q - P) = ( q0 - p0 ).
De afstand tussen de twee punten is: Abs( q0 - p0 ).
In twee-dimensionale ruimte, P = ( p0, p1 ), Q = ( q0, q1 ), en R = (Q - P) = ( q0-p0, q1-p1 ).
Het interpreteren van de twee loodrecht staat is als het loodrecht zijden van een recht-driehoek, de afstand tussen de punten overeen met de lengte van de schuine zijde van die driehoek.
De Pythagorean formule, (a*a) + (b*b) = (c*c), waar (a) en (b) zijn de lengtes van de zijden loodrecht op een recht-driehoek, en (c) is de lengte van de schuine zijde (kant beïnvloeden), kan worden gebruikt voor het bepalen van de afstand tussen de twee punten: Sqrt( Sq(q0-p0) + Sq(q1-p1) ).
In de driedimensionale ruimte, P = ( p0, p1, p2 ), Q = ( q0, q1, q2 ), en R = (Q - P) = ( q0-p0, q1-p1, q2-p2 ).
Het interpreteren van de verplaatsing ( q0-p0, q1-p1, 0 ) als het loodrecht zijden van een recht-driehoek, en met behulp van de Pythagorean formule, geeft de afstand tussen punt (P) en het punt ( q0, q1, p2 ): d01 = Sqrt( Sq(q0-p0) + Sq(q1-p1) ).
De verplaatsing ( q0-p0, q1-p1, 0 ) staat loodrecht op de verplaatsing ( 0, 0, q2-p2 ), en een andere rechts-driehoek kan worden gevormd, en de Pythagorean formule kan worden gebruikt.
Zo is de afstand van punt (P) punt (Q) wordt gegeven door: 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) ).
De methode van de uitbreiding van de afstand formule van twee-dimensionale ruimte aan drie-dimensionale ruimte kan worden toegepast bij herhaling tot uiteindelijk het bepalen van de afstand formule voor (d)-dimensionale ruimte: Sqrt( (Sq(q0-p0) + Sq(q1-p1)) + Sq(q2-p2) + ... + Sq(qd-pd) ).
De onderstaande code definieert een functie genaamd „Lengte“ die berekent de lengte van een (d)-dimensionale vector.
Wanneer een vector is de verplaatsing tussen twee punten in (d)-dimensionale ruimte, de lengte van de vector is de afstand tussen die twee punten.
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 vectoren: dot product
De „dot product“ zet (d) twee-dimensionale vectoren aan een nummer.
De onderstaande code continu een dot product van twee vectoren:
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);
}
// . . .
}
Voor het vector, (A), Length(A) = Sqrt(Dot(A,A)).
3.11 (d)-dimensionale vectoren: definitie van „parallelle“
Vectoren (A) en (B) zijn „parallel“ als aan alle volgende waar zijn:
(1) A.Length() > 0;
(2) B.Length() > 0;
(3) Abs(Dot(A,B)) = A.Length()*B.Length().
De onderstaande code bepaalt of een combinatie van de vectoren parallel (eventueel anti-gebonden).
Floating-point getallen kunnen zich fouten in de fractionele deel te wijten aan de beperkte precisie, en dus ook computer-code moet ook niet-nul-toleranties bij de vergelijking van floating-point getallen.
De code omvat bijvoorbeeld de waarden van tolerantie, maar het voorbeeld van tolerantie waarden wellicht niet geschikt zijn voor bepaalde taken.
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 vectoren: definitie van „loodrecht“
Vectoren (A) en (B) zijn „loodrecht“ als aan alle volgende waar zijn:
(1) A.Length() > 0;
(2) B.Length() > 0;
(3) Abs(Dot(A,B)) = 0.
De onderstaande code bepaalt of een combinatie van de vectoren loodrecht. Floating-point getallen kunnen zich fouten in de fractionele deel te wijten aan de beperkte precisie, en dus ook computer-code moet ook niet-nul-toleranties bij de vergelijking van floating-point getallen.
De code omvat bijvoorbeeld de waarden van tolerantie, maar het voorbeeld van tolerantie waarden wellicht niet geschikt zijn voor bepaalde taken.
public class VectorF64
{
// . . .
public static bool Perpendicular( VectorF64 a, VectorF64 b )
{
// Smallest normalized float : 1.1754943e-38
// Smallest normalized double: 2.2250738585072020e-308
double nonZeroThreshold = 1.0e-38; // conservative for double
// double: (52+1)-bit mantissa; log10(2^53)=15.95 decimal digits
double fractionalDifferenceThreshold = 1.0e-14; // conservative
if ((null == a) || (null == b))
{
return (false); // Vector is not specified.
}
if ((null == a.components) || (null == b.components))
{
return (false); // Vector is empty.
}
if (a.Dimensions() != b.Dimensions())
{
return (false); // Vectors not the same size.
}
double lengthA = a.Length();
if (lengthA <= nonZeroThreshold)
{
return (false);
}
double lengthB = b.Length();
if (lengthB <= nonZeroThreshold)
{
return (false);
}
double zeroImpliesPerpendicular =
Math.Abs( VectorF64.Dot( a, b ) ) / (lengthA * lengthB);
double absoluteDifferenceFromZero =
Math.Abs( zeroImpliesPerpendicular );
if (absoluteDifferenceFromZero <= fractionalDifferenceThreshold)
{
return (true);
}
return (false);
}
// . . .
public static void Test()
{
// . . .
// Example of testing for perpendicular vectors:
// A 6-dimensional vector:
VectorF64 vf2 = new VectorF64( 0.0, 1.0, 2.0, 0.0, 4.0, 5.0 );
vf2.WriteLine(); // ( 0, 1, 2, 0, 4, 5 )
// A 6-dimensional vector:
VectorF64 vg2 = new VectorF64( 10.0, 0.0, 0.0, -5.0, 0.0, 0.0 );
vg2.WriteLine(); // ( 10, 0, 0, -5, 0, 0 )
// Determine if the specified vectors are perpendicular
bool perpendicular = VectorF64.Perpendicular( vf2, vg2 );
Log.WriteLine( perpendicular ); // True
// Add a non-negligible displacement to a component of a vector:
vf2[0] += 1.0e-13;
vf2.WriteLine(); // ( 1E-13, 1, 2, 0, 4, 5 )
// Determine if the specified vectors are perpendicular
perpendicular = VectorF64.Perpendicular( vf2, vg2 );
Log.WriteLine( perpendicular ); // False
// . . .
}
}
3.13 Matrices
„matrix:“ Een verzameling van variabelen zodanig dat elke variabele heeft een unieke combinatie van een „rij“ en een „kolom“ naam.
„entry:“ Een variabele binnen een matrix.
Integer waarden kunnen worden gebruikt als namen van „rij“ en „kolom“ namen voor de variabelen in een matrix.
Bijvoorbeeld, als een matrix heeft (totalRows) rijen en kolommen (totalColumns), dan is de gehele waarden { 0, 1, ..., (totalRows-1) } kunnen de namen die aan de rijen, en de gehele waarden { 0, 1, ..., (totalColumns-1) } kunnen de namen die aan de kolommen.
Zo heeft een variabelen in een matrix kan worden aangegeven door het opgeven van een paar getallen, ( row, column ), met vermelding van de combinatie van rij en kolom indices in overeenstemming met de specifieke variabele.
De grootte van een matrix is opgegeven als „(totalRows) * (totalColumns)“ (of „(totalRows) door (totalColumns)).“
Deze volgorde van de afmetingen is dezelfde als de volgorde van de dimensies gebruikt worden om items in de matrix „(( row, column )).“
[Dit verdrag is een beetje jammer, want voor vele twee-dimensionale gebruikt (voorbeelden: beelden, grafieken, etc) de gebruikelijke gang van zaken is te specificeren afmetingen als „width * height“ en coördineert als „( horizontal, vertical )“ (of „( x, y )).“
Dit is het tegenovergestelde van de orde van de afmetingen en de coördinaten worden gebruikt voor de omschrijving matrices en hun inzendingen. ]
Een matrix met (totalRows) gelijk aan (totalColumns) is genoemd „plein,“ anders is de matrix is vernoemd „rechthoekig.“
Een matrix kan worden beschouwd als met een reeks van „rij-vectoren,“ waar de variabelen in elke rij worden geïnterpreteerd als behorend tot een vector.
Een matrix kan ook worden beschouwd als met een reeks „kolom vectoren,“ waar de variabelen in elke kolom worden geïnterpreteerd als behorend tot een vector.
Matrices kan vertegenwoordigen een breed scala van wiskundige relaties.
De betekenis van een matrix, en de concrete acties die kunnen worden geschikt voor de verwerking van gegevens van een matrix, hangt af van de context.
Er zijn echter fundamentele regels van de matrix rekenkundig die relevant zijn voor vele contexten, en deze fundamentele regels worden vastgelegd in een volgende sectie.
Een array en een (totalColumns) waarde, voldoende zijn om een matrix.
De array kan hebben (totalRows * totalColumns) variabelen, en de toegang aan ( row, column ) kunnen beantwoorden aan de array variabele met index ((totalColumns * row) + column).
De onderstaande code wordt een matrix, met 64-bit floating-point entries.
Een array en een (totalColumns) waarde, voldoende zijn om een matrix.
De onderstaande code is een array en een (totalColumns) waarde in een klasse, maar voor het gemak.
De code is niet bedoeld als efficiënt.
Een structuur (bijvoorbeeld: „struct“, een waarde type) die matrices van bepaalde vaste afmetingen (voorbeelden: 2*2, 3*3, of 4*4) wordt waarschijnlijk veel efficiënter dan de algemene (totalRows * totalColumns) klasse hier weergegeven.
Hoewel de onderstaande code definieert een matrix met floating-point entries, dit artikel maakt ook gebruik van matrices met integer inzendingen.
De onderstaande code kan eenvoudig worden aangepast voor de uitvoering van matrices met integer inzendingen.
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
// . . .
}
}
„nul-matrix:“ een matrix met alle items gelijk aan nul.
public class MatrixF64
{
// . . .
public static MatrixF64 Zero( int rows, int columns )
{
MatrixF64 zero = new MatrixF64( rows, columns );
// The constructor above sets all entries to zero.
return (zero);
}
// . . .
public static void Test()
{
// . . .
// An 8x2 matrix (8 rows x 2 columns) with all 16
// 64-bit floating-point entries set to zero:
MatrixF64 z = MatrixF64.Zero( 8, 2 );
z.WriteLine();
// [ 0, 0 ] // row 0
// [ 0, 0 ] // row 1
// [ 0, 0 ] // row 2
// [ 0, 0 ] // row 3
// [ 0, 0 ] // row 4
// [ 0, 0 ] // row 5
// [ 0, 0 ] // row 6
// [ 0, 0 ] // row 7
// . . .
}
}
3.14 Matrix optellen, aftrekken, vermenigvuldigen en
Matrix optellen, aftrekken, vermenigvuldigen.
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 ]
// . . .
}
}
„identiteit matrix:“ een vierkante matrix (totaal rijen is gelijk aan de totale kolommen) met de tekst gelijk is aan (1) op de diagonaal (met een rij-index gelijk aan de kolom index), en met alle andere items die gelijk is aan (0).
Als een vierkante matrix, (M), wordt vermenigvuldigd met een „identiteit matrix,“ (I), van hetzelfde aantal rijen (en kolommen), het product is gelijk aan (M).
Vermenigvuldigen (I) door (M) ook een product op de markt die gelijk is aan (M).
Zo is de „identiteit matrix“ is analoog aan het aantal „1“ voor de vermenigvuldiging van getallen (scalars).
De onderstaande code zorgt voor een identiteit matrix met een bepaald aantal rijen.
public class MatrixF64
{
// . . .
public static MatrixF64 Identity( int rows )
{
MatrixF64 identity = MatrixF64.Zero( rows, rows );
for (int i = 0; i < rows; i++)
{
identity[i, i] = 1.0;
}
return (identity);
}
// . . .
public static void Test()
{
// . . .
// Examples of multiplying by an identity matrix:
MatrixF64 identity3x3 = MatrixF64.Identity( 3 );
// [ 1, 0, 0 ]
// [ 0, 1, 0 ]
// [ 0, 0, 1 ]
MatrixF64 s3x3 =
new MatrixF64
(
3, 3,
0.0, 1.0, 2.0, // row 0
3.0, 4.0, 5.0, // row 1
6.0, 7.0, 8.0 // row 2
);
result = s3x3 * identity3x3;
result.WriteLine();
// [ 0, 1, 2 ]
// [ 3, 4, 5 ]
// [ 6, 7, 8 ]
// (the same as s3x3)
result = identity3x3 * s3x3;
result.WriteLine();
// [ 0, 1, 2 ]
// [ 3, 4, 5 ]
// [ 6, 7, 8 ]
// (the same as s3x3)
// . . .
}
}
3.15 Matrix: LU factoring; back-substitutie
„Matrix LU factoring:“ Een procedure die converteert een vierkante matrix, (M), om twee nieuwe vierkante matrices, (L) en (U), met dezelfde grootte als (M), zodanig dat (L) * (U) = (M), en dat deze matrix is (U) „bovenste driehoekige“ (alle items onder de diagonaal zijn nul), en de matrix (L) „lager“ is „driehoekig“ (alle items boven de diagonaal zijn nul).
Matrix LU factoring kunnen worden gebruikt als onderdeel van een grotere procedure voor het oplossen van een stelsel van vergelijkingen, of bij het vinden van de inverse van een matrix, of bij het vinden van de determinant van een matrix.
Matrix LU factoring is een alternatief voor de afschaffing „Gauss-Jordan procedure.
Gauss-Jordan afschaffing vereist een stelsel van vergelijkingen (A)*(x)=(b), terwijl LU factoring alleen vereist een matrix (A).
Ook na de vaststelling van de LU factoring van een matrix (A), is het erg eenvoudig om te bepalen (x) gegeven alle (b).
De procedure voor het oplossen van de vector (x) in (A)*(x)=(b), gezien (b) en de LU factoring (L)*(U)=(A), een procedure die zal worden vernoemd „back-substitutie“ hier.
De onderstaande code bevat een procedure voor de berekening van de LU factoren van een plein, niet-singuliere matrix.
De onderstaande code bevat een procedure voor het doen LU back-substitutie.
Let op: De computer berekent de onderstaande code L en U factoren van een rij-toe-versie van een bepaalde matrix.
Ook de L en U matrix resultaten zijn gebundeld in een enkele output matrix.
public class MatrixF64
{
// . . .
public static void FindLUFactorsOfARowPermutedVersionOfAnOriginalMatrix
(
MatrixF64 originalMatrix,
ref MatrixF64 LUFactorsMatrix,
ref int[] rowPermutationOfOriginalMatrix,
ref bool rowExchangeParityIsOdd
)
{
// ''LU factoring'' involves factoring a square matrix (M) in to a
// lower-triangular matrix (L), and an upper-triangular matrix (U),
// such that (L)*(U) = (M).
//
// However, the specified matrix (originalMatrix) might not be
// optimal for direct LU factoring, due to the locations of extreme
// values in the matrix. Therefore, this function instead implicitly
// factors a row-permuted version of originalMatrix, where the
// permutation of rows is determined by the values in the specified
// matrix. The specific row permutation selected by this function is
// returned in an array of row indices (rowPermutationOfOriginalMatrix).
// The resulting (L) and (U) factors are merged in to a single
// output matrix (LUFactorsMatrix). Therefore:
//
// (L part of LUFactorsMatrix) * (U part of LUFactorsMatrix)
// = (originalMatrix with rows permuted according to
// rowPermutationOfOriginalMatrix).
//
// Although factoring a row-permuted version of originalMatrix
// makes it more difficult to interpret the results of this function,
// the resulting LU factoring is likely to be more accurate. In a
// sense, this function indicates (via rowPermutationOfOriginalMatrix)
// how to permute the rows of originalMatrix to be able to produce
// the most accurate LU factoring directly, and this function produces
// that factoring (via the LUFactorsMatrix output).
//
// When using the (L) and (U) matrices (merged in to LUFactorsMatrix)
// to solve a system of equations, it is necessary to permute the rows
// of the solution vector of the system of equations according to
// rowPermutation.
//
// The output matrix (LUFactorsMatrix) is a combination of
// the (L) and (U) matrices:
//
// To form the (L) matrix from (LUFactorsMatrix):
// assume (0) for entries above the diagonal,
// assume (1) for entries on the diagonal,
// and use (LUFactorsMatrix) entries below the diagonal.
//
// To form the (U) matrix from (LUFactorsMatrix):
// assume (0) for entries below the diagonal,
// and use (LUFactorsMatrix) entries on and above the diagonal.
//
// The output array (rowPermutationOfOriginalMatrix) indicates how the
// rows of the original matrix have implicitly (not actually) been
// permuted.
//
// The output boolean (rowExchangeParityIsOdd) indicates if the parity
// of the permutation of rows is odd.
//
// This implementation is partly based on pseudo-code appearing in
// ''Introduction to Algorithms'' (Cormen; Lieserson; Rivest;
// 24th printing, 2000; MIT Press; ISBN 0-262-03141-8).
// See the section named ''Overview of LUP decomposition'', in
// chapter 31 (''Matrix Operations''). The implementation here follows
// the ''LUP-Decomposition'' pseudo-code appearing in a section named
// ''Computing an LUP decomposition''.
//
// This implementation also uses ideas found in the implementation of
// ''ludcmp'' in the book ''Numerical recipes in C : The art of
// scientific computing'' (Press; Teukolsky; Vetterling; Flannery;
// second edition; 1996 reprinting; Cambridge University Press).
// The overall loop structure here resembles that of ''ludcmp''.
// The idea of determining and caching row scale factors in advance of
// finding pivots is used here. The method in ''ludcmp'' to avoid zero
// pivot values has been modified here to handle excessively-small
// pivot values, too. HOWEVER, the row permutation indices produced by
// the following implementation is a true permutation of
// {0,1,...,(totalRows-1)}, whereas ''ludcmp'' (and ''lubksb'')
// interpret their ''permutation indices'' as swaps to perform
// (roughly, for each (i), swap b[ i ] and b[ indx[i] ]). So, the
// following implementation of LU-factoring is not directly compatible
// with ''ludcmp'' or ''lubksb''. Also, the following procedure
// does not destroy the original matrix (unlike ''ludcmp'').
double singularMatrixIfMaxRowElementIsLessThanThis = (1.0e-19);
double forcePivotsToHaveAtLeastThisAbsoluteValue = (1.0e-19);
LUFactorsMatrix = null;
rowPermutationOfOriginalMatrix = null;
rowExchangeParityIsOdd = false;
if (null == originalMatrix)
{
return; // No matrix specified.
}
if
(
(originalMatrix.totalRows <= 0)
|| (originalMatrix.totalColumns <= 0)
|| (null == originalMatrix.entries)
)
{
return; // Matrix is empty.
}
if (originalMatrix.totalRows != originalMatrix.totalColumns)
{
return; // Matrix is not square.
}
// Duplicate the original matrix
LUFactorsMatrix = new MatrixF64( originalMatrix );
// (Lines 2-3 of LUP-Decomposition)
// Initialize the row permutation array to the identity
// permutation.
rowPermutationOfOriginalMatrix = new int[LUFactorsMatrix.totalRows];
for (int i = 0; i < LUFactorsMatrix.totalRows; i++)
{
rowPermutationOfOriginalMatrix[i] = i;
}
// For each row, determine the largest absolute value of
// the row elements, and use the reciprocal as the scale for
// that row.
VectorF64 rowScales = VectorF64.Zero( LUFactorsMatrix.totalRows );
for (int i = 0; i < LUFactorsMatrix.totalRows; i++)
{
double largestElementAbsoluteValueInRow = 0.0;
for (int j = 0; j < LUFactorsMatrix.totalColumns; j++)
{
double absoluteElementValue =
Math.Abs( LUFactorsMatrix[i, j] );
if (absoluteElementValue >
largestElementAbsoluteValueInRow)
{
largestElementAbsoluteValueInRow =
absoluteElementValue;
}
}
if (largestElementAbsoluteValueInRow <
singularMatrixIfMaxRowElementIsLessThanThis)
{
return; // Matrix is singular
}
rowScales[i] = (1.0 / largestElementAbsoluteValueInRow);
}
// (Lines 4-18 of LUP-Decomposition)
// Go through the columns of the matrix.
for (int j = 0; j < LUFactorsMatrix.totalColumns; j++)
{
// (Lines 17-18 of LUP-Decomposition)
// (or equation (2.3.12) of NR)
for (int i = 0; i < j; i++)
{
double sum = LUFactorsMatrix[i, j];
for (int k = 0; k < i; k++)
{
sum -= (LUFactorsMatrix[i, k] * LUFactorsMatrix[k, j]);
}
LUFactorsMatrix[i, j] = sum;
}
// Go through all remaining rows and find the best pivot.
double largestScaledSum = 0.0;
int rowIndexOfLargestScaledSum = j;
for (int i = j; i < LUFactorsMatrix.totalRows; i++)
{
// (equation (2.3.13) of NR)
double sum = LUFactorsMatrix[i, j];
for (int k = 0; k < j; k++)
{
sum -= (LUFactorsMatrix[i, k] * LUFactorsMatrix[k, j]);
}
LUFactorsMatrix[i, j] = sum;
double scaledSum = rowScales[i] * Math.Abs( sum );
if (scaledSum >= largestScaledSum)
{
largestScaledSum = scaledSum;
rowIndexOfLargestScaledSum = i;
}
}
// If indeed we found a better pivot, then exchange rows.
if (j != rowIndexOfLargestScaledSum)
{
// (Line 12 of LUP-Decomposition)
// Exchange the row permutation indices
int tempRowIndex = rowPermutationOfOriginalMatrix[j];
rowPermutationOfOriginalMatrix[j] =
rowPermutationOfOriginalMatrix[rowIndexOfLargestScaledSum];
rowPermutationOfOriginalMatrix[rowIndexOfLargestScaledSum] =
tempRowIndex;
// (Lines 13-14 of LUP-Decomposition)
// Exchange the elements of the rows
for (int k = 0; k < LUFactorsMatrix.totalColumns; k++)
{
double temp = LUFactorsMatrix[rowIndexOfLargestScaledSum, k];
LUFactorsMatrix[rowIndexOfLargestScaledSum, k] =
LUFactorsMatrix[j, k];
LUFactorsMatrix[j, k] = temp;
}
// Exchange the row scale factors
double scaleFactor =
rowScales[rowIndexOfLargestScaledSum];
rowScales[rowIndexOfLargestScaledSum] = rowScales[j];
rowScales[j] = scaleFactor;
// Invert the overall row exchange parity
rowExchangeParityIsOdd = (!(rowExchangeParityIsOdd));
}
// Force the pivot element to have at least a certain
// absolute value.
if (Math.Abs( LUFactorsMatrix[j, j] ) <
forcePivotsToHaveAtLeastThisAbsoluteValue)
{
if (LUFactorsMatrix[j, j] < 0.0)
{
LUFactorsMatrix[j, j] =
(-(forcePivotsToHaveAtLeastThisAbsoluteValue));
}
else
{
LUFactorsMatrix[j, j] =
forcePivotsToHaveAtLeastThisAbsoluteValue;
}
}
// If not the final column, then divide all column elements
// below the diagonal by the pivot element, matrixLU[j,j].
// (Lines 15-16 of LUP-Decomposition)
if (j != (LUFactorsMatrix.totalColumns - 1))
{
double reciprocalOfPivot = (1.0 / LUFactorsMatrix[j, j]);
for (int i = (j + 1); i < LUFactorsMatrix.totalRows; i++)
{
LUFactorsMatrix[i, j] *= reciprocalOfPivot;
}
}
}
}
public static void LUBacksubstitution
(
MatrixF64 LUFactorsMatrix,
int[] rowPermutationOfOriginalMatrix,
VectorF64 givenProductVector,
ref VectorF64 solutionVector
)
{
// This implementation is based on pseudo-code appearing in
// ''Introduction to Algorithms'' (Cormen; Lieserson; Rivest;
// 24th printing, 2000; MIT Press; ISBN 0-262-03141-8).
// See the section named ''Overview of LUP decomposition'', in
// chapter 31 (''Matrix Operations''). The implementation here
// follows the ''LUP-Solve'' pseudo-code appearing in a
// section named ''forward and back substitution''.
solutionVector = null;
if (null == LUFactorsMatrix)
{
return; // No matrix specified.
}
if (null == rowPermutationOfOriginalMatrix)
{
return; // No permutation specified.
}
if (null == givenProductVector)
{
return; // No product vector specified.
}
if
(
(null == LUFactorsMatrix.entries)
|| (LUFactorsMatrix.totalRows <= 0)
|| (LUFactorsMatrix.totalColumns <= 0)
)
{
return; // Matrix is empty.
}
if (LUFactorsMatrix.totalRows !=
LUFactorsMatrix.totalColumns)
{
return; // Matrix is not square.
}
if (givenProductVector.Dimensions() !=
LUFactorsMatrix.totalRows)
{
return; // Product vector size not equal to matrix row count.
}
// Copy the product vector in to the result.
solutionVector = new VectorF64( givenProductVector );
// (See LUP-Solve, lines 2-3)
for (int i = 0; i < LUFactorsMatrix.totalRows; i++)
{
double sum =
givenProductVector[rowPermutationOfOriginalMatrix[i]];
for (int j = 0; j < i; j++)
{
sum -= (LUFactorsMatrix[i, j] * solutionVector[j]);
}
solutionVector[i] = sum;
}
// (See LUP-Solve, lines 4-5)
for (int i = (LUFactorsMatrix.totalRows - 1); i >= 0; i--)
{
double sum = solutionVector[i];
for
(
int j = (i + 1);
j < LUFactorsMatrix.totalColumns;
j++
)
{
sum -= (LUFactorsMatrix[i, j] * solutionVector[j]);
}
double diagonalElement = LUFactorsMatrix[i, i];
if (diagonalElement != 0.0)
{
solutionVector[i] = (sum / diagonalElement);
}
}
}
// . . .
public static void Test()
{
// . . .
// Example of LU factoring and back-substitution:
MatrixF64 a3x3 =
new MatrixF64
(
3, 3,
0.0, 1.0, 2.0, // row 0
3.0, 4.0, 5.0, // row 1
6.0, 7.0, 8.0 // row 2
);
int[] rowPermutationOfOriginalMatrix = null;
bool rowInterchangeParityIsOdd = false;
MatrixF64 a3x3LU = null;
MatrixF64.FindLUFactorsOfARowPermutedVersionOfAnOriginalMatrix
(
a3x3,
ref a3x3LU,
ref rowPermutationOfOriginalMatrix,
ref rowInterchangeParityIsOdd
);
a3x3LU.WriteLine();
// [ 6, 7, 8 ]
// [ 0, 1, 2 ]
// [ 0.5, 0.5, 1e-19 ]
// The matrix above is a combination of the L and U
// matrices. The L matrix is 0 above the diagonal,
// 1 on the diagonal, and the shown values below the
// diagonal. The U matrix is 0 below the diagonal,
// and the shown values on and above the diagonal.
// This LU factoring is of a ROW-PERMUTED version
// of the original matrix.
// The following L and U matrices are manually formed
// by inspecting the LU factoring matrix above.
MatrixF64 L3x3 =
new MatrixF64
(
3, 3,
1.0, 0.0, 0.0, // row 0
0.0, 1.0, 0.0, // row 1
0.5, 0.5, 1.0 // row 2
);
MatrixF64 U3x3 =
new MatrixF64
(
3, 3,
6.0, 7.0, 8.0, // row 0
0.0, 1.0, 2.0, // row 1
0.0, 0.0, 0.0 // row 2
);
// The product (L)*(U) produces a ROW-PERMUTED version
// of the original matrix:
result = L3x3 * U3x3;
result.WriteLine( 4 );
// [ 6, 7, 8 ]
// [ 0, 1, 2 ]
// [ 3, 4, 5 ]
VectorF64 b3x1 = new VectorF64( 1.0, 2.0, 3.0 );
VectorF64 x3x1 = null;
MatrixF64.LUBacksubstitution
(
a3x3LU,
rowPermutationOfOriginalMatrix,
b3x1,
ref x3x1
);
// Solution vector
x3x1.WriteLine();
// ( -0.66666667, 1, 0 )
// Both the given b3x1 and this x3x1 are relative to the
// original matrix, a3x3. The a3x3LU for a ROW-PERMUTED
// version of a3x3, but LUBacksubstitution accepts a
// non-permuted product vector (e.g., b3x1) and produces
// a non-permuted solution vector (e.g., x3x1).
// Check solution by multiplying the original matrix a3x3
// by the solution vector x3x1, to get a product vector.
VectorF64 checkb3x1 = a3x3 * x3x1;
checkb3x1.WriteLine();
// ( 1, 2, 3 )
// This matches the b3x1 vector we specified above, so
// the solution is valid.
// . . .
}
}
3.16 Matrix: determinant
„determinant:“ Een functie die converteert elke vierkante (n * n) matrix (A) aan een aantal, det(A), zodanig zijn dat:
(1) Als een matrix (B) resultaten uit te wisselen twee rijen of twee kolommen van een matrix (A), dan det(B) = (-(det(A)));
(2) Als een matrix (B) resultaten van de vermenigvuldiging van elke rij, kolom of een, van een matrix (A), door een aantal (c), dan det(B) = c * det(A);
(3) Als een matrix (B) resultaten van het toevoegen van een veelvoud van een rij naar de andere rij, of door het toevoegen van een veelvoud van een kolom naar de andere kolom, dan det(B) = det(A).
(4) Als (A) is een boven-driehoekige matrix (A[i,j]=0 voor alle (i>j)) of een kleine driehoekige matrix (A[i,j]=0 voor alle (i<j)), dan det(A) = (A[0,0] * A[1,1] * ... * A[n-1,n-1]);
(Bijvoorbeeld de determinant van een matrix is een identiteit; det(I) = 1.)
Regels (1), (2), en (3), kunnen worden gebruikt, in een proces genaamd „Gaussian uitschakeling, het omzetten van elke vierkante matrix met een driehoekige matrix.
Regel (4) kan worden gebruikt om de bepalende factor van een driehoekstransacties matrix.
De formele definitie van een „determinant:“
Laten M[n](K) aanduiding van de set van alle (n * n) matrices over het gebied (K).
De „bepalende factor“ is de unieke functie (F) zodanig zijn dat F:M[n](K) --> K met twee attributen:
(1) F is afwisselend multilinear met betrekking tot de kolommen (of rijen) en,
(2) F( I ) = 1.
Een „multilinear“ functie (D) van een (n * n) matrix (A) kan worden geschreven als: 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]] ) ), waar de som wordt genomen over alle (n^n) combinaties van 0 <= k[i] <= (n-1), en waar e[j] vertegenwoordigt rij (j) van de identiteit matrix.
Een dergelijke functie kan worden gekenmerkt door de waarden van D( e[k[0]], e[k[1]], ..., e[k[n-1]] ).
Een „afwisselend multilinear“ functie is een functie multilinear, (D), zodanig dat D( ..., e[i], ..., e[i], ... ) = 0, en D( ..., e[i], e[j], ... ) = (-( D( ..., e[j], e[i], ... ) )).
Zo ruilt de kolommen (of rijen), verandert het teken van de functie.
Ook de term D( e[k[0]], ..., e[k[n-1]] ) is nul als de combinatie van k[i] waarden niet alle unieke waarden, als het geheel van waarden is geen permutatie van de nummers {0,1,...,(n-1)}.
Verhuur D( I ) = 1 (D( e[0], e[1], ..., e[n-1] ) = 1), in combinatie met de afwisselende multilinear attributen hierboven beschreven, betekent dat alle waarden van de functie D( e[k[0]], e[k[1]], ..., e[k[n-1]] ) beurten multilinear kan worden bepaald.
De functie wordt niet-nul alleen als het geheel van waarden k[i] is een permutatie van {0,1,...,(n-1)}, en heeft een magnitude van een al is het niet-nul, en het teken zal het aantal swaps verplicht tot het uitvoeren van de permutatie aan de identiteit permutatie.
De waarde van een determinant van een (n * n) matrix (A) kan worden berekend door optelling van de producten van de vorm (sgn(p) * A[0,p[0]] * A[1,p[1]] * ... * A[(n-1),p[n-1]]) voor elke permutatie (p) van de nummers {0,1,2,...,(n-1)}.
De term (sgn(p)) betekent de „ondertekening“ (of swap-count) van permutatie (p), waar (sgn(p)) is gelijk aan (+1) als (p) is een „nog permutatie,“ en is gelijk aan (-1) als (p) is een „oneven permutatie.“
Omdat er (n!) permutaties van de getallen {0,1,2,...,(n-1)}, deze formule heeft (n!) summands.
Andere procedures voor het berekenen van de determinant (voorbeelden: de projecten die een Gaussian eliminatie, of LU factoring, of uitbreiding door minderjarigen, etc) berekenen impliciet een hiërarchie van de tussentijdse hoeveelheden die ervoor zorgen dat die procedures om de bepalende factor voor een bestelling van (n^3) stappen.
Als de determinant van een matrix nul is, de matrix niet over een inverse.
Als de determinant van een matrix is niet-nul, de matrix heeft een inverse.
Als een matrix vormt de coëfficiënten van een lineair systeem van vergelijkingen, en de determinant van de matrix nul is, dan is het systeem van vergelijkingen niet beschikt over een unieke oplossing.
Als een matrix vormt de coëfficiënten van een lineair systeem van vergelijkingen, en de determinant van de matrix is niet nul, dan is het stelsel van vergelijkingen heeft een unieke oplossing.
De onderstaande code bevat de procedures voor het berekenen van de determinant van een vierkante matrix.
Determinanten voor 1*1, 2*2, 3*3, en 4*4 matrices zijn berekend door expliciete formules.
Met behulp van expliciete formules om de determinanten voor kleine matrices is waarschijnlijk om de resultaten sneller dan met behulp van een algemene procedure om deze determinanten.
De procedure voor het algemene geval gebruikt LU factoring, maar er zijn vele andere methoden voor de berekening van determinanten.
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 vectoren: cross product
„cross product:“ een functie die verband houdt met een (d)-dimensionale ruimte dat een geordende lijst van (d-1) (d)-dimensionale vectoren, { A(0), A(1), ..., A(d-2) }, naar een (d)-dimensionale vector, zodanig dat de functie heeft de attributen hieronder:
(1) Als A(i)=0 (nul vector), voor de 0 <= i <= (d-2), dan Cross( A(0), A(1), ..., A(d-2) )=0 (nul vector).
De cross product van een geordende lijst van vectoren is nul als het vector in die lijst is nul;
(2) Indien een paar van vectoren, A(i) en A(j), voor de (i!=j) met 0 <= i <= (d-2) en 0 <= j <= (d-2), zijn parallel (Abs(Dot(A(i),A(j))) = Length(A(i)) * Length(A(j))), dan Cross( A(0), A(1), ..., A(d-2) )=0 (nul vector).
De cross product van een geordende lijst van vectoren is nul als de enige twee vectoren in die lijst zijn parallel;
(3) Als elke vector A(i) is niet nul, en als A(i) is niet parallel met A(j) voor alle (i!=j) met 0 <= i <= (d-2) en 0 <= j <= (d-2), dan B = Cross( A(0), A(1), ..., A(d-2) ) zodanig is dat Dot(A(i),B)=0 voor alle 0 <= i <= (d-2).
Als de cross product van een geordende lijst van vectoren is niet nul, dan is de cross product resultaat is dat loodrecht staat op elke vector in de volgorde lijst van vectoren;
(4) Als (c) is een numerieke waarde, dan 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) ).
Aangezien een geordende lijst van vectoren, de cross product van de lijst in volgorde van de vectoren met een vector vermenigvuldigd met een bepaalde numerieke waarde gelijk is aan de vermenigvuldiging van de cross product van de volgorde lijst van vectoren van de numerieke waarde;
(5) Aangezien twee nummers (i) en (j), waar (i!=j) en 0 <= i <= (d-2) en 0 <= j <= (d-2), en een geordende lijst van (d-1) (d)-dimensionale vectoren, { A(0), A(1), ..., A(d-2) }, en een andere volgorde lijst van (d-1) (d)-dimensionale vectoren, { B(0), B(1), ..., B(d-2) }, zodanig dat B(k) = A(k) voor alle ((k != i) && (k != j)) met 0 <= k <= (d-2), en zodanig te zijn dat B(i) = A(j) en B(j) = A(i), dan Cross( A(0), A(1), ..., A(d-2) ) = (-1) * Cross( B(0), B(1), ..., B(d-2) ).
Als de cross product van een geordende lijst van vectoren is niet nul, dan is dat cross product heeft het tegenovergestelde teken van het kruis van dat product dezelfde volgorde lijst van vectoren, met uitzondering van precies twee vectoren die worden uitgewisseld (swap).
Gelet op basis vectoren e(k), voor 0 <= k <= (d-1), van (d)-dimensionale ruimte, de cross product van een geordende lijst van (d-1) vectoren kan worden berekend door het berekenen van de determinant van de matrix (d*d) hieronder:
e(0), e(1), ..., e(d-1),
A( 0 )[0], A( 0 )[1], ..., A( 0 )[d-1],
A( 1 )[0], A( 1 )[1], ..., A( 1 )[d-1],
...,
A(d-2)[0], A(d-2)[1], ..., A(d-2)[d-1]
De determinant van deze matrix wordt niet een aantal, maar zal in plaats daarvan een vector, waar er zal worden numerieke coëfficiënten voor de (d) basis vectoren.
Voor twee-dimensionale ruimte, en de vector A = ( ax, ay ):
Cross( A )
= Determinant
(
e(x), e(y),
ax, ay
)
= ay * e(x)
- ax * e(y)
= ( ay, -ax )
Het resultaat vector loodrecht staat op de interne input vector.
Als de x-y vliegtuig bekeken is zodanig dat x toeneemt in een rechts richting en y verhogingen in een opwaartse richting, dan is de cross product resultaat is een kwart-draai „de wijzers van de klok“ ten opzichte van de input vector.
Voorbeelden van cross-producten in twee-dimensionale ruimte:
Cross( e(x) ) = (-1) * e(y)
Cross( e(y) ) = e(x)
Voor de driedimensionale ruimte, en de vectoren A = ( ax, ay, az ) en 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) )
Voorbeelden van cross-producten in de driedimensionale ruimte:
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)
Voor vier-dimensionale ruimte, en de vectoren A = ( ax, ay, az, aw ), B = ( bx, by, bz, bw ), en 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) )
Voorbeelden van cross-producten in vier-dimensionale ruimte:
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)
De onderstaande code berekent de cross-product van twee of meer vectoren.
De test voorbeelden tonen aan dat de niet-nul resultaat vectoren loodrecht op de input van alle vectoren.
public class MatrixF64
{
// . . .
public static VectorF64 Cross( params VectorF64[] parameterArray )
{
if (null == parameterArray)
{
return(new VectorF64());
}
int totalVectors = parameterArray.Length;
// The number of dimensions of the space must be one larger
// than the number of specified vectors.
int dimensions = (1 + totalVectors);
if (dimensions < 2)
{
// The cross product does not exist for space with
// fewer than two dimensions.
return (new VectorF64());
}
// All specified vectors must exist and must have a number
// of components equal to the number of dimensions of the space.
for (int k = 0; k < totalVectors; k++)
{
if (null == parameterArray[k])
{
// Vector missing.
return (new VectorF64());
}
else if (null == parameterArray[k].components)
{
// Vector is empty.
return (new VectorF64());
}
else if (parameterArray[k].components.Length != dimensions)
{
// Vector has the wrong number of components.
return (new VectorF64());
}
}
// Form a (d*d) matrix with the (d-1) supplied vectors
// forming the final (d-1) rows of the matrix.
MatrixF64 matrix = MatrixF64.Zero( dimensions, dimensions );
for (int i = 1; i < dimensions; i++)
{
VectorF64 rowVector = parameterArray[ i - 1 ];
for (int j = 0; j < dimensions; j++)
{
matrix[i, j] = rowVector[j];
}
}
VectorF64 result = VectorF64.Zero( dimensions );
// For each component of the result vector, set the first
// row of the matrix to a basis vector, compute the matrix
// determinant, and use the resulting value as the
// component of the result vector. This is inefficient,
// due to the determinant being computed (d) times.
// Clearly it would be better to have explicit formulas
// for 2-dimensional, 3-dimensional, and 4-dimensional
// cross-products.
for (int i = 0; i < dimensions; i++)
{
// Set the first row of the matrix to be the basis
// vector for the coordinate axis with index (i).
for (int j = 0; j < dimensions; j++)
{
matrix[0, j] = 0.0;
}
matrix[0, i] = 1.0;
// Compute the determinant of the modified matrix.
double determinant = matrix.Determinant();
// Set component (i) of the result to the determinant
// that was computed.
result[i] = determinant;
}
return (result);
}
// . . .
public static void Test()
{
// . . .
// Examples of cross products:
// Three-dimensional example:
VectorF64 vcp1x3a = new VectorF64( 1.0, 2.0, 3.0 );
VectorF64 vcp1x3b = new VectorF64( 3.0, 1.0, 2.0 );
VectorF64 vcp1x3 = VectorF64.Cross( vcp1x3a, vcp1x3b );
vcp1x3.WriteLine(); // ( 1, 7, -5 )
// Verify that the result is perpendicular to the original
// vectors:
Log.WriteLine( VectorF64.Perpendicular( vcp1x3, vcp1x3a ) );
Log.WriteLine( VectorF64.Perpendicular( vcp1x3, vcp1x3b ) );
// True, True
// Four-dimensional example:
VectorF64 vcp1x4a = new VectorF64( 1.0, 2.0, 3.0, 4.0 );
VectorF64 vcp1x4b = new VectorF64( 4.0, 1.0, 2.0, 3.0 );
VectorF64 vcp1x4c = new VectorF64( 3.0, 4.0, 1.0, 2.0 );
VectorF64 vcp1x4 = VectorF64.Cross( vcp1x4a, vcp1x4b, vcp1x4c );
vcp1x4.WriteLine(); // ( 4, 4, 44, -36 )
// Verify that the result is perpendicular to the original
// vectors:
Log.WriteLine( VectorF64.Perpendicular( vcp1x4, vcp1x4a ) );
Log.WriteLine( VectorF64.Perpendicular( vcp1x4, vcp1x4b ) );
Log.WriteLine( VectorF64.Perpendicular( vcp1x4, vcp1x4c ) );
// True, True, True
// . . .
}
}
3.18 Matrix: inverse
„inverse matrix:“ een vierkante matrix in verband met een ander vierkante matrix van dezelfde grootte zodanig zijn dat het product van de twee matrices is de „identiteit matrix.“
Als een matrix heeft een inverse matrix verbonden, bestaat er slechts een inverse matrix.
Als de determinant van een matrix nul is, dan is de matrix is de naam „enkelvoud,“ en de matrix niet over een bijbehorende inverse matrix.
De onderstaande code berekent de inverse van een vierkante, niet-singuliere matrix.
Inversen voor 1*1, 2*2, 3*3, en 4*4 matrices zijn berekend door expliciete formules.
Met behulp van expliciete formules kan het sneller dan met behulp van een algemene procedure.
De expliciete formules zijn ook betrouwbaar.
De procedure voor het algemene geval gebruikt LU factoring, maar er zijn vele andere methoden van computer inversen van matrices.
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 Vermenigvuldigingsmethode product van een matrix en een vector
De vermenigvuldigingsmethode product van een matrix (M) en een vector (a) is een vector (b); (M)*(a) = (b).
Matrix (M) „transformeert“ vector (a) te vector (b).
De onderstaande code berekent de vermenigvuldigingsmethode product van een matrix en een vector.
public class MatrixF64
{
// . . .
public static VectorF64 operator *( MatrixF64 m, VectorF64 v )
{
if (null == m)
{
return (null); // Matrix is not specified.
}
if
(
(null == m.entries)
|| (m.totalRows <= 0)
|| (m.totalColumns <= 0)
)
{
return (null); // Matrix is empty.
}
if (null == v)
{
// Vector is not specified.
return (VectorF64.Zero( m.totalRows ));
}
if (v.Dimensions() != m.totalColumns)
{
// The number of columns of the matrix must be equal
// to the number of rows of the vector.
return (VectorF64.Zero( m.totalRows ));
}
// The result vector has the same number of rows as the matrix.
VectorF64 result = VectorF64.Zero( m.totalRows );
// Component (i) of the result vector is equal to the
// dot product of row (i) of the matrix with the vector.
for (int i = 0; i < m.totalRows; i++)
{
double dotProduct = 0.0;
for (int j = 0; j < m.totalColumns; j++)
{
dotProduct += (m[i, j] * v[j]);
}
result[i] = dotProduct;
}
return (result);
}
// . . .
public static void Test()
{
// . . .
// Examples of a product of a matrix and a vector:
MatrixF64 t3x3 =
new MatrixF64
(
3, 3,
2.0, 5.0, 7.0, // row 0
-4.0, -1.0, 6.0, // row 1
9.0, 8.0, 3.0 // row 2
);
VectorF64 p3x1 = new VectorF64( 1.0, 2.0, 3.0 );
VectorF64 q3x1 = t3x3 * p3x1;
q3x1.WriteLine();
// ( 33, 12, 34 )
// The multiplicative product of the inverse
// of the matrix and the result vector computed
// above should reproduce the original vector.
MatrixF64 inverset3x3 = t3x3.Inverse();
VectorF64 prodinvt3x3andq3x1 = inverset3x3 * q3x1;
prodinvt3x3andq3x1.WriteLine();
// ( 1, 2, 3 )
// This vector is equal to the original vector,
// as expected.
// . . .
}
}
3.20 Met behulp van homogene matrices te vertalen homogene vectoren
„vertaling:“ Een operatie T(s), gespecificeerd door een bepaalde vector (s), dat zet het vector (r) tot een overeenkomstige vector (r + s); T(s) * (r) = (r + s).
T(-s) is het omgekeerde van T(s); Inverse(T(s)) = T(-s); T(-s) * (r) = (r - s).
Een vertaling mag worden vertegenwoordigd door een matrix, met als gevolg dat vermenigvuldigen met de vertaling door een matrix vector is het effect van het vertalen van die vector.
De bijdragen van een vertaling operatie aan de componenten van de vector resultaat niet afhangen, direct of indirect, op de onderdelen van de oorspronkelijke vector.
Zo is een vertaling matrix moet worden gestructureerd op een wijze die haar in staat stelt bij te dragen een vertaling naar het resultaat zonder te worden beïnvloed, direct of indirect, door de elementen van de oorspronkelijke vector.
Een „homogene vector“ is een vector met een laatste element dat gelijk is aan 1.
Een „homogene matrix“ is een matrix met een rij met alle items gelijk aan 0 behalve voor een definitieve inschrijving van 1.
De vermenigvuldigingsmethode product van een homogene matrix en een andere homogene matrix is een homogene matrix.
De vermenigvuldigingsmethode product van een homogene matrix en een homogene vector is een homogeen vector.
Een interessant kenmerk van dit product is dat de vermeldingen in de laatste kolom van de matrix homogeen, met uitzondering van de laatste binnenkomst, worden toegevoegd direct naar de componenten van het resultaat vector, zonder te worden beïnvloed, direct of indirect, door de componenten van de oorspronkelijke vector.
Zo is een vertaling kan worden vertegenwoordigd door een homogene matrix, en die een matrix kan worden gebruikt voor het vertalen van homogene vectoren.
Als we willen vertegenwoordigen posities en richtingen in (d)-dimensionale ruimte, en we wensen te gebruiken matrices te vertalen posities, dan homogene matrices en homogene vectoren kan worden gebruikt.
Echter, de homogene vectoren moeten (d+1) onderdelen en de homogene matrices moet ((d+1)*(d+1)) inzendingen.
Het laatste onderdeel van een homogene vector en de laatste rij van een homogene matrix met geen enkele van de (d) coördinaten die wij willen vertegenwoordigen in (d)-dimensionale ruimte, maar in plaats daarvan dienen alleen mogelijk te maken krijgen matrices te doen sommige operaties onafhankelijk van de (d) coördineren waarden.
Met behulp van homogene matrices en homogene vectoren maakt het gemakkelijker om te samengestelde transformaties en transformaties op vectoren, maar ook andere activiteiten steeds meer lastig.
Als vertegenwoordiger van standpunten of richtingen in (d)-dimensionale ruimte met behulp van homogene vectoren betekent dat de vectoren zal hebben (d+1) onderdelen, met het laatste element dat gelijk is aan 1.
Computing de stip product van twee homogene vectoren met behulp van een regelmatige dot-product functie zou impliciet de stip product van de eerste (d) onderdelen met een waarde van niet-relevante 1 te wijten aan de uiteindelijke onderdelen van 1.
Daarom, om de dot product van de betrokken (d) onderdelen van de homogene vector (met (d+1) componenten), de stip product van homogene vectoren moet worden gedaan met behulp van een speciale „homogene vector dot product“ functie die alleen maar slaat het laatste onderdeel van de vectoren.
Een andere mogelijkheid is de stip product kan worden berekend met de reguliere dot product functie met dien verstande dat een bedrag van 1 moet worden afgetrokken van het resultaat voor de interpretatie van het resultaat als met betrekking tot de (d) relevante coördinaten.
De onderstaande code laat zien hoe een matrix kan worden gevormd die vervolgens kan worden gebruikt voor het vertalen van een vector.
public class MatrixF64
{
// . . .
public static MatrixF64 FormHomogeneousTranslationMatrix( VectorF64 v )
{
// This method assumes the specified vector is a homogeneous
// vector, where the final component is not relevant to the
// translation. Thus, we form a square homogeneous matrix
// with a total number of rows equal to the number of
// components of the specified vector. We form an identity
// matrix of the required size, and copy all but the final
// component of the specified vector to the final column of
// the matrix.
if (null == v)
{
return (MatrixF64.Zero( 1, 1 )); // Vector is not specified.
}
if (v.Dimensions() <= 0)
{
return (MatrixF64.Zero( 1, 1 )); // Vector is empty.
}
MatrixF64 result =
MatrixF64.Identity( v.Dimensions() );
// Add all but the final component of the specified vector
// to the final column of the result matrix. The final
// entry of the final column of the matrix will remain one (1).
for (int i = 0; i < (result.totalRows - 1); i++)
{
result[i, (result.totalColumns - 1)] = v[i];
}
return (result);
}
// . . .
public static void Test()
{
// . . .
// Example of using a homogeneous matrix to
// translate a homogeneous vector:
// Forming a homogeneous matrix that represents
// a translation:
VectorF64 homogeneousTranslationVector4x1 =
new VectorF64( 10.0, 20.0, 30.0, 1.0 );
MatrixF64 translation4x4 =
MatrixF64.FormHomogeneousTranslationMatrix
( homogeneousTranslationVector4x1 );
translation4x4.WriteLine();
// [ 1, 0, 0, 10 ]
// [ 0, 1, 0, 20 ]
// [ 0, 0, 1, 30 ]
// [ 0, 0, 0, 1 ]
// Using the homogeneous matrix to translate
// a homogeneous vector:
VectorF64 homogeneousVector =
new VectorF64( 5.0, 6.0, 7.0, 1.0 );
VectorF64 translatedVector =
translation4x4 * homogeneousVector;
translatedVector.WriteLine();
// ( 15, 26, 37, 1 )
// The relevant vector components (5,6,7)
// were translated by (10,20,30) by the
// homogeneous translation matrix.
// . . .
}
}
3.21 (d)-dimensionale ruimte: rotaties
„rotatie in (d)-dimensionale ruimte:“ een operatie (R) door een twee-dimensionale rotatie vliegtuig (S), een rotatiehoek (t), en een rotatie middelpunt (c), dat zet een punt (p) naar een overeenkomstige punt (q), zodanig zijn dat:
(1) Length( q - c ) = Length( p - c ); de afstand tussen de nieuwe en de rotatie middelpunt hetzelfde is als de afstand tussen de oorspronkelijke en de rotatie middelpunt.
Zo komt het nieuwe punt is op het oppervlak van een bepaalde (d)-dimensionale ruimte;
(2) Laat (P) een vliegtuig dat is parallel aan de rotatie vliegtuig (S) en bevat punt (p).
Punt (q) zal worden in het vliegtuig (P).
Zo komt het nieuwe punt is op een bepaalde twee-dimensionale vlak in (d)-dimensionale ruimte;
(3) Dot( (q - c), (p - c) ) = Length( q - c ) * Length( p - c ) * Cos( t ); de stip product van de vector uit de roulatie middelpunt van de oorspronkelijke letter en de vector uit de roulatie middelpunt van het nieuwe punt is gelijk aan het product van de lengte van deze twee vectoren en de cosinus van de rotatie hoek.
Zo komt het nieuwe punt is op het oppervlak van een bepaalde (d)-dimensionale kegel;
Dwang (1) beperkt het nieuwe punt op het oppervlak van een bepaalde (d)-dimensionale ruimte.
Voorwaarde (2) verdere beperkende het nieuwe punt aan een bepaalde twee-dimensionale vlak in (d)-dimensionale ruimte.
Daarom moet het nieuwe punt moet op een bepaalde twee-dimensionale cirkel in (d)-dimensionale ruimte.
Voorwaarde (3) verdere beperkende het nieuwe punt aan een bepaalde (d)-dimensionale kegel.
Daarom moet het nieuwe punt moet op een bepaald punt (voor Cos( t ) = 1 of Cos( t ) = (-1)) of bij een van de twee punten (voor (-1) < Cos( t ) < 1).
Een laatste beperking is te kiezen voor een van de twee mogelijke punten die het gevolg zijn van de hierboven vermelde beperkingen.
Laten T(s) is een „vertaling“ operatie, waar (s) is een willekeurige vector, zodanig dat T(s) * (r) = (r + s), en Inverse(T(s)) * (r) = (r - s).
Een rotatie R(S,t,c), vlak voor S, hoek t, en middelpunt c, kan worden uitgedrukt in termen van vertaling en een vruchtwisseling met een middelpunt in de oorsprong: R(S,t,c) = T(c) * R(S,t,0) * T(-c) = T(c) * R(S,t,0) * Inverse(T(c)).
Laten R(S,t) een rotatie operatie met een rotatie centrum punt op de oorsprong van (d)-dimensionale ruimte.
Zo heeft elke willekeurige rotatie kan worden uitgedrukt als: R(S,t,c) = T(c) * R(S,t) * T(-c) = T(c) * R(S,t) * Inverse(T(c)).
Laten we als twee-dimensionale rotatie vliegtuigen (S) dat ook de oorsprong van (d)-dimensionale ruimte en zijn evenwijdig aan de coördinatie van twee assen.
Laten (a) en (b) zijn getallen die een coördineren assen in (d)-dimensionale ruimte, zodanig dat 0 <= a <= (d-1), en 0 <= b <= (d-1), en (a != b).
We definiëren R(a,b,t) moet een rotatie om de oorsprong van (d)-dimensionale ruimte, met een rotatie vlak evenwijdig aan de coördinatie van assen aangegeven door coördineren indices (a) en (b), en met een hoek (t), zodanig dat R(a,b,(pi/2)) zou zetten een positieve waarde voor de coördinatie van de coördinatie as aangegeven door (a) te coördineren een positieve waarde voor de coördinatie van as aangegeven door (b); voorbeeld: R(a,b,(pi/2)) * e(a) = e(b), waar e(k) geeft aan dat er een eenheid basis vector die overeenkomt met de coördinatie as (k).
Deze definitie biedt de definitieve disambiguation nodig zijn voor de volledige bepaling van het nieuwe punt geproduceerd door een rotatie.
Een rotatie wat beweegt punten langs de as a richting as b volgens een hoek t gelijk is aan de draaiing die beweegt punten langs de as b richting as a door een hoek (-t), dat wil zeggen, R(a,b,t) = R(b,a,(-t)).
Als (d >= 2), dan (d)-dimensionale ruimte heeft ((d*(d-1))/2) paren van coördineren assen.
Dus, voor (d >= 2), (d)-dimensionale ruimte heeft ((d*(d-1))/2) verschillende as-rotatie gebracht vliegtuigen.
(Voorbeelden: d=2: () x-y vliegtuig; d=3: (x-y vliegtuig, y-z vliegtuig, z-x vliegtuig); d=4: (x-y vliegtuig, y-z vliegtuig, z-x vliegtuig, x-w vliegtuig, z-w vliegtuig, w-y vliegtuig).)
De onderstaande code laat zien hoe een matrix kan worden gebruikt om de rotatie R(a,b,t).
Een voorbeeld is gebleken dat gebruik maakt van een rotatie-matrix in een homogene wijze te draaien een homogene vector naar een nieuwe homogene vector.
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 Voorgestelde definitie van „tegenwijzerzin“
In twee-dimensionale ruimte, de kwartaal-rotatie beurt dat het converteren van de vector e(0) te e(1) wordt beschouwd tegen.
De rotatie is: R(0,1,(pi/2)) * e(0) = e(1).
Daarom is de rotatie R(0,1,t) is tegen.
In de driedimensionale ruimte, het kwartaal-rotatie beurt dat het converteren van het vector e(0) te e(1) wordt beschouwd tegen.
De kwartaal-rotatie beurt dat het converteren van de vector e(1) te e(2) wordt ook geacht tegen.
De kwartaal-rotatie beurt dat het converteren van de vector e(2) te e(0) wordt ook geacht tegen.
Daarom is de rotaties { R(0,1,t), R(1,2,t), R(2,0,t) } zijn tegen.
De Booleaanse formule past in de definitie van tegen de twee-dimensionale en drie-dimensionale ruimte, en extrapolates het patroon voor het algemene geval van de rotatie R(a,b,t):
bool counterClockwise =
(
( (a < b) && (((a + b) % 2) == 1) )
|| ( (a > b) && (((a + b) % 2) == 0) )
);
Deze regel zou de rotaties als „tegen de“ hieronder:
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 ruimte: de klok mee rotaties
Met behulp van de voorgestelde definitie van „tegen de“ in het vorige deel van de rotaties daaronder zijn „de klok mee:“
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),
...
Bijvoorbeeld, in vier-dimensionale ruimte, de zes rotatie matrices hieronder geven de wijzers van de klok rotaties:
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 ruimte: totaal aantal as-gebonden oriëntaties
Het aantal afzonderlijke as-gebonden oriëntaties in (d)-dimensionale ruimte is:
distinctOrientations = Factorial( d ) * IntPow( 2, (d-1) );
Een argument ter ondersteuning van deze formule wordt weergegeven in onderstaande paragrafen.
Experimenteel bewijs na te gaan of de formule is ook hieronder weergegeven.
Laten { g(0), g(1), ..., g(d-1) } worden een reeks van (d) onderling loodrecht-eenheid vectoren die aan een object in (d)-dimensionale ruimte.
Deze eenheid vectoren geven aan de oriëntatie van het object in verhouding tot het coördineren van de assen (d)-dimensionale ruimte.
Laat de „standaard oriëntatie“ van het object moet zodanig zijn dat g(i) = e(i), voor alle 0 <= i <= (d-1); elke geaardheid vector van het object wordt in verband gebracht met een aparte eenheid basis vector van een coördinatenstelsel in de (d)-dimensionale ruimte.
De eerste oriëntatie vector van het object kunnen worden parallel met een van de (d) assen van de (d)-dimensionale ruimte, en kunnen in beide richtingen ten opzichte van dat as.
Zo is de eerste oriëntatie vector heeft (2 * d) opties.
De tweede oriëntatie vector van het object kunnen worden parallel aan een van de (d-1) overige assen, met twee mogelijke richtingen ten opzichte van dat as.
Daarom werd in de tweede oriëntatie vector heeft (2 * (d-1)) opties.
Op dezelfde wijze, de volgende geaardheid vector heeft (2 * (d-2)) opties.
De definitieve geaardheid vector is gedwongen om de laatste resterende as, en is gedwongen om een van de twee mogelijke richtingen langs deze as (aan het behoud van de pariteit van de set van oriëntatie-vectoren; de determinant van een vastgestelde volgorde van de vectoren geaardheid moet constant ).
Zo heeft de definitieve geaardheid vector alleen heeft (1) opti