Linear Algebra
Colin Fahey
1. Software
LinearAlgebra.zip
Lineær algebra kildekoden (C#)
19910 bytes
MD5: 11d8c8035cac30ba543e5e0b72ee9767
2. Indledning
Denne artikel beskriver, vektorer og matricer i (d)-dimensionelle rum.
3. (d)-dimensionelle rum: Attributter
3.1 Array
"array:" En samling af variabler, sådan at hver variabel har et unikt navn, og sådan, at navne kan blive tildelt en ordre.
Integer værdier kan bruges som navne til variabler i en array.
For eksempel, hvis en række indeholder (d) variabler, så integer værdier { 0, 1, 2, ..., (d-1) } kan være navne tildeles de variabler i array.
3.2 (d)-dimensional vektor
"(d)-dimensional vektor:" En vifte af (d) variabler.
"vektor komponent:" En variabel inden for en vektor.
3.3 (d)-dimensionelle vektorrum
"en-dimensionelle rum:" Det komplette sæt af værdier, der kan gemmes af en variabel.
"(d)-dimensionelle rum:" Det komplette sæt af kombinationer af værdier, der kan gemmes af en vifte af (d) variabler.
Den formelle definition på en "vektor rummet:"
Lad (T) være en grundlæggende type (eksempler: reelt tal, heltal, komplekse tal, rationelt tal, osv.).
Enhver variabel for en grundlæggende type kaldes en "skalar."
En "(T)-typen (d)-dimensional vektorrum" er det sæt (S) af (d)-dimensionelle vektorer have to operationer, vektor Ud (+), og skalar multiplikation (*), der opfylder de betingelser nedenfor:
(1) Hvis (v) og (w) er to vektorer i (S), så (v + w) er også en vektor i (S);
(2) Hvis (u), (v), og (w) er alle tre vektorer i (S), så (u + v) + w = u + (v + w);
[tilsætningsstof kommutativitet]
(3) Hvis (v) og (w) er to vektorer i (S), så (v + w) = (w + v);
[tilsætningsstof Associativitet]
(4) Der er en "nul vektor," (0) i (S), således at der for enhver vektor (v) i (S), (v + (0)) = v;
[tilsætningsstof identitet]
(5) Hvis (c) er enhver skalar af type (T), og (v) er enhver vektor i (S), så produktet (c * v) er en vektor i (S);
(6) Hvis (a), (b), og (c) er nogen scalars af type (T), og (v) og (w) er nogen vektorer i (S), så (a + b) * v = a*v + b*v, og c*(v + w) = c*v + c*w;
[multiplicative distributivity]
(7) Hvis (a) og (b) er nogen scalars af type (T), og (v) er enhver vektor i (S), så (a*b)*v = a*(b*v);
(8) Hvis "1" er en skalar af type (T) sådan, at (1*1)=1, og (v) er enhver vektor i (S), så (1*v) = v;
(9) For hver vektor (v) i (S), vektor (-1)*v = -v opfylder v + (-v) = (0);
[tilsætningsstof inverse]
3.4 (d)-dimensional vektor-kode
Nedenstående kode viser, hvordan en (d)-dimensional vektor, med 64-bit floating-point komponenter, kan gennemføres.
En matrix er tilstrækkeligt til at repræsentere en vektor.
Nedenstående kode har et array i en klasse, kun for bekvemmelighed.
Koden er ikke beregnet til at blive effektiv.
En struktur (eksempel: "struct" en værdi type) repræsenterer vektorer af et fast antal dimensioner (eksempler: 3 eller 4) må forventes at være langt mere effektiv end den generelle (d)-dimensional vektor klasse vist her.
Selv om koden nedenfor definerer en vektor med floating-point komponenter, dette dokument også gør brug af vektorer med heltal komponenter.
Nedenstående kode kan nemt blive ændret for at gennemføre vektorer med heltal komponenter.
using System;
using System.Collections.Generic;
using System.Text;
// Multidimensional vector with 64-bit floating-point components:
public class VectorF64
{
private double[] components;
public int Dimensions()
{
if (null == this.components)
{
return (0);
}
return (this.components.Length);
}
public VectorF64()
{
this.components = null;
}
public VectorF64( params double[] paramValues )
{
this.components = null;
if (null == paramValues)
{
return;
}
int dimensions = paramValues.Length;
this.components = new double[dimensions];
for (int i = 0; i < dimensions; i++)
{
this.components[i] = paramValues[i];
}
}
public VectorF64( VectorF64 other )
{
this.components = null;
if (null == other)
{
return;
}
if (null == other.components)
{
return;
}
int dimensions = other.Dimensions();
this.components = new double[dimensions];
for (int i = 0; i < dimensions; i++)
{
this.components[i] = other.components[i];
}
}
public double this[int index]
{
get
{
if (null == this.components)
{
return (0.0);
}
if
(
(index >= 0)
&& (index < this.components.Length)
)
{
return (this.components[index]);
}
return (0.0);
}
set
{
if (null == this.components)
{
return;
}
if
(
(index >= 0)
&& (index < this.components.Length)
)
{
this.components[index] = value;
}
}
}
public void Write( int precision )
{
if (null == this.components)
{
Log.Write( String.Empty + '(' + ' ' + ')' );
return;
}
int dimensions = this.Dimensions();
if (0 == dimensions)
{
Log.Write( String.Empty + '(' + ' ' + ')' );
return;
}
// Determine the largest component width in characters
// so that we can make all components an equal width.
int largestComponentWidth = 1;
for (int i = 0; i < dimensions; i++)
{
// { index [,minwidth] [:typeCode[precision]] }
// (minwidth<0) means left-justify
String text = String.Format( String.Empty + '{' + '0' + ':'
+ 'g' + precision + '}', this[i] );
if (text.Length > largestComponentWidth)
{
largestComponentWidth = text.Length;
}
}
Log.Write( '(' );
for (int i = 0; i < dimensions; i++)
{
Log.Write( ' ' );
String text =
String.Format( String.Empty + '{' + '0' + ','
+ largestComponentWidth + ':' + 'g' + precision + '}',
this[i] );
Log.Write( text );
if ((i + 1) < dimensions)
{
Log.Write( ',' );
}
else
{
Log.Write( ' ' );
}
}
Log.Write( ')' );
}
public void WriteLine( int precision )
{
this.Write( precision );
Log.WriteLine();
}
public void WriteLine()
{
const int defaultPrecision = 8;
this.Write( defaultPrecision );
Log.WriteLine();
}
// . . .
public static void Test()
{
// A 3-dimensional vector with 64-bit floating-point components:
VectorF64 v3 = new VectorF64( 0.0, 1.0, 2.0 );
v3.WriteLine(); // ( 0, 1, 2 )
// A 4-dimensional vector with 64-bit floating-point components:
VectorF64 v4 = new VectorF64( 0.0, 1.0, 2.0, 3.0 );
v4.WriteLine(); // ( 0, 1, 2, 3 )
// . . .
}
}
"(d)-dimensional nul vektor:" En vektor med alle (d) komponenter er lig med 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)-dimensionelle rum: punkt
"(d)-dimensionelle rum punkt:" En vifte af (d) variabler med særlige værdier "(koordinere værdier)."
"(d)-dimensionelle rum oprindelse:" En vifte af (d) variabler med alle værdier lig med nul.
3.6 (d)-dimensionelle vektorer: ikke-relative og relative
En (d)-dimensional vektor, som er "ikke-relative" er en (d)-dimensional vektor, som direkte repræsenterer en status eller konfiguration.
En (d)-dimensional vektor, der er "relativ" er en (d)-dimensional vektor, der repræsenterer ændringer i et sæt af komponenter.
En relativ vektor kan udgøre forskellen mellem to ikke-relativ vektorer.
Da en relativ vektor, fastsættelse af en ikke-relativ status eller konfiguration benytter denne relative vektor kræver kombinere at den relative vektor med en ikke-relativ vektor.
Ikke-relative vektorer og relative vektorer er begge vektorer.
Hvorvidt en bestemt vektor er ikke-relativ eller relative skal specificeres, hvor vektoren er defineret.
Hvis en (d)-dimensional vektor tolkes som værende ikke-relativ, så (d)-dimensional vektor kan repræsentere et punkt i (d)-dimensionelle rum.
3.7 (d)-dimensional vektor addition, subtraktion, og skalering
Vector addition, subtraktion, og skalering:
public class VectorF64
{
// . . .
public static VectorF64 operator +( VectorF64 a, VectorF64 b )
{
if ((null == a) || (null == b))
{
return (new VectorF64()); // Vector not specified.
}
if ((null == a.components) || (null == b.components))
{
return (new VectorF64()); // Vector is empty.
}
if (a.Dimensions() != b.Dimensions())
{
return (new VectorF64()); // Vectors not the same size.
}
int dimensions = a.Dimensions();
VectorF64 result = VectorF64.Zero( dimensions );
for (int i = 0; i < dimensions; i++)
{
result[i] = a[i] + b[i];
}
return (result);
}
public static VectorF64 operator -( VectorF64 a, VectorF64 b )
{
if ((null == a) || (null == b))
{
return (new VectorF64()); // Vector not specified.
}
if ((null == a.components) || (null == b.components))
{
return (new VectorF64()); // Vector is empty.
}
if (a.Dimensions() != b.Dimensions())
{
return (new VectorF64()); // Vectors not the same size.
}
int dimensions = a.Dimensions();
VectorF64 result = VectorF64.Zero( dimensions );
for (int i = 0; i < dimensions; i++)
{
result[i] = a[i] - b[i];
}
return (result);
}
public static VectorF64 operator -( VectorF64 a )
{
if (null == a)
{
return (new VectorF64()); // Vector not specified.
}
if (null == a.components)
{
return (new VectorF64()); // Vector is empty.
}
int dimensions = a.Dimensions();
VectorF64 result = VectorF64.Zero( dimensions );
for (int i = 0; i < dimensions; i++)
{
result[i] = (-( a[i] ));
}
return (result);
}
public static VectorF64 operator *( double scale, VectorF64 a )
{
if (null == a)
{
return (new VectorF64()); // Vector not specified.
}
if (null == a.components)
{
return (new VectorF64()); // Vector is empty.
}
int dimensions = a.Dimensions();
VectorF64 result = VectorF64.Zero( dimensions );
for (int i = 0; i < dimensions; i++)
{
result[i] = scale * a[i];
}
return (result);
}
public static VectorF64 operator *( VectorF64 a, double scale )
{
if (null == a)
{
return (new VectorF64()); // Vector not specified.
}
if (null == a.components)
{
return (new VectorF64()); // Vector is empty.
}
int dimensions = a.Dimensions();
VectorF64 result = VectorF64.Zero( dimensions );
for (int i = 0; i < dimensions; i++)
{
result[i] = scale * a[i];
}
return (result);
}
// . . .
public static void Test()
{
// . . .
// Examples of vector addition, subtraction, and scaling:
VectorF64 a = new VectorF64( 0.0, 1.0, 2.0, 3.0 );
a.WriteLine(); // ( 0, 1, 2, 3 )
VectorF64 b = new VectorF64( 3.0, 2.0, 1.0, 0.0 );
b.WriteLine(); // ( 3, 2, 1, 0 )
VectorF64 c = new VectorF64();
c.WriteLine(); // ( )
c = a + b;
c.WriteLine(); // ( 3, 3, 3, 3 )
c = a - b;
c.WriteLine(); // ( -3, -1, 1, 3 )
c = -b;
c.WriteLine(); // ( -3, -2, -1, 0 )
c = 3.0 * a;
c.WriteLine(); // ( 0, 3, 6, 9 )
// . . .
}
}
3.8 (d)-dimensional grundlag vektorer
Den formelle definition på "grundlag" af et vektorrum:
Lad (T) være en grundlæggende type (eksempler: reelt tal, heltal, komplekse tal, rationelt tal, osv.).
Enhver variabel for en grundlæggende type kaldes en "skalar."
Lad (V) være en "(T)-typen (d)-dimensional vektorrum."
Hvis den ikke-nul vektorer { u1, u2, ..., ud } i (V) er sådan, at enhver vektor (v) i (V) kan skrives som en "lineær kombination" af disse vektorer, v = c1*u1 + c2*u2 + ... + cd*ud, hvor { c1, c2, ..., cd } er scalars af type (T), så (V) er "kalibreret" fra vektorer { u1, u2, ..., ud }.
Et sæt af ikke-nul vektorer { u1, u2, ..., ud } at "span" vektorrum (V) er navngivet på "basis" af (V).
Et simpelt "grundlag" af en (d)-dimensional vektorrum er det sæt af (d) adskiller (d)-dimensionelle vektorer, hver med en komponent svarende til en og alle andre komponenter er lig med nul.
Sådanne grundlag vektorerne er "Orthonormal," hvilket betyder, at de er gensidigt vinkelrette "(retvinklede)," og at hver vektor har enheden længde.
Hver sådan vektor er en enhed vektor parallel til en af de (d) koordinere akserne.
Tilkendegive vilkårlig vektor som en "lineær kombination" af disse grundlag vektorer direkte; hver bestanddel af den vilkårlige vektoren er ganget med den tilsvarende grundlag vektor, og disse produkter er føjet sammen om at danne vilkårlig vektor.
Nedenstående kode viser, hvordan en vektor kan udtrykkes som en "lineær kombination" af basis vektorer.
Nedenstående kode vilkårligt definerer en Orthonormal sæt grundlag vektorer.
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 )
// . . .
}
}
Enhver vektor i (d)-dimensionelle rum kan udtrykkes som en sum af produkterne af tal og grundlag vektorer:
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)-dimensionelle rum: Afstanden mellem punkterne
Lad (P) være en (d)-dimensional vektor, der repræsenterer et punkt i (d)-dimensionelle rum.
Lad (Q) være en (d)-dimensional vektor, der repræsenterer et punkt i (d)-dimensionelle rum.
Lad (R) være en (d)-dimensional vektor, der repræsenterer ændringen i (d) koordinater til at komme fra punkt (P) til punkt (Q); R = (Q - P).
I én-dimensionelle rum, P = ( p0 ), Q = ( q0 ), og R = (Q - P) = ( q0 - p0 ).
Afstanden mellem de to punkter er: Abs( q0 - p0 ).
I to-dimensionelle rum, P = ( p0, p1 ), Q = ( q0, q1 ), og R = (Q - P) = ( q0-p0, q1-p1 ).
Tolkning de to vinkelrette fordrivelser som den vinkelrette sider af en højre-trekant, hvor afstanden mellem de punkter svarer til længden af hypotenusen af denne trekant.
De Pythagorean formel, (a*a) + (b*b) = (c*c), hvor (a) og (b) er længderne af den vinkelrette sider af en højre-trekant, og (c) er længden af hypotenusen (påvirke side), kan bruges til at bestemme afstanden mellem de to punkter: Sqrt( Sq(q0-p0) + Sq(q1-p1) ).
I tre-dimensionelle rum, P = ( p0, p1, p2 ), Q = ( q0, q1, q2 ), og R = (Q - P) = ( q0-p0, q1-p1, q2-p2 ).
Tolke forskydning ( q0-p0, q1-p1, 0 ) som den vinkelrette sider af en højre-trekant og ved at bruge Pythagorean formel, giver afstanden mellem punkt (P) og det punkt ( q0, q1, p2 ): d01 = Sqrt( Sq(q0-p0) + Sq(q1-p1) ).
Denne forskydning ( q0-p0, q1-p1, 0 ) er vinkelret på den forskydning ( 0, 0, q2-p2 ), og en anden højre-trekant kan dannes, og Pythagorean formel kan bruges igen.
Således er afstanden fra punkt (P) til punkt (Q) er givet ved: 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) ).
Den metode for at udvide afstanden formlen fra to-dimensionelle rum til tre-dimensionelle rum kan anvendes flere gange i sidste ende at bestemme afstanden formlen for (d)-dimensionelle rum: Sqrt( (Sq(q0-p0) + Sq(q1-p1)) + Sq(q2-p2) + ... + Sq(qd-pd) ).
Nedenstående kode definerer en funktion kaldet "længde," der beregner længden af en (d)-dimensional vektor.
Når en vektor repræsenterer forskydning mellem to punkter i (d)-dimensionelle rum, længden af vektoren, udgør afstanden mellem disse to punkter.
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)-dimensionelle vektorer: dot produkt
De "dot produkt" konverterer to (d)-dimensionelle vektorer til et nummer.
Nedenstående kode beregner en prik produkt af to vektorer:
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);
}
// . . .
}
For enhver vektor, (A), Length(A) = Sqrt(Dot(A,A)).
3.11 (d)-dimensionelle vektorer: definition af "parallelle"
Vektorer (A) og (B) er "parallelle" hvis alle de følgende er sandt:
(1) A.Length() > 0;
(2) B.Length() > 0;
(3) Abs(Dot(A,B)) = A.Length()*B.Length().
Nedenstående kode afgør om et par af vektorer er parallelle (muligvis anti-tilpassede).
Floating-point tal kan ophobes fejl i fraktioneret dels skyldes den begrænsede præcision, og derfor, edb-kode bør omfatte ikke-nul tolerance når man sammenligner floating-point tal.
Koden omfatter eksempelvis tolerance værdier, men eksemplet tolerance værdier måske ikke være hensigtsmæssigt for nogle opgaver.
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)-dimensionelle vektorer: definition af "vinkelret"
Vektorer (A) og (B) er "vinkelret" hvis alle de følgende er sandt:
(1) A.Length() > 0;
(2) B.Length() > 0;
(3) Abs(Dot(A,B)) = 0.
Nedenstående kode afgør om et par af vektorer er vinkelret. Floating-point tal kan ophobes fejl i fraktioneret dels skyldes den begrænsede præcision, og derfor, edb-kode bør omfatte ikke-nul tolerance når man sammenligner floating-point tal.
Koden omfatter eksempelvis tolerance værdier, men eksemplet tolerance værdier måske ikke være hensigtsmæssigt for nogle opgaver.
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 Matricer
"Matrix:" En samling af variabler, sådan at hver variabel har en unik kombination af en "række" navn og en "kolonne" navn.
"indrejse:" En variabel inden for en matrix.
Integer værdier kan bruges som "rækken" navne og "kolonne" navne til variabler i en matrix.
For eksempel, hvis en matrix har (totalRows) rækker og (totalColumns) kolonner, så integer værdier { 0, 1, ..., (totalRows-1) } kan være navne tildeles de rækker, og integer værdier { 0, 1, ..., (totalColumns-1) } kan være navne tildeles de kolonner.
Således er en variabler i en matrix kan angives ved at angive et par af heltal, ( row, column ), hvilket tyder på en kombination af række og kolonne indekser svarer til de specifikke variabel.
Størrelsen af en matrix er angivet som "(totalRows) * (totalColumns)" (eller "(totalRows) ved (totalColumns))."
Denne kendelse af dimensioner er den samme som rækkefølgen af dimensioner bruges til at angive poster i matrix "(( row, column ))."
[Denne konvention er lidt uheldigt, fordi for mange to-dimensional anvendelser (eksempler: billeder, grafer osv.) den sædvanlige konvention er at angive dimensioner som "width * height" og koordinerer som "( horizontal, vertical )" (eller "( x, y ))."
Dette er det modsatte af rækkefølgen af dimensioner og koordinater bruges til at beskrive matricer og deres poster. ]
En matrix med (totalRows) lig med (totalColumns) er opkaldt "square" ellers, matrix er opkaldt "rektangulære."
En matrix kan betragtes som indeholdende et sæt af "rækken vektorer," hvor variabler i hver række er fortolkes som tilhørende en vektor.
En matrix kan også betragtes som indeholdende et sæt "kolonne vektorer," hvor variabler i hver kolonne er fortolkes som tilhørende en vektor.
Matricer kan repræsentere en bred vifte af matematiske relationer.
Betydningen af en matrix, og de operationer, som kan være relevante for behandling poster i en matrix, afhænger af konteksten.
Der er dog grundlæggende regler i matrix regning, der er relevante i mange sammenhænge, og disse grundlæggende regler vil blive fastlagt i et senere afsnit.
En matrix og en (totalColumns) værdi er tilstrækkelige til at repræsentere en matrix.
De array kan have (totalRows * totalColumns) variabler, og optagelse på ( row, column ) kan svare til array variabel i indekset ((totalColumns * row) + column).
Nedenstående kode definerer en matrix, med 64-bit floating-point poster.
En matrix og en (totalColumns) værdi er tilstrækkelige til at repræsentere en matrix.
Nedenstående kode er en matrix og en (totalColumns) værdi er indeholdt i en klasse, kun for bekvemmelighed.
Nedenstående kode er ikke beregnet til at blive effektiv.
En struktur (eksempel: "struct" en værdi type) repræsenterer matricer af særlig faste dimensioner (eksempler: 2*2, 3*3, eller 4*4) må forventes at være langt mere effektiv end den generelle (totalRows * totalColumns) klassen vist her.
Selv om koden nedenfor definerer en matrix med floating-point poster, giver denne artikel også gør brug af matricer med heltal poster.
Nedenstående kode kan let modificeret til at gennemføre matricer med heltal poster.
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:" En matrix med alle indgange lig med 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 addition, subtraktion, og multiplikation
Matrix addition, subtraktion, og multiplikation.
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 ]
// . . .
}
}
"identitet matrix:" En kvadratisk matrix (total rækker lig samlede kolonner) med angivelser svarende til (1) om diagonal (poster med en række indeks lige til sin kolonne indekset), og med alle andre poster svarer til (0).
Hvis en kvadratisk matrix, (M), er ganget med en "identitet matrix," (I), af det samme antal rækker (og søjler), produktet er lig med (M).
Mangedobling (I) ved (M) også producerer et produkt, svarende til (M).
Således "identitet matrix" er analog til antallet "1" for en mangedobling af numre (scalars).
Nedenstående kode skaber en identitet matrix med et bestemt antal rækker.
public class MatrixF64
{
// . . .
public static MatrixF64 Identity( int rows )
{
MatrixF64 identity = MatrixF64.Zero( rows, rows );
for (int i = 0; i < rows; i++)
{
identity[i, i] = 1.0;
}
return (identity);
}
// . . .
public static void Test()
{
// . . .
// Examples of multiplying by an identity matrix:
MatrixF64 identity3x3 = MatrixF64.Identity( 3 );
// [ 1, 0, 0 ]
// [ 0, 1, 0 ]
// [ 0, 0, 1 ]
MatrixF64 s3x3 =
new MatrixF64
(
3, 3,
0.0, 1.0, 2.0, // row 0
3.0, 4.0, 5.0, // row 1
6.0, 7.0, 8.0 // row 2
);
result = s3x3 * identity3x3;
result.WriteLine();
// [ 0, 1, 2 ]
// [ 3, 4, 5 ]
// [ 6, 7, 8 ]
// (the same as s3x3)
result = identity3x3 * s3x3;
result.WriteLine();
// [ 0, 1, 2 ]
// [ 3, 4, 5 ]
// [ 6, 7, 8 ]
// (the same as s3x3)
// . . .
}
}
3.15 Matrix: LU factoring back-substitution
"Matrix LU factoring:" En procedure, der konverterer en kvadratisk matrix, (M), at to nye kvadratiske matricer, (L) og (U), der har samme størrelse som (M), således at (L) * (U) = (M), og sådan, at matrix (U) er "øverste trekantede" (alle poster under diagonal åre nul) og matrix (L) er "lavere trekantede" (alle poster over diagonal er nul).
Matrix LU factoring kan bruges som del af en større procedure til at løse et system af ligninger, eller at finde den inverse af en matrix, eller at finde den afgørende faktor for en matrix.
Matrix LU factoring er et alternativ til "Gauss-Jordan afskaffelse procedure.
Gauss-Jordan udrydde den kræver et system af ligninger (A)*(x)=(b), mens LU factoring kun kræver en matrix (A).
Også efter fastsættelsen af LU factoring af en matrix (A), er det meget let at afgøre (x) givet nogen (b).
Proceduren for at løse med henblik vektoren (x) i (A)*(x)=(b), da (b) og LU factoring (L)*(U)=(A), omfatter en procedure, der vil blive navngivet "back-substitution" her.
Nedenstående kode indeholder en procedure for beregning af LU faktorer af enhver firkantede, ikke-ental matrix.
Nedenstående kode indeholder en procedure for at gøre LU back-substitution.
Forsigtig: Computeren koden nedenfor beregner den L og U faktorer i træk-permuted version af en bestemt matrix.
Også den L og U matrix resultater er samlet i en enkelt 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:" En funktion, der konverterer enhver kvadratisk (n * n) matrix (A) til et nummer, det(A), sådan at:
(1) Hvis en matrix (B) resultater fra udveksle to rækker eller to kolonner, med en matrix (A), så det(B) = (-(det(A)));
(2) Hvis en matrix (B) resultater fra multiplicere enhver række, eller en kolonne, i en matrix (A), som en række (c), så det(B) = c * det(A);
(3) Hvis en matrix (B) resultater fra tilføje et multiplum af en række til en anden række, eller fra at tilføje et multiplum af en kolonne til en anden kolonne og derefter det(B) = det(A).
(4) Hvis (A) er en øvre-triangulære matrix (A[i,j]=0 for alle (i>j)) eller en lavere trekantede matrix (A[i,j]=0 for alle (i<j)), derefter det(A) = (A[0,0] * A[1,1] * ... * A[n-1,n-1]);
(For eksempel er det afgørende for en identitet matrix er én; det(I) = 1.)
Regler (1), (2), og (3), kan bruges i en proces kaldet "Gaussian elimination, til at konvertere en kvadratisk matrix til en trekantet matrix.
Regel (4) kan bruges til at beregne den afgørende faktor for et trekantet matrix.
Den formelle definition på en "determinant:"
Lad M[n](K) betegne det sæt af alle (n * n) matricer over området (K).
Den "afgørende faktor" er den unikke funktion (F) sådan, at F:M[n](K) --> K med to attributter:
(1) F er skiftevis multilinear med hensyn til spalter (eller rækker), og
(2) F( I ) = 1.
En "multilinear" funktion (D) af en (n * n) matrix (A) kan skrives som: 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]] ) ), hvor summen er taget over alle (n^n) kombinationer af 0 <= k[i] <= (n-1), og hvor e[j] repræsenterer rækken (j) om identiteten matrix.
En sådan funktion kan være kendetegnet ved værdierne af D( e[k[0]], e[k[1]], ..., e[k[n-1]] ).
En "skiftevis multilinear" funktion er en multilinear funktion, (D), således at D( ..., e[i], ..., e[i], ... ) = 0, og D( ..., e[i], e[j], ... ) = (-( D( ..., e[j], e[i], ... ) )).
Således swapping spalter (eller rækker) ændringer tegnet af funktion.
Også udtrykket D( e[k[0]], ..., e[k[n-1]] ) er nul, hvis kombinationen af k[i] værdier ikke involvere alle unikke værdier, og hvis det sæt af værdier er ikke en permutation af tallene {0,1,...,(n-1)}.
Udlejning D( I ) = 1 (D( e[0], e[1], ..., e[n-1] ) = 1), i kombination med skiftevis multilinear attributter er beskrevet ovenfor, betyder, at alle værdier af de skiftevis multilinear funktion D( e[k[0]], e[k[1]], ..., e[k[n-1]] ) kan bestemmes.
Funktionen vil være ikke-nul, hvis det sæt k[i] værdier er en permutation af {0,1,...,(n-1)}, og vil have en størrelsesorden af en, hvis det ikke er nul, og tegnet vil være antallet af swaps forpligtet til at konvertere permutation til identitet permutation.
Værdien af en faktor for en (n * n) matrix (A) kan beregnes ved at addere produkter i form (sgn(p) * A[0,p[0]] * A[1,p[1]] * ... * A[(n-1),p[n-1]]) for hver enkelt permutation (p) af tallene {0,1,2,...,(n-1)}.
Udtrykket (sgn(p)) betegner "underskrift" (eller swap count) permutation (p), hvor (sgn(p)) er lig med (+1) hvis (p) er "en" lige "permutation," og er lig med (-1) hvis (p) er en "ulige permutation."
Da der findes (n!) permutationer af tallene {0,1,2,...,(n-1)}, denne formel har (n!) summands.
Andre procedurer for beregning af determinant (eksempler: der involverer Gaussian afskaffelse eller LU factoring, eller udvidelse af mindreårige osv.) implicit beregne et hierarki af mellemprodukter mængder, at de procedurer for at beregne den afgørende faktor i en kendelse af (n^3) trin.
Hvis den afgørende faktor for en matrix er lig nul, matrix ikke har et omvendt.
Hvis den afgørende faktor for en matrix er ikke-nul, matrix har et omvendt.
Hvis en matrix repræsenterer de koefficienter af en lineær system af ligninger, og den afgørende faktor for, at matrix er lig nul, så systemet af ligninger ikke har en entydig løsning.
Hvis en matrix repræsenterer de koefficienter af en lineær system af ligninger, og den afgørende faktor for, at matrix er ikke-nul, så systemet af ligninger har en unik løsning.
Nedenstående kode omfatter procedurer for beregning af faktor for ethvert kvadratisk matrix.
Afgørende faktorer for 1*1, 2*2, 3*3, og 4*4 matricer er beregnet ved udtrykkelig formler.
Brug af eksplicitte formler til at beregne determinanter for små matricer er sandsynligt at beregne resultater hurtigere end at bruge en almindelig procedure for at beregne disse determinanter.
Proceduren for det almindelige tilfælde bruger LU factoring, men der er mange andre metoder til databehandling determinanter.
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)-dimensionelle vektorer: tværs produkt
"tværs produkt:" en funktion i forbindelse med en (d)-dimensionelle rum, der konverterer en bestilt liste over (d-1) (d)-dimensionelle vektorer, { A(0), A(1), ..., A(d-2) }, til en (d)-dimensional vektor, sådan at den funktion har de attributter nedenfor:
(1) Hvis A(i)=0 (nul vektor), for enhver 0 <= i <= (d-2), så Cross( A(0), A(1), ..., A(d-2) )=0 (nul vektor).
Korset produkt af en bestilt liste af vektorer er nul, hvis enhver vektor i denne liste er nul;
(2) Hvis et par af vektorer A(i) og A(j) for nogen (i!=j) med 0 <= i <= (d-2) og 0 <= j <= (d-2), er parallelle (Abs(Dot(A(i),A(j))) = Length(A(i)) * Length(A(j))), så Cross( A(0), A(1), ..., A(d-2) )=0 (nul vektor).
Korset produkt af en bestilt liste af vektorer er nul, hvis to vektorer på denne liste er parallelle;
(3) Hvis alle vektor A(i) er ikke-nul, og hvis A(i) ikke er parallel med A(j) for alle (i!=j) med 0 <= i <= (d-2) og 0 <= j <= (d-2), så B = Cross( A(0), A(1), ..., A(d-2) ) er sådan, at Dot(A(i),B)=0 for alle 0 <= i <= (d-2).
Hvis tværs produkt af en bestilt liste af vektorer er ikke-nul, derefter på tværs produkt Resultatet er vinkelret på enhver vektor i bestilt liste af vektorer;
(4) Hvis (c) er en numerisk værdi, så 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) ).
Da en bestilt liste af vektorer på tværs produkt af den bestilte liste af vektorer med nogen en vektor multipliceret med en specifik numerisk værdi svarer til at multiplicere tværs produkt af den bestilte liste af vektorer ved den numeriske værdi;
(5) Da to numre (i) og (j), hvor (i!=j) og 0 <= i <= (d-2) og 0 <= j <= (d-2), og en ordentlig liste over (d-1) (d)-dimensionelle vektorer, { A(0), A(1), ..., A(d-2) }, og en anden bestilt liste over (d-1) (d)-dimensionelle vektorer, { B(0), B(1), ..., B(d-2) }, således at B(k) = A(k) for alle ((k != i) && (k != j)) med 0 <= k <= (d-2), og sådan at B(i) = A(j) og B(j) = A(i), så Cross( A(0), A(1), ..., A(d-2) ) = (-1) * Cross( B(0), B(1), ..., B(d-2) ).
Hvis tværs produkt af en bestilt liste af vektorer er ikke-nul, og derefter at krydse produktet har den modsatte tegn på tværs produkt af samme bestilt liste af vektorer med undtagelse af nøjagtigt to vektorer bliver udvekslet (ombyttede).
Da grundlag vektorer e(k), for 0 <= k <= (d-1), (d)-dimensionelle rum, cross produkt af en bestilt liste over (d-1) vektorer kan beregnes ved beregning af determinant af (d*d) matrix nedenfor:
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]
Det afgørende i denne matrix vil ikke være et tal, men vil i stedet være en vektor, hvor der vil være numeriske koefficienter til (d) grundlag vektorer.
For to-dimensionelle rum, og vektoren A = ( ax, ay ):
Cross( A )
= Determinant
(
e(x), e(y),
ax, ay
)
= ay * e(x)
- ax * e(y)
= ( ay, -ax )
Resultatet vektor er vinkelret på det indre input vektor.
Hvis x-y plane opfattes sådan, at x stigninger i en rightward retning og y stigninger i en opadgående retning, så korset produkt Resultatet er en fjerdedel-vending "med uret" i forhold til input vektor.
Eksempler på tværs produkter i to-dimensionelle rum:
Cross( e(x) ) = (-1) * e(y)
Cross( e(y) ) = e(x)
For tre-dimensionelle rum, og vektorer A = ( ax, ay, az ) og 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) )
Eksempler på tværs produkter i tre-dimensionelle rum:
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)
For fire-dimensionelle rum, og vektorer A = ( ax, ay, az, aw ), B = ( bx, by, bz, bw ), og 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) )
Eksempler på tværs produkter i fire-dimensionelle rum:
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)
Nedenstående kode beregner cross-produkt af to eller flere vektorer.
Testen eksempler viser, at ikke-nul resultat vektorerne er vinkelret på alle de input vektorer.
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
"matrix omvendt:" En kvadratisk matrix relateret til en anden square matrix af samme størrelse sådan, at produktet af de to matricer er "identiteten matrix."
Hvis en matrix har en tilknyttet inverse matrix, er der kun én sådan inverse matrix.
Hvis den afgørende faktor for en matrix er lig nul, så den matrix er opkaldt "ental," og matrix ikke har en tilhørende inverse matrix.
Nedenstående kode beregner den inverse af nogen firkantede, ikke-ental matrix.
Værdier for 1*1, 2*2, 3*3, og 4*4 matricer er beregnet ved udtrykkelig formler.
Brug af eksplicitte formler kan være hurtigere end at bruge en almindelig procedure.
Den eksplicitte formler er også pålidelige.
Proceduren for det almindelige tilfælde bruger LU factoring, men der er mange andre metoder til databehandling værdier af matricer.
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 Multiplicative produkt af en matrix og en vektor
De multiplicative produkt af en matrix (M) og en vektor (a) er en vektor (b); (M)*(a) = (b).
Matrix (M) "omdanner" vektor (a) til vektor (b).
Nedenstående kode beregner den multiplicative produkt af en matrix og en vektor.
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 Brug af homogene matricer til at oversætte homogene vektorer
"oversættelse:" En operation T(s), der er angivet med en bestemt vektor (s), som konverterer enhver vektor (r) til en tilsvarende vektor (r + s); T(s) * (r) = (r + s).
T(-s) er omvendt T(s); Inverse(T(s)) = T(-s); T(-s) * (r) = (r - s).
En oversættelse operation kan være repræsenteret af en matrix, sådan at multiplicere oversættelse matrix af en vektor har den virkning at omsætte denne vektor.
Bidragene fra en oversættelse drift til komponenter af resultatet vektor ikke afhænger direkte eller indirekte, om de grundlæggende elementer i den oprindelige vektor.
Således er en oversættelse matrix skal være struktureret på en måde, der gør det muligt at bidrage med en oversættelse til det resultat, uden at blive ramt direkte eller indirekte af de komponenter i den oprindelige vektor.
En "homogen vektoren" er en vektor med en endelig komponent svarende til 1.
En "homogen matrix" er en matrix med en endelig rækken med alle indgange lig med 0 undtagen for en endelig optagelse af 1.
De multiplicative produkt af en homogen matrix og en anden homogen matrix er en homogen matrix.
De multiplicative produkt af en homogen matrix og en homogen vektor er et homogent vektor.
En interessant attribut af dette produkt er, at registreringerne i den sidste kolonne i homogene matrix, med undtagelse af den endelige optagelse, tilføjes direkte til komponenter af resultatet vektor, uden at blive ramt direkte eller indirekte af de komponenter i den oprindelige vektor.
Således er en oversættelse operation kan være repræsenteret af en homogen matrix, og en sådan matrix kan bruges til at oversætte homogene vektorer.
Hvis vi ønsker at repræsentere positioner og retninger i (d)-dimensionelle rum, og vi ønsker at bruge matricer til at omsætte holdninger, så homogen matricer og homogene vektorer kan bruges.
Men den homogene vektorer skal have (d+1) komponenter, og den homogene matricer skal have ((d+1)*(d+1)) poster.
Den sidste del af en homogen vektor og den sidste række i en homogen matrix ikke svarer til nogen af de (d) koordinater, som vi ønsker at repræsentere i (d)-dimensionelle rum, men i stedet kun tjener til at gøre det muligt at få matricer til at gøre nogle operationer uafhængigt af (d) koordinere værdier.
Brug af homogene matricer og homogene vektorer gør det lettere at sammensatte omdannelser og anvende omdannelser med vektorer, men andre operationer bliver mere ubelejligt.
Varetagelse af positioner eller retninger i (d)-dimensionelle rum ved hjælp af homogene vektorer betyder, at vektorer vil have (d+1) komponenter, med den endelige element lig med 1.
Computing punktummet produkt af to homogene vektorer ved hjælp af en regelmæssig dot-produkt funktion ville kombinere implicit punktummet produkt af den første (d) komponenter med en irrelevant værdi af 1 grund til endelige komponenter i 1.
Derfor, for at få dot produktet af de relevante (d) komponenter i homogene vektor (med (d+1) komponenter), dot produkt af homogene vektorer bør ske ved hjælp af en særlig "homogent vektor dot produkt" funktion, der bare hopper den sidste komponent af vektorer.
Alternativt kan den dot produkt kan beregnes med den regelmæssige dot produkt fungere med den forståelse, at et beløb på 1 bør trækkes fra resultatet før fortolke resultatet i relation til (d) relevante koordinater.
Nedenstående kode viser, hvordan en matrix kan dannes, som derefter kan bruges til at oversætte en vektor.
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)-dimensionelle rum: rotationer
"rotation i (d)-dimensionelle rum:" En operation (R) angivet af en to-dimensional rotation plane (S), en rotationsvinkel (t), og en rotation midterste punkt (c), der konverterer ethvert punkt (p) til et tilsvarende punkt (q), sådan at:
(1) Length( q - c ) = Length( p - c ); afstanden mellem det nye punkt og rotation midterste punkt er den samme som afstanden mellem den oprindelige punkt og rotation midterste punkt.
Derfor vil de nye punkt er på overfladen af en særlig (d)-dimensionelle rum;
(2) Lad (P) være en planet, der er parallel med rotation plane (S) og indeholder punkt (p).
Punkt (q) vil blive i planet (P).
Derfor vil de nye punkt er, om en bestemt todimensionalt fly i (d)-dimensionelle rum;
(3) Dot( (q - c), (p - c) ) = Length( q - c ) * Length( p - c ) * Cos( t ); dot produktet af vektoren fra rotation midterste punkt til det oprindelige punkt og vektoren fra rotation midterste punkt til det nye punkt er lig med produktet af længderne af de to vektorer og cosine af rotation vinkel.
Derfor vil de nye punkt er på overfladen af en særlig (d)-dimensional kegle;
Constraint (1) begrænser nyt punkt til overfladen af en særlig (d)-dimensionelle rum.
Betingelse (2) yderligere begrænser nyt punkt til et særligt to-dimensional plane i (d)-dimensionelle rum.
Derfor er det nye punkt skal være på en bestemt todimensionel cirkel i (d)-dimensionelle rum.
Betingelse (3) yderligere begrænser nyt punkt til et bestemt (d)-dimensional kegle.
Derfor er det nye punkt skal være på et bestemt punkt (for Cos( t ) = 1 eller Cos( t ) = (-1)) eller på en af to mulige punkter (for (-1) < Cos( t ) < 1).
Et sidste hindring er forpligtet til at vælge en af de to mulige punkter, som skyldes de begrænsninger, nævnt ovenfor.
Lad T(s) repræsenterer en "oversættelse" operation, hvor (s) er en vilkårlig vektor, således at T(s) * (r) = (r + s), og Inverse(T(s)) * (r) = (r - s).
En rotation R(S,t,c), for plane S, vinkel t, og midterste punkt c, kan udtrykkes i oversættelse og en rotation med et midtpunkt på oprindelse: R(S,t,c) = T(c) * R(S,t,0) * T(-c) = T(c) * R(S,t,0) * Inverse(T(c)).
Lad R(S,t) repræsenterer en rotation med en rotation midterste punkt på oprindelsen af (d)-dimensionelle rum.
Således enhver vilkårlig rotation kan udtrykkes som: R(S,t,c) = T(c) * R(S,t) * T(-c) = T(c) * R(S,t) * Inverse(T(c)).
Lad os overveje to-dimensional rotation fly (S) at omfatte oprindelsen af (d)-dimensionelle rum og er parallel til to koordinere akserne.
Lad (a) og (b) være heltal, der repræsenterer koordinere akserne i (d)-dimensionelle rum, sådan at 0 <= a <= (d-1), og 0 <= b <= (d-1), og (a != b).
Vi definerer R(a,b,t) at være en rotation om oprindelsen af (d)-dimensionelle rum, med en rotation plan parallelt at koordinere akserne angivet ved koordinere indekser (a) og (b), og med en vinkel (t), således at R(a,b,(pi/2)) vil konvertere en positiv koordinere værdi for koordinere akse angivet ved (a) til en positiv koordinere værdi for koordinere akse angives ved (b); eksempel: R(a,b,(pi/2)) * e(a) = e(b), hvor e(k) indikerer en enhed grundlag vektor svarer til at koordinere akse (k).
Denne definition giver den endelige disambiguation nødvendige for at fastlægge de nye punkt fremstillet af en rotation.
En rotation, der flytter punkter langs aksen a mod akse b ifølge en vinkel t svarer til en rotation, der bevæger sig langs aksen b mod akse a ved en vinkel (-t); dvs. R(a,b,t) = R(b,a,(-t)).
Hvis (d >= 2), så (d)-dimensionelle rum har ((d*(d-1))/2) par koordinere akserne.
Således har for (d >= 2), (d)-dimensionelle rum har ((d*(d-1))/2) særskilt akse alliancefrie rotation fly.
(Eksempler: d=2: (x-y plane); d=3: (x-y fly, y-z fly, z-x plane); d=4: (x-y fly, y-z fly, z-x fly, x-w fly, z-w fly, w-y plane).)
Nedenstående kode viser, hvordan en matrix kan bruges til at repræsentere rotation R(a,b,t).
Et eksempel er vist, at bruger en rotation matrix i en homogen måde at rotere en homogen vektor til en ny homogen vektor.
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 Foreslåede definition af "mod uret"
I to-dimensionelle rum, det kvartal-vending rotation, der kan konvertere vektor e(0) til e(1) anses mod uret.
Den rotation er: R(0,1,(pi/2)) * e(0) = e(1).
Derfor er rotation R(0,1,t) er mod uret.
I tre-dimensionelle rum, det kvartal-vending rotation, der kan konvertere vektor e(0) til e(1) anses mod uret.
Kvarteret-vending rotation, der kan konvertere vektor e(1) til e(2) anses også for mod uret.
Kvarteret-vending rotation, der kan konvertere vektor e(2) til e(0) anses også for mod uret.
Derfor er rotationer { R(0,1,t), R(1,2,t), R(2,0,t) } er mod uret.
De Boolesk formlen nedenfor passer definitionen af mod uret for to-dimensional og tre-dimensionelle rum, og ekstrapolere mønstret for den generelle fald i sædskiftet R(a,b,t):
bool counterClockwise =
(
( (a < b) && (((a + b) % 2) == 1) )
|| ( (a > b) && (((a + b) % 2) == 0) )
);
Denne regel vil definere rotationer nedenfor som "mod uret:"
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)-dimensionelle rum: uret rotationer
Brug af den foreslåede definition af "mod uret" i forrige afsnit, rotationer nedenfor er "uret:"
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),
...
For eksempel, i fire-dimensionelle rum, de seks rotation matricer nedenfor repræsenterer uret rotationer:
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)-dimensionelle rum: totalt antal akse alliancefrie retningslinjer
Antallet af særskilt akse alliancefrie retningslinjer i (d)-dimensionelle rum er:
distinctOrientations = Factorial( d ) * IntPow( 2, (d-1) );
Et argument, der understøtter denne formel vises i afsnittene nedenfor.
Eksperimentel beviser kontrollere formlen er også præsenteret nedenfor.
Lad { g(0), g(1), ..., g(d-1) } være et sæt af (d) gensidigt vinkelrette enhed vektorer knyttet til et objekt i (d)-dimensionelle rum.
Disse enhed vektorer angive orientering af objektet i forhold til koordinere akserne i (d)-dimensionelle rum.
Lad "default orientering" af objektet være sådan, at g(i) = e(i) for alle 0 <= i <= (d-1); hver orientering vektor af objektet er forbundet med en særlig enhed grundlag vektor af et koordinatsystem i (d)-dimensionelle rum.
Den første orientering vektor af objektet kan være parallel med nogen af de (d) akser af (d)-dimensionelle rum, og kan være i et af de to retninger i forhold til denne akse.
Således er den første orientering vektor har (2 * d) muligheder.
Den anden orientering vektor af objektet kan være parallel med nogen af de (d-1) resterende akser, med mulighed for to retninger i forhold til denne akse.
Således er den anden orientering vektor har (2 * (d-1)) muligheder.
På samme måde vil det næste orientering vektor har (2 * (d-2)) muligheder.
Den endelige orientering vektor er tvunget til det endelige resterende akse, og er tvunget til en af de to mulige retninger langs denne akse (for at bevare paritet af det sæt orientering vektorer; den afgørende faktor for en bestilt sæt af den orientering vektorer skal være konstant ).
Således er de endelige orientering vektor kun har (1) mulighed.
Det samlede antal kombinationer er givet ved produktet: ((2 * d) * (2 * (d-1)) * ... * (2 * 3) * (2 * 2) * (1)) = (2^(d-1)) * (d*(d-1)*(d-2)*...*2*1) = (2^(d-1)) * Factorial(d).
I afsnittene nedenfor stede beroligende eksperimentelle bevismateriale til støtte for formlen for det samlede antal akse alliancefrie retningslinjer, som et objekt kan være.
Den orientering af et objekt kan være repræsenteret af en matrix (M) hvor vektorerne { g(0), g(1), ..., g(d-1) } er de kolonner af den matrix: