Linearna algebra
Colin Fahey
1. Software
LinearAlgebra.zip
Linearna algebra source code (C#)
19910 bytes
MD5: 11d8c8035cac30ba543e5e0b72ee9767
2. Uvod
U ovom se članku opisuje vektore i matrice u (d)-dimenzionalan prostor.
3. (d)-dimenzionalan prostor: atribute
3.1 Array
"polje:" Zbirka varijabli takvih da se svaka varijabla ima jedinstveno ime, i da takva imena mogu biti dodijeljena jedna narudžba.
Integer vrijednosti može se koristiti kao imena za varijable u niz.
Na primjer, ako polje sadrži (d) varijable, onda je cijeli broj vrijednosti { 0, 1, 2, ..., (d-1) } mogu biti nazivi dodjeljuju na varijable u polje.
3.2 (d)-dimenzionalan vektor
"(d)-dimenzionalan vektor:" niz (d) varijabli.
"vektorska komponenta:" varijablu unutar vektor.
3.3 (d)-dimenzionalan vektorski prostor
"jedne dimenzije prostora:" potpuni skup vrijednosti koje mogu biti pohranjene od strane jedne varijable.
"(d)-dimenzionalan prostor:" potpuni skup kombinacije vrijednosti koje mogu biti pohranjene po jedan niz (d) varijabli.
Formalna definicija jedan "vektorski prostor:"
Neka (T) biti osnovna tipa (primjeri: pravi broj, cijeli broj, kompleksni broj, racionalan broj, itd.).
Bilo koje varijable je osnovni tip se zove jedan "skalarni."
A "(T)-type (d)-dimenzionalan vektorski prostor" je skup (S) od (d)-dimenzionalan vektori imaju dvije operacije, vektor Osim (+), i skalar multiplication (*), zadovoljavanju uvjeta u nastavku:
(1) Ako (v) i (w) su bilo koje dvije vektori u (S), zatim (v + w) je također vektor u (S);
(2) Ako (u), (v), i (w) su vektori u bilo koje tri (S), zatim (u + v) + w = u + (v + w);
[additive commutativity]
(3) Ako (v) i (w) su bilo koje dvije vektori u (S), zatim (v + w) = (w + v);
[additive associativity]
(4) Postoji jedna "nula vektor," (0), u (S), tako da za bilo koji vektor (v) u (S), (v + (0)) = v;
[additive identiteta]
(5) Ako je bilo koji skalarni (c) tipa (T), i (v) je bilo koji vektor u (S), zatim proizvod (c * v) je vektor u (S);
(6) Ako (a), (b), i (c) su bilo koje scalars tipa (T), i (v) i (w) su bilo koji vektori u (S), zatim (a + b) * v = a*v + b*v, i c*(v + w) = c*v + c*w;
[multiplikativne distributivity]
(7) Ako (a) i (b) su scalars bilo kojeg tipa (T), i (v) je bilo koji vektor u (S), zatim (a*b)*v = a*(b*v);
(8) Ako "1" je skalar tipa (T) takva da (1*1)=1, i (v) je bilo koji vektor u (S), zatim (1*v) = v;
(9) Za svaki vektor (v) u (S), vektor (-1)*v = -v zadovoljava v + (-v) = (0);
[additive inverse]
3.4 (d)-dimenzionalan vektor kod
Kod ispod pokazuje kako (d)-dimenzionalan vektor, s 64-bitnim komponentama pokretni zarez, može biti implementirana.
Niz je dovoljno da predstavljaju vektor.
Kod ispod ima niz sadržane u razredu, samo za praktičnost.
Kod ne namjerava biti učinkovita.
A structure (primjer: "struct", upišite vrijednost) predstavljaju vektore a fiksni broj dimenzija (primjeri: 3 ili 4) je vjerojatno da će biti puno efikasniji od općeg (d)-dimenzionalan vektor klasa je ovdje prikazano.
Iako je kod ispod definira vektor s pokretni zarez komponente, ovaj dokument također koristi vektore s cijeli broj komponenti.
Kod ispod lako mogu biti izmijenjene kako bi se proveo sa vektori cijeli broj komponenti.
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) nula-dimenzionalan vektor:" vektor sa svim (d) komponente jednaka nuli.
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)-dimenzionalan prostor: boda
"(d)-dimenzionalan prostor point:" niz (d) varijabli s određene vrijednosti "(koordinirati vrijednosti)."
"(d)-dimenzionalan prostor porijekla:" niz (d) varijabli sa svih vrijednosti jednaka nuli.
3.6 (d)-dimenzionalan vektori: non-relativna i relativna
A (d)-dimenzionalan vektor koji je "non-relativna" je (d)-dimenzionalan vektor koji izravno predstavlja status ili oblik.
A (d)-dimenzionalan vektor koji je "relativna" je (d)-dimenzionalan vektor koji predstavlja promjene na set komponenti.
Relativni vektor može predstavljati razliku između dva relativna ne-vektori.
S obzirom na relativnu vektor, koji određuju ne-relativni status ili konfiguracije koje koristite u odnosu vektor zahtijeva da se kombiniranjem relativna vektor s non-relativna vektor.
Non-vektori relativne i relativna vektori su i vektori.
Da li se određeni vektor je non-relativna ili relativna mora biti naveden kada je vektor je definiran.
Ako (d)-dimenzionalan vektor se tumačiti kao da su ne-relativna, onda (d)-dimenzionalan vektor može predstavljati točku u (d)-dimenzionalan prostor.
3.7 (d)-dimenzionalan vektor addition, subtraction, i skaliranje
Vector addition, subtraction, i skaliranje:
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)-dimenzionalan temelju vektori
Formalne definicije na "temelju" jedan vektorski prostor:
Neka (T) biti osnovna tipa (primjeri: pravi broj, cijeli broj, kompleksni broj, racionalan broj, itd.).
Bilo koje varijable je osnovni tip se zove jedan "skalarni."
Neka (V) biti "(T)-type (d)-dimenzionalan vektorski prostor."
Ako je različit od nule vektori { u1, u2, ..., ud } u (V) su takve da je svaki vektor (v) u (V) se može zapisati kao "linearna kombinacija" od one vektore, v = c1*u1 + c2*u2 + ... + cd*ud, gdje { c1, c2, ..., cd } su scalars tipa (T), zatim (V) je "spojeni" by vektori { u1, u2, ..., ud }.
Bilo koji skup ne-nula vektori { u1, u2, ..., ud } da "span" vektorski prostor (V) se zove "temelju" (V).
Jedan jednostavan "osnovu" jednog (d)-dimenzionalan vektorski prostor je skup (d) različita (d)-dimenzionalan vektore, svaki s jedne komponenta jednaka jedan i sve ostale komponente jednaka nuli.
Takva osnova vektori su "orthonormal," što znači da su međusobno-okomica "(ortogonalni)" i da je svaki vektor je jedinica duljine.
Svaki takav vektor je jedinični vektor paralelan jednoj od koordinatnih osi (d).
Izražavanje proizvoljnog vektor kao "linearna kombinacija" tih osnova je vektore direktne; svaku komponentu od proizvoljan vektor je pomnožen s odgovarajućim temelju vektor, i tih proizvoda zbroje da se formira na proizvoljnim vektor.
Kod ispod pokazuje kako se vektor može se izraziti kao "linearna kombinacija" temelju vektori.
Kod ispod proizvoljno definira jedan orthonormal set temelju vektori.
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 )
// . . .
}
}
Bilo koji vektor u (d)-dimenzionalan prostor može se izraziti kao zbroj proizvode brojevi i vektori osnovu:
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)-dimenzionalan prostor: udaljenost između točaka
Neka (P) biti (d)-dimenzionalan vektor koji predstavlja točku u (d)-dimenzionalan prostor.
Neka (Q) biti (d)-dimenzionalan vektor koji predstavlja točku u (d)-dimenzionalan prostor.
Neka (R) biti (d)-dimenzionalan vektor koji predstavlja promjenu u (d) koordinate kako bi dobili (P) iz točke do točke (Q); R = (Q - P).
U jednom-dimenzionalan prostor, P = ( p0 ), Q = ( q0 ), i R = (Q - P) = ( q0 - p0 ).
Udaljenost između dvije točke je: Abs( q0 - p0 ).
U dvodimenzionalni prostor, P = ( p0, p1 ), Q = ( q0, q1 ), i R = (Q - P) = ( q0-p0, q1-p1 ).
Tumačenje dva okomita raseljavanja kao okomita strane pravo-trokut, udaljenost između točaka odgovara na dužini od hypotenuse tog trokuta.
The Pythagorean formula, (a*a) + (b*b) = (c*c), gdje (a) i (b) su duljina od okomica strane pravo-trokut, i (c) je duljina od hypotenuse (izvrtati strani), može se koristiti da biste utvrdili je udaljenost između dvije točke: Sqrt( Sq(q0-p0) + Sq(q1-p1) ).
U trodimenzionalnom prostoru, P = ( p0, p1, p2 ), Q = ( q0, q1, q2 ), i R = (Q - P) = ( q0-p0, q1-p1, q2-p2 ).
Tumačenje displacement ( q0-p0, q1-p1, 0 ) kao okomita strane pravo-trokut, i koristeći Pythagorean formula, daje udaljenost između točke (P) i točka ( q0, q1, p2 ): d01 = Sqrt( Sq(q0-p0) + Sq(q1-p1) ).
The displacement ( q0-p0, q1-p1, 0 ) je okomita na displacement ( 0, 0, q2-p2 ), a drugi s desna trokut može biti formirana, a Pythagorean formula se može koristiti ponovno.
Ovako, udaljenost od točke do točke (P) je (Q) given by: 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) ).
The method of produžuje na udaljenost od formula dvodimenzionalni prostor na trodimenzionalni prostor može se primijeniti više puta kako bi konačno utvrdili udaljenost formula za (d)-dimenzionalan prostor: Sqrt( (Sq(q0-p0) + Sq(q1-p1)) + Sq(q2-p2) + ... + Sq(qd-pd) ).
Kod ispod definira funkcija zove "dužina" koja računa na dužinu od (d)-dimenzionalan vektor.
Kad vektor predstavlja pomak između dvije točke u (d)-dimenzionalan prostor, u dužini od vektor predstavlja udaljenost između ta dva boda.
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)-dimenzionalan vektori: skalarni produkt
Je "skalarni produkt" dva pretvara (d)-dimenzionalan vektori na broj.
Kod ispod jednog računa skalarni produkt dva vektori:
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);
}
// . . .
}
Za bilo koji vektor, (A), Length(A) = Sqrt(Dot(A,A)).
3.11 (d)-dimenzionalan vektori: definicija "parallel"
Vektori (A) i (B) su "paralelno" ako se sve sljedeće su istinite:
(1) A.Length() > 0;
(2) B.Length() > 0;
(3) Abs(Dot(A,B)) = A.Length()*B.Length().
Kod ispod određuje ako par vektori su paralelno (eventualno anti-uskladio).
Floating-point brojevi mogu akumulirati pogreške u djelomični dio s obzirom na ograničene preciznosti, i, dakle, računalo kod bi trebao uključivati ne-nula tolerances Kada uspoređujete s pomičnim-point brojevima.
Kod uključuje primjer tolerancije vrijednosti, ali je primjer tolerancije vrijednosti ne bi mogli biti prikladni za neke zadatke.
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)-dimenzionalan vektori: definicija "okomica"
Vektori (A) i (B) su "okomica" ako sve sljedeće su istinite:
(1) A.Length() > 0;
(2) B.Length() > 0;
(3) Abs(Dot(A,B)) = 0.
Kod ispod određuje ako par su vektori okomit. Floating-point brojevi mogu akumulirati pogreške u djelomični dio s obzirom na ograničene preciznosti, i, dakle, računalo kod bi trebao uključivati ne-nula tolerances Kada uspoređujete s pomičnim-point brojevima.
Kod uključuje primjer tolerancije vrijednosti, ali je primjer tolerancije vrijednosti ne bi mogli biti prikladni za neke zadatke.
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 Matrice
"matrica:" A collection of varijabli takvih da se svaka varijabla ima jedinstvenu kombinaciju za "redom" ime i ime "stupca."
"unos:" varijablu unutar matrica.
Integer vrijednosti može se koristiti kao "red" i imena "stupca" imena za varijable u matrica.
Na primjer, ako je matrica (totalRows) retke i stupce (totalColumns), onda je cijeli broj vrijednosti { 0, 1, ..., (totalRows-1) } mogu biti nazivi dodjeljuju na redaka, i cijeli vrijednosti { 0, 1, ..., (totalColumns-1) } mogu biti dodijeljena imena stupaca.
Dakle, varijable u matrica može biti specificiran po izrijekom par integers, ( row, column ), ukazuju na kombinaciju redak i stupac indeksi odgovaraju specifičnim varijable.
Veličina jedna matrica je naveden kao "(totalRows) * (totalColumns)" (ili "(totalRows) by (totalColumns))."
Ovaj redoslijed je dimenzija ista kao i redoslijed dimenzije koristiti da biste naveli unose u matrica "(( row, column ))."
[Ovaj ugovor je pomalo nesretan, jer je za mnoge dvodimenzionalni koristi (primjeri: slike, grafikoni, etc) uobičajene konvencije je navesti dimenzije kao "width * height" i koordinate kao "( horizontal, vertical )" (ili "( x, y ))."
To je suprotno od reda dimenzija i koordinata koristi za opisivanje matrice i njihova unosa. ]
A matrica s (totalRows) jednaka (totalColumns) se zove "trg," na neki drugi način, matrica se zove "pravokutni."
A matrica može se smatrati sadrži skup "vektori redom," gdje su varijable u svaki red se tumačiti kao da pripadaju jednoj vektor.
A matrica također može smatrati sadrži niz "stupac vektori," gdje su varijable u svakom stupcu se tumačiti kao da pripadaju jednoj vektor.
Matrice može zastupati širok izbor matematičkih odnosa.
Značenje od matrica i operacije koje bi mogle biti prikladne za obradu entries od matrica, ovisi o kontekstu.
Međutim, postoje osnovna pravila matrica aritmetička koje su relevantne za mnogim kontekstima, i tih osnovnih pravila biti će definirani u kasniji dio.
Razvrstati i (totalColumns) vrijednosti su dovoljne da predstavljaju matrica.
U polje može imati (totalRows * totalColumns) varijabli, i ulazak u ( row, column ) možete odgovarati na varijabla u indeks ((totalColumns * row) + column).
Kod ispod definira matrica, s 64-bitnim plutajuće-točka unosa.
Razvrstati i (totalColumns) vrijednosti su dovoljne da predstavljaju matrica.
Kod ispod je razvrstati i (totalColumns) vrijednosti sadržane u razredu, samo za praktičnost.
The code ispod ne namjerava biti učinkovita.
A structure (primjer: "struct", upišite vrijednost) predstavljaju matrice određenih fiksnih dimenzija (primjeri: 2*2, 3*3, ili 4*4) je vjerojatno da će biti puno efikasniji od općeg (totalRows * totalColumns) klasa je ovdje prikazano.
Iako je kod ispod definira matrica s pomičnim-točka unosa, u ovom članku također omogućuje korištenje matrice s cijeli broj unosa.
Kod ispod mogu biti lako modificirana kako bi se proveo cijeli matrice s entries.
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
// . . .
}
}
"nulta matrica:" A matrica sa svim stavkama jednaka nuli.
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 Matrica addition, subtraction, multiplication
Matrica addition, subtraction, multiplication.
public class MatrixF64
{
// . . .
public static MatrixF64 operator +( MatrixF64 a, MatrixF64 b )
{
if ((null == a) || (null == b))
{
return (null);
}
MatrixF64 result = MatrixF64.Zero( a.totalRows, a.totalColumns );
if
(
(null == a.entries)
|| (null == b.entries)
|| (a.totalRows != b.totalRows)
|| (a.totalColumns != b.totalColumns)
)
{
return (result);
}
for (int i = 0; i < a.totalRows; i++)
{
for (int j = 0; j < a.totalColumns; j++)
{
result[i, j] = a[i, j] + b[i, j];
}
}
return (result);
}
public static MatrixF64 operator -( MatrixF64 a, MatrixF64 b )
{
if ((null == a) || (null == b))
{
return (null);
}
MatrixF64 result = MatrixF64.Zero( a.totalRows, a.totalColumns );
if
(
(null == a.entries)
|| (null == b.entries)
|| (a.totalRows != b.totalRows)
|| (a.totalColumns != b.totalColumns)
)
{
return (result);
}
for (int i = 0; i < a.totalRows; i++)
{
for (int j = 0; j < a.totalColumns; j++)
{
result[i, j] = a[i, j] - b[i, j];
}
}
return (result);
}
public static MatrixF64 operator -( MatrixF64 a )
{
if (null == a)
{
return (null);
}
MatrixF64 result = MatrixF64.Zero( a.totalRows, a.totalColumns );
if (null == a.entries)
{
return (result);
}
for (int i = 0; i < a.totalRows; i++)
{
for (int j = 0; j < a.totalColumns; j++)
{
result[i, j] = (-( a[i, j] ));
}
}
return (result);
}
public static MatrixF64 operator *( double scale, MatrixF64 a )
{
if (null == a)
{
return (null);
}
MatrixF64 result = MatrixF64.Zero( a.totalRows, a.totalColumns );
if (null == a.entries)
{
return (result);
}
for (int i = 0; i < a.totalRows; i++)
{
for (int j = 0; j < a.totalColumns; j++)
{
result[i, j] = scale * a[i, j];
}
}
return (result);
}
public static MatrixF64 operator *( MatrixF64 a, double scale )
{
if (null == a)
{
return (null);
}
MatrixF64 result = MatrixF64.Zero( a.totalRows, a.totalColumns );
if (null == a.entries)
{
return (result);
}
for (int i = 0; i < a.totalRows; i++)
{
for (int j = 0; j < a.totalColumns; j++)
{
result[i, j] = scale * a[i, j];
}
}
return (result);
}
public static MatrixF64 operator *( MatrixF64 a, MatrixF64 b )
{
if ((null == a) || (null == b))
{
return (null);
}
MatrixF64 result = MatrixF64.Zero( a.totalRows, b.totalColumns );
// The number of columns of the first matrix must be equal
// to the number of rows of the second matrix.
if
(
(null == a.entries)
|| (null == b.entries)
|| (a.totalColumns != b.totalRows)
)
{
return (result);
}
// Entry (i,j) of the result is equal to the dot product
// of row (i) of (a) and column (j) of (b).
for (int i = 0; i < a.totalRows; i++)
{
for (int j = 0; j < b.totalColumns; j++)
{
double dotProduct = 0.0;
for (int k = 0; k < a.totalColumns; k++)
{
dotProduct += (a[i, k] * b[k, j]);
}
result[i, j] = dotProduct;
}
}
return (result);
}
// . . .
public static void Test()
{
// . . .
// Examples of matrix addition, subtraction, and multiplication:
MatrixF64 a =
new MatrixF64
(
3, 5,
0.0, 1.0, 2.0, 3.0, 4.0, // row 0
5.0, 6.0, 7.0, 8.0, 9.0, // row 1
10.0, 11.0, 12.0, 13.0, 14.0 // row 2
);
MatrixF64 b =
new MatrixF64
(
3, 5,
4.0, -3.0, 2.0, -1.0, 0.0, // row 0
-9.0, 8.0, -7.0, 6.0, -5.0, // row 1
14.0, -13.0, 12.0, -11.0, 10.0 // row 2
);
MatrixF64 c =
new MatrixF64
(
5, 3,
0.0, 1.0, 2.0, // row 0
3.0, 4.0, 5.0, // row 1
6.0, 7.0, 8.0, // row 2
9.0, 10.0, 11.0, // row 3
12.0, 13.0, 14.0 // row 4
);
MatrixF64 result = new MatrixF64();
result = a + b;
result.WriteLine();
// [ 4, -2, 4, 2, 4 ]
// [ -4, 14, 0, 14, 4 ]
// [ 24, -2, 24, 2, 24 ]
result = a - b;
result.WriteLine();
// [ -4, 4, 0, 4, 4 ]
// [ 14, -2, 14, 2, 14 ]
// [ -4, 24, 0, 24, 4 ]
result = -a;
result.WriteLine();
// [ 0, -1, -2, -3, -4 ]
// [ -5, -6, -7, -8, -9 ]
// [ -10, -11, -12, -13, -14 ]
result = 3.0 * a;
result.WriteLine();
// [ 0, 3, 6, 9, 12 ]
// [ 15, 18, 21, 24, 27 ]
// [ 30, 33, 36, 39, 42 ]
result = a * c; // (3x5) * (5x3) = (3x3)
result.WriteLine();
// [ 90, 100, 110 ]
// [ 240, 275, 310 ]
// [ 390, 450, 510 ]
result = c * a; // (5x3) * (3x5) = (5x5)
result.WriteLine();
// [ 25, 28, 31, 34, 37 ]
// [ 70, 82, 94, 106, 118 ]
// [ 115, 136, 157, 178, 199 ]
// [ 160, 190, 220, 250, 280 ]
// [ 205, 244, 283, 322, 361 ]
// . . .
}
}
"identitet matrica:" A kvadratna matrica (ukupni broj redaka iznosi ukupno stupci) s entries jednaka (1) na dijagonale (entries s nizom svojih indeksa jednaka stupca indeksa), i sa svim drugim entries jednaka (0).
Ako kvadratna matrica, (M), pomnožen je jedan "identitet matrica," (I), na isti broj redaka (i stupci), proizvod je jednaka (M).
Množenjem (I) by (M) i proizvodi proizvod jednaka (M).
Ovako, "identitet matrica" je analogan na broj "1" za množenje brojeva (scalars).
Kod ispod stvara identitet matrica s određenog broja redaka.
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-supstitucija
"Matrica LU faktore:" postupak koji pretvara trga matrica, (M), dva nova trga matrices, (L) i (U), istu veličinu kao (M), kao da (L) * (U) = (M), i da takva matrica (U) je "gornji trokutastim" (all entries ispod dijagonale su nula) i matrica (L) je "niži trokutastim" (all entries iznad dijagonale su nula).
Matrica LU faktoring može se koristiti kao dio veće postupak riješiti sustav jednadžbi, ili pronaći inverse od matrica, ili pronaći determinanta jedne matrice.
Matrica LU faktoring je alternativa za uklanjanje "Gauss-Jordan postupak.
Gauss-Jordan uklanjanje zahtijeva sustav jednadžbi (A)*(x)=(b), dok LU faktoriziranje samo zahtijeva matrica (A).
Također, nakon utvrđivanja LU faktoring je matrica (A), vrlo je lako odrediti (x) dati bilo koju (b).
Postupak riješiti za vektor (x) u (A)*(x)=(b), given (b) i LU faktoriziranje (L)*(U)=(A), uključuje postupak koji će se zove "back-supstitucija" ovdje.
Kod ispod uključuje postupak za računanje na LU faktora bilo trgu, ne-singularnih matrica.
Kod ispod uključuje procedure za obavljanje LU back-zamjena.
Oprez: U računalo kod ispod računa, L i U čimbenici za redom-permuted verziju određenog matrica.
Također, L i U matrica Rezultati su kombinirani u jednoj izlazna matrica.
public class MatrixF64
{
// . . .
public static void FindLUFactorsOfARowPermutedVersionOfAnOriginalMatrix
(
MatrixF64 originalMatrix,
ref MatrixF64 LUFactorsMatrix,
ref int[] rowPermutationOfOriginalMatrix,
ref bool rowExchangeParityIsOdd
)
{
// ''LU factoring'' involves factoring a square matrix (M) in to a
// lower-triangular matrix (L), and an upper-triangular matrix (U),
// such that (L)*(U) = (M).
//
// However, the specified matrix (originalMatrix) might not be
// optimal for direct LU factoring, due to the locations of extreme
// values in the matrix. Therefore, this function instead implicitly
// factors a row-permuted version of originalMatrix, where the
// permutation of rows is determined by the values in the specified
// matrix. The specific row permutation selected by this function is
// returned in an array of row indices (rowPermutationOfOriginalMatrix).
// The resulting (L) and (U) factors are merged in to a single
// output matrix (LUFactorsMatrix). Therefore:
//
// (L part of LUFactorsMatrix) * (U part of LUFactorsMatrix)
// = (originalMatrix with rows permuted according to
// rowPermutationOfOriginalMatrix).
//
// Although factoring a row-permuted version of originalMatrix
// makes it more difficult to interpret the results of this function,
// the resulting LU factoring is likely to be more accurate. In a
// sense, this function indicates (via rowPermutationOfOriginalMatrix)
// how to permute the rows of originalMatrix to be able to produce
// the most accurate LU factoring directly, and this function produces
// that factoring (via the LUFactorsMatrix output).
//
// When using the (L) and (U) matrices (merged in to LUFactorsMatrix)
// to solve a system of equations, it is necessary to permute the rows
// of the solution vector of the system of equations according to
// rowPermutation.
//
// The output matrix (LUFactorsMatrix) is a combination of
// the (L) and (U) matrices:
//
// To form the (L) matrix from (LUFactorsMatrix):
// assume (0) for entries above the diagonal,
// assume (1) for entries on the diagonal,
// and use (LUFactorsMatrix) entries below the diagonal.
//
// To form the (U) matrix from (LUFactorsMatrix):
// assume (0) for entries below the diagonal,
// and use (LUFactorsMatrix) entries on and above the diagonal.
//
// The output array (rowPermutationOfOriginalMatrix) indicates how the
// rows of the original matrix have implicitly (not actually) been
// permuted.
//
// The output boolean (rowExchangeParityIsOdd) indicates if the parity
// of the permutation of rows is odd.
//
// This implementation is partly based on pseudo-code appearing in
// ''Introduction to Algorithms'' (Cormen; Lieserson; Rivest;
// 24th printing, 2000; MIT Press; ISBN 0-262-03141-8).
// See the section named ''Overview of LUP decomposition'', in
// chapter 31 (''Matrix Operations''). The implementation here follows
// the ''LUP-Decomposition'' pseudo-code appearing in a section named
// ''Computing an LUP decomposition''.
//
// This implementation also uses ideas found in the implementation of
// ''ludcmp'' in the book ''Numerical recipes in C : The art of
// scientific computing'' (Press; Teukolsky; Vetterling; Flannery;
// second edition; 1996 reprinting; Cambridge University Press).
// The overall loop structure here resembles that of ''ludcmp''.
// The idea of determining and caching row scale factors in advance of
// finding pivots is used here. The method in ''ludcmp'' to avoid zero
// pivot values has been modified here to handle excessively-small
// pivot values, too. HOWEVER, the row permutation indices produced by
// the following implementation is a true permutation of
// {0,1,...,(totalRows-1)}, whereas ''ludcmp'' (and ''lubksb'')
// interpret their ''permutation indices'' as swaps to perform
// (roughly, for each (i), swap b[ i ] and b[ indx[i] ]). So, the
// following implementation of LU-factoring is not directly compatible
// with ''ludcmp'' or ''lubksb''. Also, the following procedure
// does not destroy the original matrix (unlike ''ludcmp'').
double singularMatrixIfMaxRowElementIsLessThanThis = (1.0e-19);
double forcePivotsToHaveAtLeastThisAbsoluteValue = (1.0e-19);
LUFactorsMatrix = null;
rowPermutationOfOriginalMatrix = null;
rowExchangeParityIsOdd = false;
if (null == originalMatrix)
{
return; // No matrix specified.
}
if
(
(originalMatrix.totalRows <= 0)
|| (originalMatrix.totalColumns <= 0)
|| (null == originalMatrix.entries)
)
{
return; // Matrix is empty.
}
if (originalMatrix.totalRows != originalMatrix.totalColumns)
{
return; // Matrix is not square.
}
// Duplicate the original matrix
LUFactorsMatrix = new MatrixF64( originalMatrix );
// (Lines 2-3 of LUP-Decomposition)
// Initialize the row permutation array to the identity
// permutation.
rowPermutationOfOriginalMatrix = new int[LUFactorsMatrix.totalRows];
for (int i = 0; i < LUFactorsMatrix.totalRows; i++)
{
rowPermutationOfOriginalMatrix[i] = i;
}
// For each row, determine the largest absolute value of
// the row elements, and use the reciprocal as the scale for
// that row.
VectorF64 rowScales = VectorF64.Zero( LUFactorsMatrix.totalRows );
for (int i = 0; i < LUFactorsMatrix.totalRows; i++)
{
double largestElementAbsoluteValueInRow = 0.0;
for (int j = 0; j < LUFactorsMatrix.totalColumns; j++)
{
double absoluteElementValue =
Math.Abs( LUFactorsMatrix[i, j] );
if (absoluteElementValue >
largestElementAbsoluteValueInRow)
{
largestElementAbsoluteValueInRow =
absoluteElementValue;
}
}
if (largestElementAbsoluteValueInRow <
singularMatrixIfMaxRowElementIsLessThanThis)
{
return; // Matrix is singular
}
rowScales[i] = (1.0 / largestElementAbsoluteValueInRow);
}
// (Lines 4-18 of LUP-Decomposition)
// Go through the columns of the matrix.
for (int j = 0; j < LUFactorsMatrix.totalColumns; j++)
{
// (Lines 17-18 of LUP-Decomposition)
// (or equation (2.3.12) of NR)
for (int i = 0; i < j; i++)
{
double sum = LUFactorsMatrix[i, j];
for (int k = 0; k < i; k++)
{
sum -= (LUFactorsMatrix[i, k] * LUFactorsMatrix[k, j]);
}
LUFactorsMatrix[i, j] = sum;
}
// Go through all remaining rows and find the best pivot.
double largestScaledSum = 0.0;
int rowIndexOfLargestScaledSum = j;
for (int i = j; i < LUFactorsMatrix.totalRows; i++)
{
// (equation (2.3.13) of NR)
double sum = LUFactorsMatrix[i, j];
for (int k = 0; k < j; k++)
{
sum -= (LUFactorsMatrix[i, k] * LUFactorsMatrix[k, j]);
}
LUFactorsMatrix[i, j] = sum;
double scaledSum = rowScales[i] * Math.Abs( sum );
if (scaledSum >= largestScaledSum)
{
largestScaledSum = scaledSum;
rowIndexOfLargestScaledSum = i;
}
}
// If indeed we found a better pivot, then exchange rows.
if (j != rowIndexOfLargestScaledSum)
{
// (Line 12 of LUP-Decomposition)
// Exchange the row permutation indices
int tempRowIndex = rowPermutationOfOriginalMatrix[j];
rowPermutationOfOriginalMatrix[j] =
rowPermutationOfOriginalMatrix[rowIndexOfLargestScaledSum];
rowPermutationOfOriginalMatrix[rowIndexOfLargestScaledSum] =
tempRowIndex;
// (Lines 13-14 of LUP-Decomposition)
// Exchange the elements of the rows
for (int k = 0; k < LUFactorsMatrix.totalColumns; k++)
{
double temp = LUFactorsMatrix[rowIndexOfLargestScaledSum, k];
LUFactorsMatrix[rowIndexOfLargestScaledSum, k] =
LUFactorsMatrix[j, k];
LUFactorsMatrix[j, k] = temp;
}
// Exchange the row scale factors
double scaleFactor =
rowScales[rowIndexOfLargestScaledSum];
rowScales[rowIndexOfLargestScaledSum] = rowScales[j];
rowScales[j] = scaleFactor;
// Invert the overall row exchange parity
rowExchangeParityIsOdd = (!(rowExchangeParityIsOdd));
}
// Force the pivot element to have at least a certain
// absolute value.
if (Math.Abs( LUFactorsMatrix[j, j] ) <
forcePivotsToHaveAtLeastThisAbsoluteValue)
{
if (LUFactorsMatrix[j, j] < 0.0)
{
LUFactorsMatrix[j, j] =
(-(forcePivotsToHaveAtLeastThisAbsoluteValue));
}
else
{
LUFactorsMatrix[j, j] =
forcePivotsToHaveAtLeastThisAbsoluteValue;
}
}
// If not the final column, then divide all column elements
// below the diagonal by the pivot element, matrixLU[j,j].
// (Lines 15-16 of LUP-Decomposition)
if (j != (LUFactorsMatrix.totalColumns - 1))
{
double reciprocalOfPivot = (1.0 / LUFactorsMatrix[j, j]);
for (int i = (j + 1); i < LUFactorsMatrix.totalRows; i++)
{
LUFactorsMatrix[i, j] *= reciprocalOfPivot;
}
}
}
}
public static void LUBacksubstitution
(
MatrixF64 LUFactorsMatrix,
int[] rowPermutationOfOriginalMatrix,
VectorF64 givenProductVector,
ref VectorF64 solutionVector
)
{
// This implementation is based on pseudo-code appearing in
// ''Introduction to Algorithms'' (Cormen; Lieserson; Rivest;
// 24th printing, 2000; MIT Press; ISBN 0-262-03141-8).
// See the section named ''Overview of LUP decomposition'', in
// chapter 31 (''Matrix Operations''). The implementation here
// follows the ''LUP-Solve'' pseudo-code appearing in a
// section named ''forward and back substitution''.
solutionVector = null;
if (null == LUFactorsMatrix)
{
return; // No matrix specified.
}
if (null == rowPermutationOfOriginalMatrix)
{
return; // No permutation specified.
}
if (null == givenProductVector)
{
return; // No product vector specified.
}
if
(
(null == LUFactorsMatrix.entries)
|| (LUFactorsMatrix.totalRows <= 0)
|| (LUFactorsMatrix.totalColumns <= 0)
)
{
return; // Matrix is empty.
}
if (LUFactorsMatrix.totalRows !=
LUFactorsMatrix.totalColumns)
{
return; // Matrix is not square.
}
if (givenProductVector.Dimensions() !=
LUFactorsMatrix.totalRows)
{
return; // Product vector size not equal to matrix row count.
}
// Copy the product vector in to the result.
solutionVector = new VectorF64( givenProductVector );
// (See LUP-Solve, lines 2-3)
for (int i = 0; i < LUFactorsMatrix.totalRows; i++)
{
double sum =
givenProductVector[rowPermutationOfOriginalMatrix[i]];
for (int j = 0; j < i; j++)
{
sum -= (LUFactorsMatrix[i, j] * solutionVector[j]);
}
solutionVector[i] = sum;
}
// (See LUP-Solve, lines 4-5)
for (int i = (LUFactorsMatrix.totalRows - 1); i >= 0; i--)
{
double sum = solutionVector[i];
for
(
int j = (i + 1);
j < LUFactorsMatrix.totalColumns;
j++
)
{
sum -= (LUFactorsMatrix[i, j] * solutionVector[j]);
}
double diagonalElement = LUFactorsMatrix[i, i];
if (diagonalElement != 0.0)
{
solutionVector[i] = (sum / diagonalElement);
}
}
}
// . . .
public static void Test()
{
// . . .
// Example of LU factoring and back-substitution:
MatrixF64 a3x3 =
new MatrixF64
(
3, 3,
0.0, 1.0, 2.0, // row 0
3.0, 4.0, 5.0, // row 1
6.0, 7.0, 8.0 // row 2
);
int[] rowPermutationOfOriginalMatrix = null;
bool rowInterchangeParityIsOdd = false;
MatrixF64 a3x3LU = null;
MatrixF64.FindLUFactorsOfARowPermutedVersionOfAnOriginalMatrix
(
a3x3,
ref a3x3LU,
ref rowPermutationOfOriginalMatrix,
ref rowInterchangeParityIsOdd
);
a3x3LU.WriteLine();
// [ 6, 7, 8 ]
// [ 0, 1, 2 ]
// [ 0.5, 0.5, 1e-19 ]
// The matrix above is a combination of the L and U
// matrices. The L matrix is 0 above the diagonal,
// 1 on the diagonal, and the shown values below the
// diagonal. The U matrix is 0 below the diagonal,
// and the shown values on and above the diagonal.
// This LU factoring is of a ROW-PERMUTED version
// of the original matrix.
// The following L and U matrices are manually formed
// by inspecting the LU factoring matrix above.
MatrixF64 L3x3 =
new MatrixF64
(
3, 3,
1.0, 0.0, 0.0, // row 0
0.0, 1.0, 0.0, // row 1
0.5, 0.5, 1.0 // row 2
);
MatrixF64 U3x3 =
new MatrixF64
(
3, 3,
6.0, 7.0, 8.0, // row 0
0.0, 1.0, 2.0, // row 1
0.0, 0.0, 0.0 // row 2
);
// The product (L)*(U) produces a ROW-PERMUTED version
// of the original matrix:
result = L3x3 * U3x3;
result.WriteLine( 4 );
// [ 6, 7, 8 ]
// [ 0, 1, 2 ]
// [ 3, 4, 5 ]
VectorF64 b3x1 = new VectorF64( 1.0, 2.0, 3.0 );
VectorF64 x3x1 = null;
MatrixF64.LUBacksubstitution
(
a3x3LU,
rowPermutationOfOriginalMatrix,
b3x1,
ref x3x1
);
// Solution vector
x3x1.WriteLine();
// ( -0.66666667, 1, 0 )
// Both the given b3x1 and this x3x1 are relative to the
// original matrix, a3x3. The a3x3LU for a ROW-PERMUTED
// version of a3x3, but LUBacksubstitution accepts a
// non-permuted product vector (e.g., b3x1) and produces
// a non-permuted solution vector (e.g., x3x1).
// Check solution by multiplying the original matrix a3x3
// by the solution vector x3x1, to get a product vector.
VectorF64 checkb3x1 = a3x3 * x3x1;
checkb3x1.WriteLine();
// ( 1, 2, 3 )
// This matches the b3x1 vector we specified above, so
// the solution is valid.
// . . .
}
}
3.16 Matrix: determinanta
"Određujem:" A funkcija pretvara da je bilo koji trg (n * n) matrica (A) na broj, det(A), tako da:
(1) Ako je matrica (B) rezultate iz razmjenu dva reda, ili dva stupca, a matrica (A), zatim det(B) = (-(det(A)));
(2) Ako je matrica (B) rezultate iz množenjem bilo koji redak, stupac ili bilo koju, a matrica (A), a broj (c), zatim det(B) = c * det(A);
(3) Ako matrica (B) rezultate od dodavanja više od jednog reda u drugi red, ili s dodavanjem više od jednog stupca u drugi stupac, a zatim det(B) = det(A).
(4) Ako (A) je gornji-trokutastim matrica (A[i,j]=0 za sve (i>j)) ili nižim trokutastim-matrica (A[i,j]=0 za sve (i<j)), zatim det(A) = (A[0,0] * A[1,1] * ... * A[n-1,n-1]);
(Na primjer, jedna od odrednica identiteta matrica je jedan; det(I) = 1.)
Pravila (1), (2), i (3), može se koristiti u procesu pod nazivom "Gaussian eliminacija, na bilo koji trg matrica pretvori u trokutastim matrica.
Vladavina (4) se može koristiti kako bi izračunali je determinanta od trokutastim matrica.
Formalna definicija jedna "odrednica:"
Neka M[n](K) označuju skup svih (n * n) matrice nad poljskim (K).
The "determinanta" je jedinstvena funkcija (F) takva da F:M[n](K) --> K sa dva atributa:
(1) F je izmjenična multilinear u odnosu na stupce (ili retke), i,
(2) F( I ) = 1.
A "multilinear" funkcija (D) od (n * n) matrica (A) može biti napisan kao: 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]] ) ), gdje sam je preuzela sve (n^n) kombinacije 0 <= k[i] <= (n-1), i gdje e[j] predstavlja red (j) od matrice identiteta.
Takva funkcija se može karakterizira vrijednosti D( e[k[0]], e[k[1]], ..., e[k[n-1]] ).
An "izmjenične multilinear" funkcija je funkcija multilinear, (D), kao da D( ..., e[i], ..., e[i], ... ) = 0, i D( ..., e[i], e[j], ... ) = (-( D( ..., e[j], e[i], ... ) )).
Ovako, swapanju stupce (ili retke) mijenja predznak funkcije.
Također, pojam D( e[k[0]], ..., e[k[n-1]] ) je nula ako se kombinacija k[i] vrijednosti ne uključuje sve jedinstvene vrijednosti, ako je skup vrijednosti nije permutacijom od brojeva {0,1,...,(n-1)}.
Najam D( I ) = 1 (D( e[0], e[1], ..., e[n-1] ) = 1), u kombinaciji s izmjenične multilinear atributa je gore opisano, znači da su sve vrijednosti od izmjenične multilinear funkcija D( e[k[0]], e[k[1]], ..., e[k[n-1]] ) može biti određena.
U funkciji će biti različit od nule samo ako je skup vrijednosti k[i] je permutacijom od {0,1,...,(n-1)}, i imat će veličine jednog ako je različit od nule, a znak će biti broj swaps dužan da se pretvori u permutacijom na identitet permutacijom.
Vrijednost jedne determinanta od (n * n) matrica (A) može biti kompjutorska dodajući se proizvodi u obliku (sgn(p) * A[0,p[0]] * A[1,p[1]] * ... * A[(n-1),p[n-1]]) za svaki permutacijom (p) na brojeve {0,1,2,...,(n-1)}.
Pojam (sgn(p)) označava "potpis" (ili swap count) s permutacijom (p), gdje (sgn(p)) je jednaka (+1) ako (p) je "čak i permutacijom," te je jednaka (-1) ako (p) je "ak permutacijom."
Budući da postoji (n!) permutations na brojeve {0,1,2,...,(n-1)}, ova formula je (n!) summands.
Ostale procedure za računanje, determinanta (primjeri: one uključuju Gaussian uklanjanje, ili LU faktore, ili za proširenje po maloljetnike, etc) implicitno compute jedan hijerarhiju srednje veličine koje omogućuju one postupke kako bi izračunali su odrednica u poretku (n^3) korake.
Ako je determinanta jedne matrice je nula, matrica nema inverse.
Ako je determinanta jedne matrice je različit od nule, matrica ima inverse.
Ako je matrica predstavlja koeficijenti od linearnih sustava jednadžbi, i da je determinanta matrice nula, tada je sistem jednadžbi nema jedinstveno rješenje.
Ako je matrica predstavlja koeficijenti od linearnih sustava jednadžbi, i da je determinanta matrice se ne-nula, onda sustav jednadžbi ima jedinstveno rješenje.
Kod ispod uključuje postupke za računanje je determinanta bilo kvadratnih matrica.
Odrednice za 1*1, 2*2, 3*3, i 4*4 matrice su kompjutorska by eksplicitno formule.
Korištenje eksplicitno formule izračunati da je za male determinante matrice je vjerojatno da izračunali rezultate brže nego pomoću općeg postupka izračunali da one determinante.
Postupak u slučaju opće koristi LU faktore, ali postoji mnogo drugih metoda za računanje determinante.
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)-dimenzionalan vektori: vektorski produkt
"vektorski produkt:" A funkcije povezane s (d)-dimenzionalan prostor koji pretvara jedan naredio popis (d-1) (d)-dimenzionalan vektori, { A(0), A(1), ..., A(d-2) }, na (d)-dimenzionalan vektor, tako da funkcija ima atribute u nastavku:
(1) Ako A(i)=0 (nula vektor), za bilo koju 0 <= i <= (d-2), zatim Cross( A(0), A(1), ..., A(d-2) )=0 (nula vektor).
Križ proizvod je naredio popis vektori je nula ako se bilo koji vektor u toj listi je nula;
(2) Ako je bilo par vektori, A(i) i A(j), za bilo koju (i!=j) s 0 <= i <= (d-2) i 0 <= j <= (d-2), su paralelno (Abs(Dot(A(i),A(j))) = Length(A(i)) * Length(A(j))), zatim Cross( A(0), A(1), ..., A(d-2) )=0 (nula vektor).
Križ proizvod je naredio popis vektori je nula ako se bilo koje dvije vektori u toj listi su paralelni;
(3) Ako se svaki vektor A(i) je različit od nule, i ako A(i) nije paralelno sa A(j) za sve (i!=j) s 0 <= i <= (d-2) i 0 <= j <= (d-2), zatim B = Cross( A(0), A(1), ..., A(d-2) ) je takva da Dot(A(i),B)=0 za sve 0 <= i <= (d-2).
Ako vektorski produkt je naredio popis vektori je različit od nule, onda je vektorski produkt rezultat je okomit na svaki vektor u naredio popis vektori;
(4) Ako (c) je numerička vrijednost, tada 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) ).
S obzirom na jednom naredio popis vektori, vektorski produkt, naredio popis vektore s bilo kojeg jedan vektor pomnožen s određenim numerička vrijednost je ekvivalentna množenjem vektorski produkt je naredio popis vektori od numerička vrijednost;
(5) S obzirom na dvije brojevi (i) i (j), gdje (i!=j) i 0 <= i <= (d-2) i 0 <= j <= (d-2), i jednom je naredio popis (d-1) (d)-dimenzionalan vektori, { A(0), A(1), ..., A(d-2) }, a drugi je naredio popis (d-1) (d)-dimenzionalan vektori, { B(0), B(1), ..., B(d-2) }, kao da B(k) = A(k) za sve ((k != i) && (k != j)) s 0 <= k <= (d-2), i takva da B(i) = A(j) i B(j) = A(i), zatim Cross( A(0), A(1), ..., A(d-2) ) = (-1) * Cross( B(0), B(1), ..., B(d-2) ).
Ako vektorski produkt je naredio popis vektori je različit od nule, a zatim da vektorski produkt je suprotno znak križa proizvod je naredio da isti popis vektore s izuzetkom točno dva vektori se razmjenjuju (swapati).
S obzirom na temelju vektori e(k), za 0 <= k <= (d-1), (d) dimenzije prostora, vektorski produkt je naredio popis (d-1) vektori mogu biti kompjutorska by computing je determinanta od (d*d) matrica ispod:
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]
Je determinanta matrice od koje neće biti broj, ali će umjesto toga biti vektor, gdje će biti numerički koeficijenti na temelju (d) vektori.
Za dvodimenzionalni prostor, a vektor A = ( ax, ay ):
Cross( A )
= Determinant
(
e(x), e(y),
ax, ay
)
= ay * e(x)
- ax * e(y)
= ( ay, -ax )
Rezultat je vektor okomit na jednu ulazni vektor.
Ako x-y plane gleda kao da x povećava u rightward smjeru i y povećava u jednom smjeru prema gore, zatim križ proizvod rezultat je četvrtina-skrenite "u smjeru kazaljke na satu" u odnosu na ulazni vektor.
Primjeri križ proizvodi u dvodimenzionalni prostor:
Cross( e(x) ) = (-1) * e(y)
Cross( e(y) ) = e(x)
Za trodimenzionalni prostor, a vektori A = ( ax, ay, az ) i B = ( bx, by, bz ):
Cross( A, B )
= Determinant
(
e(x), e(y), e(z),
ax, ay, az,
bx, by, bz
)
= (ay*bz - az*by) * e(x)
+ (az*bx - ax*bz) * e(y)
+ (ax*by - ay*bx) * e(z)
= ( (ay*bz - az*by), (az*bx - ax*bz), (ax*by - ay*bx) )
Primjeri križ proizvodi u trodimenzionalni prostor:
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)
Za četiri-dimenzionalan prostor, a vektori A = ( ax, ay, az, aw ), B = ( bx, by, bz, bw ), i C = ( cx, cy, cz, cw ):
Cross( A, B, C )
= Determinant
(
e(x), e(y), e(z), e(w),
ax, ay, az, aw,
bx, by, bz, bw,
cx, cy, cz, cw
)
= ( ay*bz*cw - ay*bw*cz + az*bw*cy - az*by*cw + aw*by*cz - aw*bz*cy) * e(x)
+ (-ax*bz*cw + ax*bw*cz - az*bw*cx + az*bx*cw - aw*bx*cz + aw*bz*cx) * e(y)
+ ( ax*by*cw - ax*bw*cy + ay*bw*cx - ay*bx*cw + aw*bx*cy - aw*by*cx) * e(z)
+ (-ax*by*cz + ax*bz*cy - ay*bz*cx + ay*bx*cz - az*bx*cy + az*by*cx) * e(w)
= ( ( ay*bz*cw - ay*bw*cz + az*bw*cy - az*by*cw + aw*by*cz - aw*bz*cy),
(-ax*bz*cw + ax*bw*cz - az*bw*cx + az*bx*cw - aw*bx*cz + aw*bz*cx),
( ax*by*cw - ax*bw*cy + ay*bw*cx - ay*bx*cw + aw*bx*cy - aw*by*cx),
(-ax*by*cz + ax*bz*cy - ay*bz*cx + ay*bx*cz - az*bx*cy + az*by*cx) )
Primjeri križ proizvodi u četiri dimenzije prostora:
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)
Kod ispod računa križa-proizvod od dvije ili više vektori.
The test primjeri pokazuju da je rezultat različit od nule su vektori okomit na sve ono što se umeće u vektore.
public class MatrixF64
{
// . . .
public static VectorF64 Cross( params VectorF64[] parameterArray )
{
if (null == parameterArray)
{
return(new VectorF64());
}
int totalVectors = parameterArray.Length;
// The number of dimensions of the space must be one larger
// than the number of specified vectors.
int dimensions = (1 + totalVectors);
if (dimensions < 2)
{
// The cross product does not exist for space with
// fewer than two dimensions.
return (new VectorF64());
}
// All specified vectors must exist and must have a number
// of components equal to the number of dimensions of the space.
for (int k = 0; k < totalVectors; k++)
{
if (null == parameterArray[k])
{
// Vector missing.
return (new VectorF64());
}
else if (null == parameterArray[k].components)
{
// Vector is empty.
return (new VectorF64());
}
else if (parameterArray[k].components.Length != dimensions)
{
// Vector has the wrong number of components.
return (new VectorF64());
}
}
// Form a (d*d) matrix with the (d-1) supplied vectors
// forming the final (d-1) rows of the matrix.
MatrixF64 matrix = MatrixF64.Zero( dimensions, dimensions );
for (int i = 1; i < dimensions; i++)
{
VectorF64 rowVector = parameterArray[ i - 1 ];
for (int j = 0; j < dimensions; j++)
{
matrix[i, j] = rowVector[j];
}
}
VectorF64 result = VectorF64.Zero( dimensions );
// For each component of the result vector, set the first
// row of the matrix to a basis vector, compute the matrix
// determinant, and use the resulting value as the
// component of the result vector. This is inefficient,
// due to the determinant being computed (d) times.
// Clearly it would be better to have explicit formulas
// for 2-dimensional, 3-dimensional, and 4-dimensional
// cross-products.
for (int i = 0; i < dimensions; i++)
{
// Set the first row of the matrix to be the basis
// vector for the coordinate axis with index (i).
for (int j = 0; j < dimensions; j++)
{
matrix[0, j] = 0.0;
}
matrix[0, i] = 1.0;
// Compute the determinant of the modified matrix.
double determinant = matrix.Determinant();
// Set component (i) of the result to the determinant
// that was computed.
result[i] = determinant;
}
return (result);
}
// . . .
public static void Test()
{
// . . .
// Examples of cross products:
// Three-dimensional example:
VectorF64 vcp1x3a = new VectorF64( 1.0, 2.0, 3.0 );
VectorF64 vcp1x3b = new VectorF64( 3.0, 1.0, 2.0 );
VectorF64 vcp1x3 = VectorF64.Cross( vcp1x3a, vcp1x3b );
vcp1x3.WriteLine(); // ( 1, 7, -5 )
// Verify that the result is perpendicular to the original
// vectors:
Log.WriteLine( VectorF64.Perpendicular( vcp1x3, vcp1x3a ) );
Log.WriteLine( VectorF64.Perpendicular( vcp1x3, vcp1x3b ) );
// True, True
// Four-dimensional example:
VectorF64 vcp1x4a = new VectorF64( 1.0, 2.0, 3.0, 4.0 );
VectorF64 vcp1x4b = new VectorF64( 4.0, 1.0, 2.0, 3.0 );
VectorF64 vcp1x4c = new VectorF64( 3.0, 4.0, 1.0, 2.0 );
VectorF64 vcp1x4 = VectorF64.Cross( vcp1x4a, vcp1x4b, vcp1x4c );
vcp1x4.WriteLine(); // ( 4, 4, 44, -36 )
// Verify that the result is perpendicular to the original
// vectors:
Log.WriteLine( VectorF64.Perpendicular( vcp1x4, vcp1x4a ) );
Log.WriteLine( VectorF64.Perpendicular( vcp1x4, vcp1x4b ) );
Log.WriteLine( VectorF64.Perpendicular( vcp1x4, vcp1x4c ) );
// True, True, True
// . . .
}
}
3.18 Matrix: inverse
"matrica inverse:" A kvadratna matrica se odnose na neki drugi trg matrica od iste veličine kao da je proizvod od dvije matrice je "matrica identiteta."
Ako je pridružena matrica matrica inverse, postoji samo jedna takva matrica inverse.
Ako je determinanta jedne matrice nula, tada je matrica se zove "u jednini," i matrica nema pridružena matrica inverse.
Kod ispod računa u bilo koje inverse square, non-singularnih matrica.
Obrnuto za 1*1, 2*2, 3*3, i 4*4 matrice su kompjutorska by eksplicitno formule.
Korištenje eksplicitno formule može biti brže nego koristeći opći postupak.
Na izričiti formule su i pouzdane.
Postupak u slučaju opće koristi LU faktore, ali postoje mnoge druge metode računanja obrnuto od matrice.
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 Varijabli multiplikativne proizvoda od matrica i vektor
The multiplikativne proizvoda od matrica (M) i (a) vektor je vektor (b); (M)*(a) = (b).
Matrica (M) "transformira" vektor (a) na vektor (b).
Kod ispod računa na multiplikativne proizvoda od matrica i 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 Korištenje homogena matrices prevoditi homogena vektori
"prijevod:" an operation T(s), određuje se određeni vektor (s), da se pretvara bilo koji vektor (r) na odgovarajuću vektor (r + s); T(s) * (r) = (r + s).
T(-s) je inverse of T(s); Inverse(T(s)) = T(-s); T(-s) * (r) = (r - s).
A prijevod operacija može biti zastupljena je matrica, kao da se množenjem matrica od prevođenje vektor ima učinak prevodio da vektor.
The contributions of prijevod rad na komponente rezultat vektor ne ovise, izravno ili neizravno, na komponente izvorne vektor.
Dakle, prijevod matrica mora biti strukturiran na način koji omogućuje da se priložiti prijevod na rezultat bez da su pogođeni, izravno ili neizravno, od strane komponente originalni vektor.
A "homogena vektor" je vektor sa završnom komponenta jednaka 1.
A "homogena matrica" je matrica s konačni red s all entries jednaka 0 osim za konačni ulazak 1.
The multiplikativne proizvoda od homogena matrica i drugi homogena matrica je matrica homogena.
The multiplikativne proizvoda od homogena matrica i homogena vektor je vektor homogena.
Jedan interesantan atribut ovog proizvoda je da se unose u posljednjem stupcu na homogena matrica, s izuzetkom konačne unos, dodao je izravno na komponente rezultat vektor, bez da su pogođeni, izravno ili neizravno, od strane komponente originalni vektor.
Dakle, prijevod operacija može biti zastupljena je homogena matrica, i takva matrica može se koristiti za prevođenje homogena vektori.
Ako želimo predstavljaju pozicije i smjera u (d)-dimenzionalan prostor, a mi želimo koristiti matrice prevoditi pozicije, a zatim istovrstan homogeneous vektori matrice i može se koristiti.
Međutim, homogena vektori moraju imati (d+1) komponente, i homogena matrice moraju imati ((d+1)*(d+1)) entries.
Krajnji komponenta je homogena vektor i konačni red je homogena matrica ne odgovaraju bilo kojem (d) koordinate koje želimo predstaviti u (d)-dimenzionalan prostor, ali umjesto služe samo kako bi ga moguće dobiti matrices učiniti neke operacije neovisno od (d) koordinirati vrijednosti.
Korištenje homogena matrice i vektore homogena olakšava compound transformacije i transformacije se primjenjuju na vektore, ali druge radnje postati više inconvenient.
Predstavljanje pozicije ili smjera u (d)-dimenzionalan prostor koristeći homogena vektore znači da su vektori će imati (d+1) komponentama, sa završnom komponenta jednaka 1.
Računala i skalarni produkt dva homogena vektori koristeći obične dot-product funkcija bi kombinirati implicitno je skalarni produkt od prvih (d) komponente s jedne nevažne vrijednost 1 zbog konačne komponente 1.
Stoga, da bi se skalarni produkt od relevantnih (d) komponente od homogena vektor (s (d+1) komponente), dot product of homogeneous vektori bi trebao biti ispunjavanja koristeći poseban "homogena vektor skalarni produkt" funkcija koje jednostavno preskace konačne komponenta u vektore.
Alternativno, skalarni produkt mogao biti kompjutorska s redovnim skalarni produkt funkcija s razumijevanjem da je iznos od 1 bi trebalo oduzeti od rezultat prije prevođenja kao rezultat koji se odnose na relevantne (d) koordinate.
Kod ispod pokazuje kako se matrica može biti formirana da mogu onda se prevesti 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)-dimenzionalan prostor: ophodnje
"rotacija u (d)-dimenzionalan prostor:" an operation (R) naveden je dvodimenzionalni rotacija ravnine (S), rotation angle (t), i rotacija centra mjesta (c), koji pretvara bilo koju točku (p) na odgovarajuće mjesto (q), tako da:
(1) Length( q - c ) = Length( p - c ), udaljenost između nove točke i rotacija centra mjesta je ista kao udaljenost između izvorne točke i rotacija centra mjesta.
Ovako, novi je točka na površini određenoj (d)-dimenzionalan sferi;
(2) Neka (P) biti ravnina koja je paralelno s rotacija ravnine (S) i sadrži točku (p).
Point (q) će biti u ravnini (P).
Ovako, novi točka je na određeni dvodimenzionalni plane u (d)-dimenzionalan prostor;
(3) Dot( (q - c), (p - c) ) = Length( q - c ) * Length( p - c ) * Cos( t ), skalarni produkt, vektorski od centra rotacije točke do točke i originalni vektor od centra rotacije točke na novu stvar je jednaka proizvodu, dužina od ta dva vektore i kosinusna od rotacije kut.
Ovako, novi je točka na površini određenoj (d) dimenzije membrana;
Ograničenje (1) ograničava novi pokažite na površini određenoj (d)-dimenzionalan sferi.
Stanje (2) dalje constrains novu točku na određenu dvodimenzionalni plane u (d)-dimenzionalan prostor.
Stoga, nova točka mora biti na određenoj dvodimenzionalni krug u (d)-dimenzionalan prostor.
Stanje (3) dalje constrains novu točku na određenu (d) dimenzije membrana.
Stoga, nova točka mora biti na jednoj određenoj točki (za Cos( t ) = 1 ili Cos( t ) = (-1)) ili na jedan od dva moguća bodova (za (-1) < Cos( t ) < 1).
Jedan final ograničenje je dužan da se odabere jedan od dva moguća boda, da rezultat od navedenog ograničenja.
Neka T(s) predstavlja "prijevod" operacija, gdje (s) je proizvoljan vektor, tako da T(s) * (r) = (r + s), i Inverse(T(s)) * (r) = (r - s).
A rotation R(S,t,c), za S ravnine, kut t, i centra mjesta c, može se izražava u smislu prevođenja operacija i rotacija s centra točka na podrijetla: R(S,t,c) = T(c) * R(S,t,0) * T(-c) = T(c) * R(S,t,0) * Inverse(T(c)).
Neka rotacija R(S,t) predstavlja rad s centra rotacije točke na podrijetlo (d)-dimenzionalan prostor.
Ovako, bilo koje proizvoljne rotacije može se izraziti kao: R(S,t,c) = T(c) * R(S,t) * T(-c) = T(c) * R(S,t) * Inverse(T(c)).
Let us razmislite dvodimenzionalni rotacija ravnine (S) da uključite podrijetlo (d)-dimenzionalan prostor te se paralelno na dvije coordinate axes.
Neka (a) i (b) biti integers koje predstavljaju coordinate axes u (d)-dimenzionalan prostor, kao da 0 <= a <= (d-1), i 0 <= b <= (d-1), i (a != b).
We define R(a,b,t) biti rotacija o podrijetlu (d)-dimenzionalan prostor, s rotacijom plane parallel to coordinate axes naznačeni su koordinirati indeksi (a) i (b), i sa kuta (t), kao da bi se pretvoriti R(a,b,(pi/2)) pozitivno koordinirati vrijednost za koordinirati osi naznačeni su (a) na pozitivne vrijednosti za koordinatno coordinate axis naznačeni su (b); primjer: R(a,b,(pi/2)) * e(a) = e(b), gdje e(k) označava jedinicu na temelju odgovarajuće koordinatni vektor osi (k).
Ova definicija daje završnu Wikicitata koje je potrebno kako bi u potpunosti utvrditi nove točke koje proizvodi jedan vrtnja.
A što se miče rotacija točke uzduž osi a prema osi b prema kuta t je ekvivalent za smjenjivanje što se miče točke uzduž osi b prema osi a po jedan kut (-t); tj. R(a,b,t) = R(b,a,(-t)).
Ako (d >= 2), zatim (d)-dimenzionalan prostor je ((d*(d-1))/2) parova od koordinatnih osi.
Ovako, za (d >= 2), (d)-dimenzionalan prostor je ((d*(d-1))/2) razlikuje os-uskladio rotacija ravnine.
(Primjeri: d=2: x-y plane (); d=3: (x-y plane, y-z plane, z-x plane); d=4: (x-y plane, y-z plane, z-x plane, x-w plane, z-w plane, w-y plane).)
Kod ispod pokazuje kako se matrica može se koristiti da zastupaju rotation R(a,b,t).
Jedan primjer je pokazala da koristi matrica rotacije u homogena način to okretati se homogena vektor u novi vektor homogena.
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 Predložena definicija "suprotno kazaljci sata"
U dvodimenzionalni prostor, četvrtina-turn rotation da može pretvoriti u vektor e(0) da e(1) smatra suprotno kazaljci sata.
Vrtnja je: R(0,1,(pi/2)) * e(0) = e(1).
Stoga, vrtnja R(0,1,t) je suprotno kazaljci sata.
U trodimenzionalnom prostoru, četvrtina-turn rotation da može pretvoriti u vektor e(0) da e(1) smatra suprotno kazaljci sata.
The četvrtina-turn rotation da može pretvoriti u vektor e(1) da e(2) je također smatraju suprotno kazaljci sata.
The četvrtina-turn rotation da može pretvoriti u vektor e(2) da e(0) je također smatraju suprotno kazaljci sata.
Stoga, rotacije { R(0,1,t), R(1,2,t), R(2,0,t) } se suprotno kazaljci sata.
The Booleova formula ispod odgovara definiciji suprotno kazaljci sata za dvodimenzionalni i trodimenzionalni prostor, i extrapolates uzorak za opće slučaju vrtnja R(a,b,t):
bool counterClockwise =
(
( (a < b) && (((a + b) % 2) == 1) )
|| ( (a > b) && (((a + b) % 2) == 0) )
);
Ovo pravilo bi se definirati kao "suprotno kazaljci sata" ophodnje u nastavku:
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)-dimenzionalan prostor: rotacija u smjeru kazaljke na satu
Korištenje predloženih definicija "suprotno kazaljci sata" u prethodnom odjeljku, rotacije su "u smjeru kazaljke na satu" u nastavku:
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),
...
Na primjer, u četiri dimenzije prostora, šest rotacija matrice ispod predstavljaju rotacije u smjeru kazaljke na satu:
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)-dimenzionalan prostor: ukupni broj osi-uskladio ORIJENTACIJE
Broj razlikuje os-uskladio orijentacije u (d)-dimenzionalan prostor je:
distinctOrientations = Factorial( d ) * IntPow( 2, (d-1) );
An argument podupiranju ovog formula se pojavljuje u paragrafima koji slijede.
Eksperimentalne provjere dokaza formula je također prezentiran u nastavku.
Neka { g(0), g(1), ..., g(d-1) } biti skup međusobno (d)-vektori okomit jedinica priključena na jedan objekt u (d)-dimenzionalan prostor.
Ove jedinice vektori ukazuju na orijentaciju objekta u odnosu na koordinatni osovine od (d)-dimenzionalan prostor.
Neka je "zadana orientation" of objekt biti takva da g(i) = e(i), za sve 0 <= i <= (d-1); orientation vektor za svaki objekt je povezana s jedinicom temelju vektor razlikuje se od koordinatni sustav u (d)-dimenzionalan prostor.
Prvi vektor orientation of objekt može biti paralelno s bilo kojim od (d) osovine od (d)-dimenzionalan prostor, a može biti u jednom od dva smjera u odnosu na koje osi.
Ovako, prvi orientation vektor ima (2 * d) opcije.
Drugi vektor orientation of objekt može biti paralelno na bilo koji od preostalih (d-1) osovine, sa dva moguća smjera u odnosu na koje osi.
Dakle, drugi vektor orijentiran je (2 * (d-1)) opcije.
Na isti način, pored orijentacije vektor ima (2 * (d-2)) opcije.
Krajnji orientation vektor je prisilio na konačno preostale osi, te je prisilio na jedan od dva moguća smjera koji uz os (čuvati u paritet, skup orientation vektori, determinanta je naredio set od orijentacije vektori mora biti konstanta ).
Dakle, konačna opredjeljenja vektor samo je (1) opciju.
Ukupni broj kombinacija je given by proizvoda: ((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).
U paragrafima koji slijede predstaviti ohrabrujuća eksperimentalni dokazi u prilog formula za ukupni broj osi-uskladio orijentacije u kojem jedan objekt može biti.
Orijentacija jedan objekt može biti zastupljena je matrica (M) gdje vektori { g(0), g(1), ..., g(d-1) } su stupci od matrica:
M = g(0)[ 0 ], g(1)[ 0 ], ..., g(d-2)[ 0 ], g(d-1)[ 0 ],
g(0)[ 1 ], g(1)[ 1 ], ..., g(d-2)[ 1 ], g(d-1)[ 1 ],
g(0)[ 2 ], &nb