Lineární algebra
Colin Fahey
1. Software
LinearAlgebra.zip
Lineární algebry zdrojový kód (C#)
19910 bytes
MD5: 11d8c8035cac30ba543e5e0b72ee9767
2. Úvod
Tento článek popisuje, vektory a maticemi v (d)-rozměrném prostoru.
3. (d)-rozměrném prostoru: atributy
3.1 Array
"pole:" Soubor proměnných taková, že každá proměnná má jedinečný název, a taková, že jména mohou být přiřazeny objednávky.
Celočíselné hodnoty mohou být použity jako názvy proměnných v poli.
Například v případě, že pole obsahuje (d) proměnných, pak je celé hodnoty { 0, 1, 2, ..., (d-1) } mohou být přiřazeny názvy proměnných v poli.
3.2 (d)-dimenzionální vektorové
"(d)-dimenzionální vektor:" pole na (d) proměnných.
"vektorové složky:" Proměnná ve vektoru.
3.3 (d)-dimenzionální vektorový prostor
"jedno-rozměrném prostoru:" Kompletní soubor hodnot, které mohou být uloženy do jedné proměnné.
"(d)-rozměrném prostoru:" Kompletní sada kombinace hodnot, které mohou být uloženy na několik (d) proměnných.
Formální definice "vektorového prostoru:"
Nechte (T) být základní typy (příklady: reálné číslo, celé číslo, komplexní číslo, racionální číslo, atd.).
Každá proměnná ze základního typu se označuje jako "skalární."
A "(T)-typ (d)-dimenzionální vektorový prostor" je sada (S) na (d)-rozměrné vektory, které mají dvě operace, vektorové navíc (+), a skalární násobení (*), splňující níže uvedené podmínky:
(1) Je-li (v) a (w) se dvěma vektory v (S), pak (v + w) je také vektor v (S);
(2) Je-li (u), (v), a (w) jsou všechny tři vektory v (S), pak (u + v) + w = u + (v + w);
[přídatné látky komutativita]
(3) Je-li (v) a (w) se dvěma vektory v (S), pak (v + w) = (w + v);
[přídatné látky asociativnost]
(4) Je "nula vektoru," (0), v (S), takové, že pro každého vektoru (v) v (S), (v + (0)) = v;
[přídatné látky identitu]
(5) Pokud je (c) jakékoli skalární typu (T), a (v) je jakýkoli vektor v (S), potom produkt (c * v) je vektor v (S);
(6) Pokud (a), (b), a (c) jsou všechny skaláry typu (T), a (v) a (w) jsou všechny vektory v (S), pak (a + b) * v = a*v + b*v, a c*(v + w) = c*v + c*w;
[multiplikativní distributivnost]
(7) Je-li (a) a (b) jsou všechny skaláry typu (T), a (v) je jakýkoli vektor v (S), pak (a*b)*v = a*(b*v);
(8) Je-li "1" je skalární typu (T) taková, že (1*1)=1, a (v) je jakýkoli vektor v (S), pak (1*v) = v;
(9) Pro každý vektor (v) v (S), vektoru (-1)*v = -v splňuje v + (-v) = (0);
[Opačné]
3.4 (d)-rozměrné vektorové kód
Následující kód ukazuje, jak (d)-dimenzionální vektor, s 64-bitovou plovoucí-bod, komponenty, může být provedena.
Pole je dostatečné k tomu, aby představují vektor.
Následující kód má pole obsažené ve třídě, pouze pro pohodlí.
Tento kód není určen jako účinná.
Strukturou (příklad: "struct", jejichž hodnota je typu) zastupující vektory z pevné linky na rozměry (příklady: 3 nebo 4), které by mohly být mnohem účinnější než obecný (d)-dimenzionální vektor třídy naleznete zde.
Ačkoliv níže uvedený kód definuje vektor s plovoucí-bod složek se v tomto dokumentu rovněž využívá vektory se celé složky.
Následující kód může být snadno upraven k provádění vektory se celé složky.
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)-dimenzionálním nulový vektor:" A vektor se všemi (d) komponenty se rovná nule.
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)-rozměrném prostoru: body
"(d)-rozměrném prostoru bod:" V současnosti je přijímána řada (d) proměnné s konkrétní hodnoty "(koordinovat hodnoty)."
"(d)-rozměrném prostoru původu:" pole na (d) všechny hodnoty proměnné se rovnají nule.
3.6 (d)-rozměrné vektory: non-relativní a relativní
A (d)-dimenzionální vektor, který je "ne-relativní" je (d)-dimenzionální vektor, který přímo reprezentuje stav nebo konfigurace.
A (d)-dimenzionální vektor, který je "relativní" je (d)-dimenzionální vektor, které představuje změny v souboru složek.
Relativní vektor může představovat rozdíl mezi dvěma ne-relativní vektory.
Vzhledem k relativní vektor, se stanoví mimo-relativní postavení a konfigurace pomocí tohoto relativního vektoru vyžaduje kombinování, že relativní vektor s non-relativní vektor.
Non-relativní vektory a relativní vektory jsou vektory.
Ať už se jedná o konkrétní vektor, není-relativní nebo relativní musí být upřesněno, kdy vektor je definován.
Pokud se (d)-dimenzionální vektor je interpretována tak, aby ne-relativní, pak je (d)-dimenzionální vektor může představovat místo v (d)-rozměrném prostoru.
3.7 (d)-rozměrné vektorové sčítání, odečítání, a škálování
Vektorové sčítání, odečítání, a škálování:
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)-dimenzionálním základě vektorů
Formální definice "na" bázi vektorového prostoru:
Nechte (T) být základní typy (příklady: reálné číslo, celé číslo, komplexní číslo, racionální číslo, atd.).
Každá proměnná ze základního typu se označuje jako "skalární."
Nechte (V) být "(T)-typ (d)-dimenzionální vektorový prostor."
Pokud non-nulové vektory { u1, u2, ..., ud } v (V) jsou takové, že každý vektor (v) v (V) lze zapsat jako "lineární kombinace" těchto vektorů, v = c1*u1 + c2*u2 + ... + cd*ud, kde { c1, c2, ..., cd } jsou skaláry typu (T), pak (V) je "rozložené" vektory { u1, u2, ..., ud }.
Každá sada non-nulové vektory { u1, u2, ..., ud } "span," že vektorový prostor (V) je pojmenován "na" bázi (V).
Jeden jednoduchý "základ" pro (d)-dimenzionální vektorový prostor je množina (d) odlišné (d)-rozměrné vektory, každý s jedním složky ve výši jedné a všechny ostatní součásti se rovná nule.
Takový základ vektory jsou "ortonormální," což znamená, že jsou-vzájemně kolmé "(ortogonální)" a že každý vektor má délku jednotky.
Každý takový vektor je vektor paralelní jednotky do jedné z (d) koordinuje os.
Vyjádření libovolného vektoru jako "lineární kombinace" těchto vektorů je základ přímo; každou část libovolného vektoru se násobí odpovídající základ vektor, a tyto produkty se sčítají za účelem vytvoření libovolného vektoru.
Následující kód ukazuje, jak se vektor lze vyjádřit jako "lineární kombinaci" na základě vektorů.
Následující kód libovolně definuje ortonormální soubor základě vektory.
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 )
// . . .
}
}
Jakýkoliv vektor v (d)-rozměrném prostoru lze vyjádřit jako součet těchto produktů na základě čísel a vektorů:
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)-rozměrném prostoru: vzdálenost mezi body
Nechte (P) být (d)-dimenzionální vektor, které představuje místo v (d)-rozměrném prostoru.
Nechte (Q) být (d)-dimenzionální vektor, které představuje místo v (d)-rozměrném prostoru.
Nechte (R) být (d)-dimenzionální vektor, který představuje změnu v (d) souřadnice se dostat z bodu (P) bodu (Q); R = (Q - P).
V jednom-rozměrném prostoru, P = ( p0 ), Q = ( q0 ), a R = (Q - P) = ( q0 - p0 ).
Vzdálenost mezi oběma body je: Abs( q0 - p0 ).
Ve dvou-rozměrném prostoru, P = ( p0, p1 ), Q = ( q0, q1 ), a R = (Q - P) = ( q0-p0, q1-p1 ).
Význam dvou kolmých posuvu jako kolmé břehy právo-trojúhelník, vzdálenost mezi body odpovídá délce části přepona tohoto trojúhelníku.
Na Pythagorean vzorec, (a*a) + (b*b) = (c*c), kde (a) a (b) jsou délky kolmicí stranách právo-trojúhelník, a (c) je doba trvání přepona (pořádání straně), může být použita pro určení vzdálenosti mezi dvěma body: Sqrt( Sq(q0-p0) + Sq(q1-p1) ).
Ve tří-rozměrném prostoru, P = ( p0, p1, p2 ), Q = ( q0, q1, q2 ), a R = (Q - P) = ( q0-p0, q1-p1, q2-p2 ).
Význam vysídlení ( q0-p0, q1-p1, 0 ) jako kolmé břehy právo-trojúhelník, a pomocí Pythagorean vzorec, dává vzdálenost mezi bodem (P) a místo ( q0, q1, p2 ): d01 = Sqrt( Sq(q0-p0) + Sq(q1-p1) ).
Posunutí ( q0-p0, q1-p1, 0 ) je kolmo k posuvu ( 0, 0, q2-p2 ), a jiné právo-trojúhelník může být formována, a Pythagorean vzorec lze použít znovu.
Proto je vzdálenost od bodu (P) bodu (Q) je dána: 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) ).
Způsob prodlužuje vzdálenost vzorec ze dvou-rozměrném prostoru až tří-rozměrném prostoru lze aplikovat opakovaně, aby nakonec určit vzdálenost vzorec pro (d)-rozměrném prostoru: Sqrt( (Sq(q0-p0) + Sq(q1-p1)) + Sq(q2-p2) + ... + Sq(qd-pd) ).
Následující kód definuje funkci s názvem "Délka" který počítá délku a (d)-dimenzionální vektor.
Pokud je vektor reprezentuje posunutí mezi dvěma body na (d)-rozměrném prostoru, délky vektoru představuje vzdálenost mezi těmito dvěma body.
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)-rozměrné vektory: skalární součin
"Tečka produktu" převádí dva (d)-rozměrné vektory na řadu.
Následující kód spočítáme tečkou produktu dvou vektorů:
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);
}
// . . .
}
Pro každého vektoru, (A), Length(A) = Sqrt(Dot(A,A)).
3.11 (d)-rozměrné vektory: definice "paralelní"
Vektory (A) a (B) jsou "rovnoběžné," pokud všechny tyto jsou pravdivé:
(1) A.Length() > 0;
(2) B.Length() > 0;
(3) Abs(Dot(A,B)) = A.Length()*B.Length().
Následující kód určuje, je-li pár vektory jsou rovnoběžné (případně anti-sladěna).
Plovoucí-bod čísla mohou hromadit chyby v frakční část vzhledem k omezené přesnosti, a proto, počítačový kód by měl obsahovat non-nulové tolerance při porovnání plovoucí-bod čísel.
Tento kodex obsahuje například hodnoty tolerance, ale například toleranci hodnot nemusí být vhodné pro některé úkoly.
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)-rozměrné vektory: definice "kolmá"
Vektory (A) a (B) jsou "kolmé," jestliže všechny následující jsou pravdivé:
(1) A.Length() > 0;
(2) B.Length() > 0;
(3) Abs(Dot(A,B)) = 0.
Následující kód určuje, je-li pár vektory jsou kolmé. Plovoucí-bod čísla mohou hromadit chyby v frakční část vzhledem k omezené přesnosti, a proto, počítačový kód by měl obsahovat non-nulové tolerance při porovnání plovoucí-bod čísel.
Tento kodex obsahuje například hodnoty tolerance, ale například toleranci hodnot nemusí být vhodné pro některé úkoly.
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 Matice
"matice:" Soubor proměnných taková, že každá proměnná má jedinečnou kombinaci "po" sobě jména, příjmení a názvu "sloupce."
"vstup:" Proměnná v matrice.
Celočíselné hodnoty mohou být použity jako "řádek" jména a názvy "sloupců" pro proměnné v matice.
Například, je-li matice má (totalRows) řádků a sloupců (totalColumns), pak je celé hodnoty { 0, 1, ..., (totalRows-1) } mohou být názvy přiřazeny řádky, a integer hodnoty { 0, 1, ..., (totalColumns-1) } mohou být přiřazeny názvy sloupců.
Což znamená, že proměnné v matici, může být stanoveno, že se stanoví dvojice celá čísla, ( row, column ), což naznačuje, kombinace řádku a sloupci indexy odpovídající konkrétní proměnnou.
Velikost matice je uvedena jako "(totalRows) * (totalColumns)" (nebo "(totalRows) o (totalColumns))."
Toto pořadí rozměry jsou stejné jako pořadí na rozměry použitých upřesnit zápisy v matrice "(( row, column ))."
[Tato úmluva je poněkud nešťastné, protože pro mnohé dva-dimenzionálním použití (příklady: obrázky, grafy, atd.), obvyklé konvence je specifikovat rozměry jako "width * height" a souřadnice jako "( horizontal, vertical )" (nebo "( x, y ))."
To je přesný opak toho řádově rozměry a souřadnice se používá k označení matice a jejich zápisy. ]
Matice s (totalRows) rovná (totalColumns) se jmenuje "náměstí," jinak matrice se jmenuje "obdélníkové."
Matice může být považováno za obsahující soubor "řádek vektory," kde se proměnné v každém řádku jsou interpretovány jako něco, co patří k vektoru.
Matice se lze také obsahuje soubor "sloupce vektory," kde se proměnné v každém sloupci jsou interpretována jako něco, co patří k vektoru.
Matice může reprezentovat celou řadu matematických vztahů.
Význam matice a operace, které by mohly být vhodné pro zpracování záznamů z matice, závisí na kontextu.
Nicméně, existují základní pravidla, na matrice aritmetika, které jsou významné pro mnoho kontextů, a tyto základní pravidla budou stanoveny v následující části.
Pole a (totalColumns) hodnoty jsou dostatečné k tomu, aby představují matice.
Poli může mít (totalRows * totalColumns) proměnných, a vstup na ( row, column ) může odpovídat na pole proměnné v indexu ((totalColumns * row) + column).
Následující kód definuje matici, s 64-bitovou plovoucí-bod zápisy.
Pole a (totalColumns) hodnoty jsou dostatečné k tomu, aby představují matice.
Následující kód je pole a (totalColumns) hodnoty obsažené ve třídě, pouze pro pohodlí.
Následující kód je nelze považovat za efektivní.
Strukturou (příklad: "struct", jejichž hodnota je typu) reprezentující matice na konkrétní pevné rozměry (příklady: 2*2, 3*3, nebo 4*4), které by mohly být mnohem účinnější než obecný (totalRows * totalColumns) třídy naleznete zde.
Ačkoliv níže uvedený kód definuje matice s plovoucí-bod, poznámek, tento článek také využívá matice se celé položky.
Následující kód může být snadno upraven k provádění matice se celé položky.
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
// . . .
}
}
"nulové matice:" matice se všechny záznamy rovnající se nule.
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 navíc, odčítání, násobení a
Matrix sčítání, odečítání a násobení.
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 ]
// . . .
}
}
"identity matrix:" čtvercové matice (celkem řádků se rovná celkové sloupce) se záznamy rovná (1) na diagonální (položky s sebou index se rovná jeho sloupci indexu), a se všemi ostatními zápisy rovná (0).
Pokud čtvercové matice, (M), je vynásobeno "identity matice," (I), na stejný počet řádků (a sloupce), výrobek je rovna (M).
Vynásobením (I) o (M) také vyrábí produkt rovná (M).
Proto "totožnost matice" je analogický k počtu "1" pro násobení čísel (skaláry).
Následující kód vytváří identitu matice s určeným počtem řádků.
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 faktoring, back-substituce
"Matrix LU faktoring:" postup, který převádí čtvercové matice, (M), aby dva nové čtvercové matice, (L) a (U), které mají stejnou velikost jako (M), jako že (L) * (U) = (M), a takové, že matice (U) je "horní trojúhelníková" (všechny položky jsou nižší, než je úhlopříčka nula), a matrice (L) je "nižší trojúhelníkové" (všechny položky nad diagonální jsou nulové).
Matrix LU factoringu je možné použít jako součást větší postup pro řešení systému rovnic, či k nalezení na inverzní matice, nebo k nalezení faktorem matice.
Matrix LU Factoring je alternativa k "Gauss-Jordan odstranění postupem.
Gauss-Jordan odstranění vyžaduje soustava rovnic (A)*(x)=(b), že LU faktoring, pouze vyžaduje, aby matice (A).
Také po stanovení LU factoringu na matice (A), že je velmi snadné zjistit, (x) žádný (b).
Postup při řešení pro vektorové (x) v (A)*(x)=(b), vzhledem k tomu, (b) a LU faktoring (L)*(U)=(A), zahrnuje postup, který bude pojmenován "back-substituční" zde.
Následující kód obsahuje postup pro nalezení LU faktorů na každém náměstí, ne-singulární matice.
Následující kód obsahuje postup pro děláš LU back-substituce.
Upozornění: Počítač kód spočítáme na L a U faktorů na řadě-permutované verzi určité matrice.
Také se L a U matice výsledky jsou kombinovány se do jediného výstupu matice.
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 Matice: determinant
"Určuje se:" A funkci, která převádí žádné náměstí (n * n) matice (A) na řadu, det(A), jako že:
(1) Je-li matice (B) výsledky ve výměně dvou řadách, nebo dva sloupy, na matice (A), pak det(B) = (-(det(A)));
(2) Je-li matice (B) vyplývá z vynásobením každém řádku, sloupci nebo žádné, na matice (A), řadou (c), pak det(B) = c * det(A);
(3) Je-li matice (B) vyplývá z přidáním několika jednoho řádku do druhého řádku, nebo z přidání několika jednoho sloupce do druhého sloupce, pak det(B) = det(A).
(4) Pokud (A) je-horní trojúhelníková matice (A[i,j]=0 pro všechny (i>j)) nebo menší-trojúhelníkové matice (A[i,j]=0 pro všechny (i<j)), potom det(A) = (A[0,0] * A[1,1] * ... * A[n-1,n-1]);
(Například, je určujícím faktorem pro totožnosti matice je jeden, det(I) = 1.)
Pravidla (1), (2), a (3), mohou být používány v procesu nazvaného "Gaussian odstranění, převést jakékoliv čtvercové matice na trojúhelníkové matice.
Článek (4) mohou být použity pro vypočítání faktorem trojúhelníkové matice.
Formální definice "určující:"
Nechte M[n](K) označuje soubor všech matic nad (n * n) oblasti (K).
Na "faktorem" je unikátní funkce (F) taková, že F:M[n](K) --> K se dvěma atributy:
(1) F je střídavý multilineární s ohledem na sloupců (nebo řádků), a
(2) F( I ) = 1.
A "multilineární" funkce (D) o (n * n) matice (A) lze zapsat jako: 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]] ) ), kde součet se bere přes všechny (n^n) kombinace 0 <= k[i] <= (n-1), a kde e[j] představuje řadu (j) o totožnosti matice.
Tyto funkce mohou být charakterizovány hodnotami D( e[k[0]], e[k[1]], ..., e[k[n-1]] ).
"Střídavým multilineární" funkce je multilineární funkce, (D), jako že D( ..., e[i], ..., e[i], ... ) = 0, a D( ..., e[i], e[j], ... ) = (-( D( ..., e[j], e[i], ... ) )).
Proto vyměňování sloupců (nebo řádků) změní znaménko funkce.
Také termín D( e[k[0]], ..., e[k[n-1]] ) je nula v případě, že kombinace k[i] hodnoty nezahrnuje všechny jedinečných hodnot, je-li soubor hodnot není permutací z čísla {0,1,...,(n-1)}.
Pronájem D( I ) = 1 (D( e[0], e[1], ..., e[n-1] ) = 1), v kombinaci se střídavým multilineární atributů popsaných výše, znamená, že všechny hodnoty na střídavý multilineární funkce D( e[k[0]], e[k[1]], ..., e[k[n-1]] ) mohou být stanoveny.
Funkce bude non-nula pouze v případě, že soubor k[i] hodnot je permutace na {0,1,...,(n-1)}, a bude mít velikost jednoho, pokud se nejedná-nula, a toto označení bude počet swapy požadované pro převod permutace na totožnost permutace.
Hodnota určujícím faktorem pro (n * n) matice (A) se může přidáním do produktů formě (sgn(p) * A[0,p[0]] * A[1,p[1]] * ... * A[(n-1),p[n-1]]) pro každého a každý permutace (p) na čísla {0,1,2,...,(n-1)}.
Pojem označuje (sgn(p)) "podpisu" (nebo odkládací count) z permutace (p), kde (sgn(p)) je rovna (+1) li (p) je "ještě permutace," a je rovna (-1) li (p) je "lichá permutace."
Protože tam jsou (n!) permutace čísel {0,1,2,...,(n-1)}, tento vzorec má (n!) summands.
Další postupy k výpočtu faktorem (příklady: těch, které zahrnují Gaussian odstranění, nebo LU faktoring, nebo rozšíření o nezletilé, atd.) implicitně počítají hierarchie množství meziproduktů, které umožňují tyto postupy spočítat determinant v pořadí (n^3) kroky.
Pokud je určujícím faktorem pro matice je nula, matice nemá inverzní.
Pokud je určujícím faktorem pro matice je non-nula, inverzní matice má.
Pokud matice představuje koeficienty lineární systém rovnic, a určujícím, že matice je nulová, pak je soustava rovnic nemá jedinečné řešení.
Pokud matice představuje koeficienty lineární systém rovnic, a určujícím, že matice je non-nula, pak systém rovnic má jedinečné řešení.
Následující kód obsahuje postupy k výpočtu faktorem jakékoliv čtvercové matice.
Faktory pro 1*1, 2*2, 3*3, a 4*4 matice jsou počítány pomocí explicitního vzorce.
Použití explicitní vzorce pro výpočet rozhodujícím pro malé matice je pravděpodobně počítají výsledky rychleji než za použití obecný postup pro výpočet těchto faktorů.
Postup pro obecné použití LU případě faktoringu, ale existuje mnoho dalších způsobů počítání determinantů.
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)-rozměrné vektory: vektorový součin
"vektorový součin:" A funkci spojenou s (d)-rozměrném prostoru, který převádí objednané seznam (d-1) (d)-rozměrné vektory, { A(0), A(1), ..., A(d-2) }, na (d)-dimenzionální vektor, takové, že funkce má atributy níže:
(1) Je-li A(i)=0 (nulový vektor), na každém 0 <= i <= (d-2), pak Cross( A(0), A(1), ..., A(d-2) )=0 (nulový vektor).
Kříž produkt objednané seznam vektorů je nulový, pokud jakýkoli vektor v tomto seznamu je nulová;
(2) Pokud se některá dvojice vektorů, A(i) a A(j), pro jakýkoli (i!=j) s 0 <= i <= (d-2) a 0 <= j <= (d-2), jsou paralelní (Abs(Dot(A(i),A(j))) = Length(A(i)) * Length(A(j))), pak Cross( A(0), A(1), ..., A(d-2) )=0 (nulový vektor).
Kříž produkt objednané seznam vektorů je nulový, jestliže žádné dva vektory v tomto seznamu jsou paralelní;
(3) Je-li každý vektor A(i) je non-nula, a pokud A(i) není souběžně s A(j) pro všechny (i!=j) s 0 <= i <= (d-2) a 0 <= j <= (d-2), pak B = Cross( A(0), A(1), ..., A(d-2) ) je taková, že Dot(A(i),B)=0 pro všechny 0 <= i <= (d-2).
V případě, že výrobek na objednané seznam vektorů je non-nula, pak vektorový součin Výsledkem je kolmá k vektoru v každé objednané seznam vektorů;
(4) Pokud (c) je číselná hodnota, pak 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) ).
Vzhledem k tomu, objednané seznam vektorů, kříž výrobku z objednaných seznam vektorů se některý vektor vynásobený konkrétní číselná hodnota je rovna součinu křížovém z objednaných seznam vektorů by číselná hodnota;
(5) Vzhledem k tomu, dvě čísla (i) a (j), kde (i!=j) a 0 <= i <= (d-2) a 0 <= j <= (d-2), a objednané seznam (d-1) (d)-rozměrné vektory, { A(0), A(1), ..., A(d-2) }, a další objednané seznam (d-1) (d)-rozměrné vektory, { B(0), B(1), ..., B(d-2) }, jako že B(k) = A(k) pro všechny ((k != i) && (k != j)) s 0 <= k <= (d-2), a takové, že B(i) = A(j) a B(j) = A(i), pak Cross( A(0), A(1), ..., A(d-2) ) = (-1) * Cross( B(0), B(1), ..., B(d-2) ).
V případě, že výrobek na objednané seznam vektorů je non-nula, pak to vektorový součin má opačné znaménko na kříži produktu téhož objednat seznamu vektorů s výjimkou přesně dva vektory byly vyměněny (vyměnil).
Vzhledem k tomu, základ vektorů e(k), na 0 <= k <= (d-1), na (d)-rozměrném prostoru, příčný produkt objednané seznam (d-1) vektory mohou být počítány pomocí výpočtu prvkem, který rozhoduje o (d*d) matice níže:
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]
Na faktorem, že matice nebude číslo, ale namísto toho vektoru, kde bude numerických koeficientů na základě (d) vektory.
U dvou-rozměrném prostoru, a vektoru A = ( ax, ay ):
Cross( A )
= Determinant
(
e(x), e(y),
ax, ay
)
= ay * e(x)
- ax * e(y)
= ( ay, -ax )
Výsledkem je vektor kolmý k jednotné vstupní vektor.
Je-li x-y rovině je zobrazena tak, že x zvyšuje ve pravicový směr a y zvyšuje směrem vzhůru, pak vektorový součin výsledek je o čtvrt-zase "po směru hodinových ručiček" vzhledem k vstupní vektor.
Příklady z kříže výrobky ve dvou-rozměrném prostoru:
Cross( e(x) ) = (-1) * e(y)
Cross( e(y) ) = e(x)
Pro tří-rozměrném prostoru, a vektory A = ( ax, ay, az ) a 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) )
Příklady z křížových produktů ve tří-rozměrném prostoru:
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)
U čtyř-dimenzionálním prostoru, a vektory A = ( ax, ay, az, aw ), B = ( bx, by, bz, bw ), a 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) )
Příklady z křížových produktů ve čtyř-dimenzionálním prostoru:
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)
Následující kód spočítáme kříž-výrobků ze dvou nebo více vektory.
Zkušební příklady ukazují, že ne-nulový výsledek vektory jsou kolmé na všechny vstupní vektory.
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: inverzní
"inverzní matice:" čtvercové matice ve spojení s jiným čtvercové matice stejné velikosti takové, že součin dvou matic je "matice identity."
Pokud matice je spojena inverzní matice, existuje jen jedna taková inverzní matice.
Pokud je určujícím faktorem pro matice je nulová, potom se nazývá "singulární" matice, matice a nemá související inverzní matice.
Následující kód spočítáme inverzní na každém náměstí, ne-singulární matice.
Inverses pro 1*1, 2*2, 3*3, a 4*4 matice jsou počítány pomocí explicitního vzorce.
Použití explicitní vzorce by mohly být rychlejší než používání obecné řízení.
Výslovného vzorce jsou také spolehlivé.
Postup pro obecné použití LU případě faktoringu, ale existuje řada jiných metod výpočetní inverses matic.
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 Multiplikativní výrobku na matice a vektorové
Na multiplikativní produkt a matice (M) a vektorové (a) je vektor (b); (M)*(a) = (b).
Matrix (M) "transformuje" vektor (a) k vektoru (b).
Následující kód spočítáme na multiplikativní výrobku na matice a 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 Použití homogenní matice přeložit homogenní vektory
"překlad:" Operace T(s), která určí zejména vektorové (s), že převede jakékoli vektorové (r) na odpovídající vektorové (r + s); T(s) * (r) = (r + s).
T(-s) je inverzí T(s); Inverse(T(s)) = T(-s); T(-s) * (r) = (r - s).
Překlad provoz může být reprezentován jako matice, taková, že vynásobením překlad matice a vektoru má za následek, že vektorový překládáním.
Příspěvky na provoz překlad do složky na výsledek vektor nejsou závislé, přímo nebo nepřímo, v jehož složky původního vektoru.
Proto překlad matice musí být strukturována takovým způsobem, který umožní, aby přispěla překladem do výsledku, aniž by byly ovlivněny, přímo nebo nepřímo, jehož složky původního vektoru.
"Homogenní vektoru" je vektor s konečným částky rovnající se 1.
"Homogenní hmota" je matice s konečným řadě se všechny záznamy rovnající se 0 výjimkou závěrečné vstup 1.
Na multiplikativní produkt a homogenní matice a další homogenní hmota je homogenní matice.
Na multiplikativní produkt a homogenní matrici a homogenní vektoru je homogenní vektor.
Jeden zajímavý atribut tohoto produktu je skutečnost, že zápisy do posledního sloupce na homogenní matice, s výjimkou závěrečné vstupu, se přidávají přímo do složky na výsledek vektor, aniž by byly ovlivněny, přímo nebo nepřímo, jehož složky původní vektor.
Proto překlad provoz může být reprezentován jako homogenní matice, a tato matice je možné použít na překlad homogenní vektory.
Pokud chceme zastupovat míst a směrů v (d)-rozměrném prostoru, a my jsme chtěli využít matice přeložit místa, pak matice homogenní a homogenní vektory mohou být použity.
Nicméně, homogenní vektory musí mít (d+1) komponentů a homogenní matice musí mít ((d+1)*(d+1)) zápisy.
Konečný složkou homogenní vektoru a závěrečná řada homogenní matice, neodpovídají žádné z (d) souřadnice že chceme zastupovat v (d)-rozměrném prostoru, ale místo toho slouží jen proto, aby bylo možné se dostat matice dělat některé operace nezávisle na (d) koordinovat hodnot.
Pomocí matice homogenní a homogenní vektory usnadňuje složené transformace a proměny vztahuje na vektory, ale další operace se staly nevyhovující.
Zastupující míst nebo směrů v (d)-rozměrném prostoru pomocí homogenní vektory znamená, že vektory budou mít (d+1) komponenty, s konečným částky rovnající se 1.
Computing absolutně homogenní produkt dvou vektorů pomocí pravidelných dot-funkce výrobku bude kombinovat implicitně tečkou výrobku z prvních (d) komponenty s irelevantní hodnota 1 konečném důsledku složek 1.
Proto, aby získali dot produktu na příslušné (d) složek z homogenní vektorové (s (d+1) komponenty) se tečka na výrobek homogenní vektory by mělo být provedeno pomocí speciální "homogenní vektorové dot" funkce "výrobku," že prostě vynechává závěrečné části vektory.
Další možností je, že dot produktu by mohl být vypočítán s pravidelnou dot funkce výrobku s tím, že částka na 1 by měla být odečtena od výsledku před interpretaci výsledku jako vztahující se k (d) příslušné souřadnice.
Následující kód ukazuje, jak se matice mohou být tvořeny, který pak může být použita k přeložit 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)-rozměrném prostoru: střídání
"rotace v (d)-rozměrném prostoru:" Operace (R) specifikovány ve dvou-dimenzionálním otočení roviny (S), což je rotace úhel (t) a rotace středu (c), že převádí každém místě (p) k odpovídajícímu bodu (q), jako že:
(1) Length( q - c ) = Length( p - c ); vzdálenost mezi novým bodem a rotace středu je stejná jako vzdálenost mezi původním místě a střídání centrum bod.
Proto byl nový bod, který na povrchu určitého (d)-dimenzionální sféry;
(2) Nechť (P) být rovinou, která je rovnoběžná s rovinou rotace (S) a obsahuje bod (p).
Bod (q) bude v rovině (P).
Proto byl nový bod, který na určitou dva-dimenzionální rovinou v (d)-rozměrném prostoru;
(3) Dot( (q - c), (p - c) ) = Length( q - c ) * Length( p - c ) * Cos( t ); tečka součin vektoru od středu otáčení k výchozímu bodu a vektoru od středu otáčení na nový bod, který se rovná součinu délky těchto dvou vektorů a cosinus rotace úhlu.
Proto byl nový bod, který na povrchu určitého (d)-dimenzionálním kužel;
Constraint (1) omezuje nový bod na povrchu určitého (d)-rozměrné koule.
Podmínka (2) dále omezuje nový bod na konkrétní dva-dimenzionální rovinou v (d)-rozměrném prostoru.
Proto se nový bod, musí být v určitém dva-dimenzionálním kruhu v (d)-rozměrném prostoru.
Podmínka (3) dále omezuje nový bod na konkrétní (d)-dimenzionálním kužel.
Proto se nový bod, musí být na jednom určitém místě (pro Cos( t ) = 1 nebo Cos( t ) = (-1)) nebo na jednom ze dvou možných bodů (pro (-1) < Cos( t ) < 1).
Moje poslední překážkou je povinen vybrat si jednu ze dvou možných míst, které vyplývají z omezení uvedených výše.
Nechte T(s) představuje "překlad" provozu, kde je (s) libovolný vektorový, jako že T(s) * (r) = (r + s), a Inverse(T(s)) * (r) = (r - s).
Rotační R(S,t,c), na rovině S, úhel t a středu c, mohou být vyjádřeny v překladu operací a rotace se středem v místě původu: R(S,t,c) = T(c) * R(S,t,0) * T(-c) = T(c) * R(S,t,0) * Inverse(T(c)).
Nechte R(S,t) představuje střídání provoz s rotací středu na původ (d)-rozměrném prostoru.
Proto libovolného střídání může být vyjádřen jako: R(S,t,c) = T(c) * R(S,t) * T(-c) = T(c) * R(S,t) * Inverse(T(c)).
Uvažujme dva-dimenzionálním střídání rovin (S), které zahrnují původ (d)-rozměrném prostoru a jsou rovnoběžné se dvěma koordinuje os.
Nechte (a) a (b) být celá čísla, které představují koordinuje os v (d)-rozměrném prostoru, jako že 0 <= a <= (d-1), a 0 <= b <= (d-1), a (a != b).
Definujeme R(a,b,t) být rotace kolem původu (d)-rozměrném prostoru, s rotací roviny rovnoběžné, aby koordinovaly koordinovaly osy označeny indexy (a) a (b), a s úhlem (t), jako že by R(a,b,(pi/2)) převést pozitivní koordinovat hodnoty pro koordinaci osy označeny (a) ke kladnému koordinovat hodnoty pro koordinaci osy označeny (b); příklad: R(a,b,(pi/2)) * e(a) = e(b), kde e(k) ukazuje na základě vektorové jednotky odpovídající koordinaci osy (k).
Tato definice se stanoví konečné dvojsmyslů potřebné k určení nový bod, který je produkován buď otočení.
Rotační, že se pohybuje podél osy bodů a k ose b podle úhlu t je rovnocenný s rotací, který se pohybuje podél osy bodů b k ose a o úhlu (-t), tj. R(a,b,t) = R(b,a,(-t)).
Pokud (d >= 2), pak (d)-rozměrném prostoru má ((d*(d-1))/2) párů koordinuje os.
Proto, pro (d >= 2), (d)-rozměrném prostoru má ((d*(d-1))/2) odlišná osa-sladěna střídání rovin.
(Příklady: d=2: (x-y rovinou); d=3: (x-y rovině, y-z rovina, rovina z-x); d=4: (x-y rovině, y-z rovině, z-x rovině, x-w rovině, z-w rovina, rovina w-y).)
Následující kód ukazuje, jak se matice mohou být použity k zastupování střídání R(a,b,t).
Příkladem je prokázáno, že využívá otáčení matice v k homogenním způsobem střídat homogenní vektoru do nového homogenní 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 Navrhovaná definice "proti"
Ve dvou-rozměrném prostoru, čtvrtletí,-zase střídání, které může převést vektorové e(0) k e(1) je považována proti.
Otočení je: R(0,1,(pi/2)) * e(0) = e(1).
Proto, střídání R(0,1,t) je proti.
Ve tří-rozměrném prostoru, čtvrtletí,-zase střídání, které může převést vektorové e(0) k e(1) je považována proti.
Čtvrtletí-zase střídání, které může převést vektorové e(1) k e(2) je také považována proti.
Čtvrtletí-zase střídání, které může převést vektorové e(2) k e(0) je také považována proti.
Proto se rotací { R(0,1,t), R(1,2,t), R(2,0,t) } jsou proti.
Booleova níže uvedeného vzorce zapadá vymezení proti dobu dvou a tří-dimenzionálním-rozměrném prostoru, a extrapolaci toků pro obecný případ rotace R(a,b,t):
bool counterClockwise =
(
( (a < b) && (((a + b) % 2) == 1) )
|| ( (a > b) && (((a + b) % 2) == 0) )
);
Toto pravidlo by se definovat jako "proti" střídání níže:
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)-rozměrném prostoru: rotace ve směru hodinových ručiček
Použití navrhované definice "proti" v předchozím oddílu, střídání níže jsou "ve směru hodinových ručiček:"
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),
...
Například ve čtyř-dimenzionálním prostoru, šest střídání matice níže představují rotací ve směru hodinových ručiček:
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)-rozměrném prostoru: celkový počet osa-sladěna směry
Počet odlišná osa-sladěna směry v (d)-rozměrném prostoru je:
distinctOrientations = Factorial( d ) * IntPow( 2, (d-1) );
Argument na podporu tohoto vzorce se zobrazí v odstavcích níže.
Experimentální ověření důkazů vzorce je rovněž uveden níže.
Nechte { g(0), g(1), ..., g(d-1) } být soubor (d)-vzájemně kolmé jednotkové vektory připojuje se k objektu v (d)-rozměrném prostoru.
Tyto jednotkové vektory ukazují orientaci objektu vzhledem ke koordinaci osy z (d)-rozměrném prostoru.
Nechť "výchozí orientace" objektu je taková, že g(i) = e(i) pro všechny 0 <= i <= (d-1); každý orientaci vektoru z objektů je spojena s výraznou jednotka základě vektoru z souřadnicového systému v (d)-rozměrném prostoru.
První orientace vektoru z objektů může být rovnoběžná s některou z (d) osy z (d)-rozměrném prostoru, a může být v obou dvou směrech vzhledem k ose.
Proto byl první orientaci vektoru má (2 * d) možnosti.
Druhá orientace vektoru z objektů může být rovnoběžná s žádným z (d-1) zbývajících osách, s dvou možných směrech vzhledem k ose.
Proto se druhý orientaci vektoru má (2 * (d-1)) možnosti.
Ve stejné cestě, další orientaci vektoru má (2 * (d-2)) možnosti.
Konečný orientaci vektoru je nucen na konečnou zbývající osy, a je nucen k jednomu ze dvou možných směrech podél této osy (k zachování parity ze souboru orientace vektorů, přičemž určujícím faktorem pro objednané soubor orientaci vektorů musí být konstantní ).
Tudíž konečné orientaci vektoru pouze (1) možnost.
Celkový počet kombinací je dán produkt: ((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).
Odstavců níže přítomné uklidnit a experimentální důkazy o tom, že vzorec pro celkový počet osa-sladěna směry, které v objektu může být.
Orientace na objekt může být reprezentován jako matice (M) kde vektory { g(0), g(1), ..., g(d-1) } jsou sloupce v matici: