Линейна Алгебра
Colin Fahey
1. Софтуер
LinearAlgebra.zip
Линейна алгебра изходния код (C#)
19910 bytes
MD5: 11d8c8035cac30ba543e5e0b72ee9767
2. Въведение
В тази статия се описва вектори и матрици (d) в тримерно пространство.
3. (d) тримерно пространство: атрибути
3.1 Масив
"ред:" А събирането на променливите, така че всяка променлива е уникално име, както и такива, че имената могат да бъдат възложени поръчка.
Цяло число стойности могат да се използват като имена на променливите в един масив.
Така например, ако един масив съдържа (d) променливи, тогава цялото { 0, 1, 2, ..., (d-1) } стойности могат да бъдат възложени на имената на променливите в масив.
3.2 (d) тримерно векторно
"(d) тримерно векторно:" пълен спектър от (d) променливи.
"вектор компонента:" променлив вектор в едно.
3.3 (d) тримерно векторно пространство
"едно тримерно пространство:" Пълният набор от ценности, които могат да бъдат съхранявани от една променлива.
"(d) тримерно пространство:" Пълният набор от комбинации от стойности, които могат да бъдат съхранени от един ред на (d) променливи.
Формалната дефиниция на "векторно пространство:"
Нека (T) бъде основен вид (примери: реален брой, число, комплексно число, рационално число, и т.н.).
Всяка променлива на основен вид се нарича един "scalar."
А "(T) тип (d) тримерно векторно пространство" е набор (S) на (d) тримерно вектори като двете операции, вектор допълнение (+) и scalar размножаване (*), отговарящи на условията по-долу:
(1) Ако (v) и (w) са носители на две по-(S), а после (v + w) също е вектор в (S);
(2) Ако (u), (v) и (w) са всякакви носители в три (S), а после (u + v) + w = u + (v + w);
[добавка Комутативност]
(3) Ако (v) и (w) са носители на две по-(S), а после (v + w) = (w + v);
[добавка асоциативност]
(4) Налице е "нулев вектор," (0), в (S), така че за всеки вектор (v) в (S), (v + (0)) = v;
[добавка за самоличност]
(5) Ако (c) е всяко scalar от типа (T) и (v) е всеки вектор в (S), тогава продуктът (c * v) е вектор в (S);
(6) Ако (a), (b) и (c) са всякакъв вид scalars на (T) и (v) и (w) са всички носители на (S), а после (a + b) * v = a*v + b*v и c*(v + w) = c*v + c*w;
[Мултипликативна distributivity]
(7) Ако (a) и (b) са всякакъв вид scalars на (T) и (v) е всеки вектор в (S), а после (a*b)*v = a*(b*v);
(8) Ако "1" е scalar от типа (T) такава, че (1*1)=1 и (v) е всеки вектор в (S), а после (1*v) = v;
(9) За всеки вектор (v) в (S), вектора (-1)*v = -v отговаря v + (-v) = (0);
[добавка обратния]
3.4 (d) тримерно векторно код
Кода по-долу показва как един (d) тримерно векторно, с 64-битов плаваща точка компоненти, могат да бъдат приложени.
Една матрица е достатъчно, за да представлява вектор.
Кода по-долу е един масив, съдържащи се в един клас, само за удобство.
Кодът не е предназначено да бъде ефективна.
Структура (например: "struct", чиято стойност е тип), представляващи носители на определен брой размери (примери: 3 или 4) е вероятно да бъде много по-ефективно, отколкото общите (d) тримерно векторно клас е показано тук.
Макар и кода по-долу дефинира вектор с плаваща точка компоненти, този документ да се възползва от вектори с цялото компоненти.
Кода по-долу, могат лесно да бъдат модифицирани за изпълнение на вектори с цялото компоненти.
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) тримерно нулев вектор:" A вектор с всички компоненти (d) равна на нула.
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) тримерно пространство: точки
"(d) тримерно пространство точка:" пълен спектър от (d) променливи с конкретни стойности "(координират стойности)."
"(d) тримерно пространство произход:" пълен спектър от (d) променливи с всички стойности, равни на нула.
3.6 (d) тримерно вектори: не са относителни и в относително
А (d) тримерно векторно, че "не" е "роднина" е (d) тримерно векторно, които пряко да представлява състояние или конфигурация.
А (d) тримерно векторно, че е "роднина" е (d) тримерно векторно, който представлява промени набор от компоненти.
Относително вектор може да представлява разликата между два не-роднина носители.
Като се има предвид относително вектор, определяне не-роднина статус или конфигурация, че използването на относителния вектор изисква комбиниране на този относителен вектор с не-роднина вектор.
Не-относителна вектори и относителната вектори са носители.
Дали дадена вектор не е роднина или роднина трябва да бъде уточнен, когато се дефинира вектор.
Ако (d) тримерно векторно се тълкува като не са относителни, тогава (d) тримерно векторно може да представлява една точка в (d) тримерно пространство.
3.7 (d) тримерно векторно Освен това, логаритмична, и мащабиране
Вектор Освен това, логаритмична, и мащабиране:
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) тримерно база вектори
Формалната дефиниция на "базата" на векторно пространство:
Нека (T) бъде основен вид (примери: реален брой, число, комплексно число, рационално число, и т.н.).
Всяка променлива на основен вид се нарича един "scalar."
Нека (V) бъде "(T) тип (d) тримерно векторно пространство."
Ако не-нула вектори { u1, u2, ..., ud } в (V) са такива, че всеки вектор (v) в (V) могат да бъдат написани като "линейна комбинация" от тези вектори, v = c1*u1 + c2*u2 + ... + cd*ud, където { c1, c2, ..., cd } са scalars от типа (T), а после (V) е "обхватът" с вектори { u1, u2, ..., ud }.
Всеки набор от не-нула вектори { u1, u2, ..., ud }, които "покриват" векторно пространство (V) се нарича "основа" на (V).
Един прост "основа" на (d) тримерно векторно пространство е съвкупност от (d) различим (d) тримерно вектори, всяка от които с един компонент, равен на една и всички останали компоненти, равна на нула.
Тези вектори са "ортонормирана" база, което означава, че те са взаимно-перпендикулярни "(ортогонална)" и че всеки вектор е единица дължина.
Всеки такъв вектор е единица вектор паралел с една от (d) координира оси.
Изразяване на произволен вектор като "линейна комбинация" на базата на тези вектори е пряка; всеки компонент на произволен вектор е умножен по съответната база вектор, и тези продукти се добавят заедно да формират произволен вектор.
Кода по-долу показва как един вектор може да се изрази като "линейна комбинация" на базата носители.
Кода по-долу произволно определя един набор от ортонормирана база вектори.
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 )
// . . .
}
}
Всеки вектор в тримерно пространство (d) може да се изрази като сума на продуктите на номера, както и база вектори:
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) тримерно пространство: разстояние между точките
Нека (P) бъде (d) тримерно векторно, който представлява една точка в (d) тримерно пространство.
Нека (Q) бъде (d) тримерно векторно, който представлява една точка в (d) тримерно пространство.
Нека (R) бъде (d) тримерно векторно, която представлява промяната в (d) координати, за да получите от (P) точка към точка (Q); R = (Q - P).
В един тримерно пространство, P = ( p0 ), Q = ( q0 ) и R = (Q - P) = ( q0 - p0 ).
Разстоянието между две точки е: Abs( q0 - p0 ).
В двумерните пространство, P = ( p0, p1 ), Q = ( q0, q1 ) и R = (Q - P) = ( q0-p0, q1-p1 ).
Тълкуване на две перпендикулярни премествания на перпендикулярни страни на десния триъгълник, разстоянието между точките съответства на дължината на hypotenuse на този триъгълник.
В Pythagorean формула, (a*a) + (b*b) = (c*c), където (a) и (b) са дължини на страни перпендикулярни на десния триъгълник, както и (c) е дължината на hypotenuse (skew странични), могат да бъдат използвани за определяне на разстоянието между две точки: Sqrt( Sq(q0-p0) + Sq(q1-p1) ).
В триизмерното пространство, P = ( p0, p1, p2 ), Q = ( q0, q1, q2 ) и R = (Q - P) = ( q0-p0, q1-p1, q2-p2 ).
Тълкуване на разселването ( q0-p0, q1-p1, 0 ) като перпендикулярни страни на десния триъгълник, както и използването на Pythagorean формула, дава разстоянието между точка (P) и точката ( q0, q1, p2 ): d01 = Sqrt( Sq(q0-p0) + Sq(q1-p1) ).
В разселване ( q0-p0, q1-p1, 0 ) е перпендикулярна на разселването ( 0, 0, q2-p2 ), и друг десен триъгълник могат да се образуват и Pythagorean формула може да се използва отново.
По този начин разстоянието от точка (P) точка (Q) се дава с: 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) ).
Методът за удължаване на разстоянието от формулата двумерен пространство за тримерно пространство могат да бъдат прилагани многократно в крайна сметка да се определи разстоянието формула за (d) тримерно пространство: Sqrt( (Sq(q0-p0) + Sq(q1-p1)) + Sq(q2-p2) + ... + Sq(qd-pd) ).
Кода по-долу дефинира функция, наречена "дължина," която computes дължината на (d) тримерно векторно.
Когато един вектор представлява преместване между две точки в (d) тримерно пространство, дължината на вектор представлява разстоянието между тези две точки.
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) тримерно вектори: точкова продукт
В "точка продукт" конвертира две (d) тримерно носители на определен брой.
Кода по-долу computes една точка от два продукта вектори:
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);
}
// . . .
}
За всеки вектор, (A), Length(A) = Sqrt(Dot(A,A)).
3.11 (d) тримерно вектори: определение на "паралелни"
Вектори (A) и (B) са "паралелни," ако всеки един от следните са верни:
(1) A.Length() > 0;
(2) B.Length() > 0;
(3) Abs(Dot(A,B)) = A.Length()*B.Length().
Кода по-долу, ако прецени, двойка вектори са паралелни (евентуално анти-изравни).
Плаващ точка номера могат да се натрупват грешки в дробна част се дължи на ограничена точност, и, следователно, компютърен код трябва да включва не-нула отклонения при сравняването на плаваща точка числа.
Кодът включва например толерантност стойности, но например толерантност стойности може да не е подходящ за някои задачи.
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) тримерно вектори: определение на "перпендикулярна"
Вектори (A) и (B) са "перпендикулярни," ако всички изброени по-долу са верни:
(1) A.Length() > 0;
(2) B.Length() > 0;
(3) Abs(Dot(A,B)) = 0.
Кода по-долу, ако прецени, двойка вектори са перпендикулярни. Плаващ точка номера могат да се натрупват грешки в дробна част се дължи на ограничена точност, и, следователно, компютърен код трябва да включва не-нула отклонения при сравняването на плаваща точка числа.
Кодът включва например толерантност стойности, но например толерантност стойности може да не е подходящ за някои задачи.
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 Вероятностите
"Матрица:" A събиране на променливите, така че всяка променлива е уникално съчетание на "реда" и име на "колона" име.
"влизане:" променлива в матрицата.
Цяло число стойности могат да бъдат използвани като "ред" имената и "колоната" за имената на променливите в матрицата.
Например, ако една матрица е (totalRows) редове и колони (totalColumns), тогава цялото { 0, 1, ..., (totalRows-1) } стойности могат да бъдат имената, възложени на редове, както и цялото { 0, 1, ..., (totalColumns-1) } стойности могат да бъдат имената, възложени на колони.
Така, един променливи в една матрица може да бъде определен като посочите една двойка числа, ( row, column ), посочвайки комбинация от ред и колона показатели, отговарящи на специфичните променлива.
Размерът на матрицата е посочена като "(totalRows) * (totalColumns)" (или "(totalRows) от (totalColumns))."
Това порядъка на размерите е същият като размери от порядъка на Използва се за задаване на вписвания в матрица "(( row, column ))."
[Тази конвенция е малко неблагоприятни, тъй като за много двумерен употреба (примери: изображения, графики и т.н.) обичайните конвенцията е да се определят размерите "width * height," както и координати, както "( horizontal, vertical )" (или "( x, y ))."
Това е обратното на реда, с размери и координати, използвана за описание на вероятностите и техните вписванията. ]
Матрица с (totalRows) равен (totalColumns) се нарича "квадратна," в противен случай, матрицата се нарича "правоъгълна."
Матрица може да се разглежда като съдържа набор от "ред носители," когато променливи във всеки ред се тълкува като принадлежащи към един вектор.
Матрица може да се счита, съдържащ набор от "колона носители," когато променливи във всяка колона се тълкува като принадлежащи към един вектор.
Вероятностите може да представлява голямо разнообразие от математически отношения.
Същността на една матрица, както и дейностите, които могат да бъдат подходящи за обработка на записите на матрица, зависи от контекста.
Все пак има основни правила на матрицата аритметиката, които са подходящи за много контексти, както и на тези основни правила, ще бъдат определени в един следващ раздел.
Множество и (totalColumns) стойност са достатъчни, за да представляват матрица.
В масив може да има (totalRows * totalColumns) променливи, както и влизането в ( row, column ) да отговаря на множество променливи в индекса ((totalColumns * row) + column).
Кода по-долу дефинира матрица, с 64-битов плаваща точка вписванията.
Множество и (totalColumns) стойност са достатъчни, за да представляват матрица.
Кода по-долу е един масив, както и (totalColumns) стойност, съдържаща се в един клас, само за удобство.
Кода по-долу не е предназначено да бъде ефективна.
Структура (например: "struct", чиято стойност е тип), представляващи матрици на специално фиксирани размери (примери: 2*2, 3*3, или 4*4) е вероятно да бъде много по-ефективно, отколкото общите (totalRows * totalColumns) клас е показано тук.
Макар че определя кода по-долу една матрица с плаваща точка вписвания, в тази статия също се възползва от матрици с цялото вписванията.
Кода по-долу могат лесно да бъдат модифицирани за изпълнение на матрици с цялото вписванията.
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
// . . .
}
}
"нулева матрица:" матрица с всички вписвания, равна на нула.
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 Матрицата Освен това, логаритмична, и размножаване
Матрицата Освен това, логаритмична, и размножаване.
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 ]
// . . .
}
}
"самоличността матрица:" квадратна матрица (общо редове се равнява на общия колони) по вписванията, равна на (1) по диагонал (записи с ред индекс, равен на своята колона на индекса), както и с всички други вписвания, равна на (0).
Ако един квадратен матрица, (M), се умножава с "лична матрица," (I), от същия брой на редовете (и колони), продуктът е равна на (M).
Умножаване (I) от (M) също така да произвежда продукт, равен на (M).
По този начин, "идентичността матрица" е аналогична на броя "1" за размножаването на номера (scalars).
Кода по-долу създава идентичност матрица с определен брой редове.
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 Матрицата: LU факторинг; назад-заместваща
"Матрицата LU факторинг:" процедура, която конвертира квадратна матрица, (M), до два нови квадратни матрици, (L) и (U), като със същия размер като (M), така че (L) * (U) = (M), както и такива, че матрицата е (U) "горната триъгълна" (всички вписвания под диагонал са нула), както и матрица (L) е "долна триъгълна" (всички вписвания над диагонала са нула).
Матрицата LU факторинг може да се използва като част от по-голяма процедура за решаване на система от уравнения, или да намерят обратния на матрица, или да намери определящ фактор на матрица.
Матрицата LU факторинг е алтернатива на "Gauss-Jordan елиминиране на процедурата.
Gauss-Jordan елиминиране изисква система от уравнения (A)*(x)=(b), че LU факторинг само изисква матрица (A).
Освен това, след определяне на LU факторинг на матрица (A), че е много лесно да се определи (x) даде (b).
Условията и редът за разрешаване за векторно (x) в (A)*(x)=(b), като се има предвид (b) и LU факторинг (L)*(U)=(A), свързано с процедурата, която ще се нарича "задна заместване" тук.
Кода по-долу включва процедура за изчисляване на LU фактори за всяко поле, които не са в единствено матрица.
Кода по-долу включва процедура за извършване на LU задна заместване.
Внимание: Компютърът кода по-долу computes на L и U фактори на реда-permuted версия на определена матрица.
Също така, L и U матрица резултати са комбинирани в един единствен изход матрица.
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 Матрицата: определящият
"определящ фактор:" A функция, която конвертира всякакви квадратни (n * n) матрица (A) до номер, det(A), например, че:
(1) Ако матрицата (B) резултатите от обмена на два реда, или от две колони, на матрица (A), а после det(B) = (-(det(A)));
(2) Ако матрица (B) резултати от умножаване на всеки ред или колона, на матрица (A), от редица (c), а после det(B) = c * det(A);
(3) Ако матрицата (B) резултати от добавяне на повече от един ред към друг ред, или от добавяне на повече от една колона на друга колона, а след това det(B) = det(A).
(4) Ако (A) е-горна триъгълна матрица (A[i,j]=0 за всички (i>j)) или по-ниска-триъгълна матрица (A[i,j]=0 за всички (i<j)), а после det(A) = (A[0,0] * A[1,1] * ... * A[n-1,n-1]);
(Например, определящ фактор за идентичност е една матрица; det(I) = 1.)
Правила (1), (2) и (3), могат да бъдат използвани, в процес на елиминиране на име "Gaussian, да конвертирате всеки квадратен матрица с триъгълна матрица.
Правило (4) може да се използва за изчисляване на определящ фактор на триъгълна матрица.
Формалната дефиниция на "определящ фактор:"
Нека M[n](K) означаваме множеството на всички (n * n) матрици над полето (K).
В "определящ фактор" е уникална функция (F) такава, че F:M[n](K) --> K с две качества:
(1) F е променлив multilinear във връзка с колони (или редове); и,
(2) F( I ) = 1.
А "multilinear" функция (D) на (n * n) матрица (A) може да бъде написано като: 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]] ) ), когато сумата е поела всички (n^n) комбинации от 0 <= k[i] <= (n-1), и когато e[j] представлява ред (j) на самоличността матрица.
Тази функция може да се характеризира със стойностите на D( e[k[0]], e[k[1]], ..., e[k[n-1]] ).
Една "променлив multilinear" функция е multilinear функция, (D), така че D( ..., e[i], ..., e[i], ... ) = 0 и D( ..., e[i], e[j], ... ) = (-( D( ..., e[j], e[i], ... ) )).
Така, гащеризон вместо колони (или редовете) промени в знак на функция.
Освен това, понятието D( e[k[0]], ..., e[k[n-1]] ) е нула, ако комбинацията от k[i] ценностите не включва всички уникални стойности; ако набор от ценности, не е permutation на номерата {0,1,...,(n-1)}.
Отдаване под наем D( I ) = 1 (D( e[0], e[1], ..., e[n-1] ) = 1), в комбинация с променлив multilinear атрибути, описани по-горе, означава, че всички стойности на променливите multilinear функция D( e[k[0]], e[k[1]], ..., e[k[n-1]] ) може да се определи.
Функцията ще бъде не-нула, само ако набор от ценности k[i] е permutation на {0,1,...,(n-1)}, и ще има размера на една, ако тя не е нула, и ще се подпише на броя на суапове, необходими за превръщане на permutation на permutation идентичност.
Стойността на един определящ фактор за един (n * n) матрица (A) могат да бъдат изчислени чрез добавяне на продукти от формата (sgn(p) * A[0,p[0]] * A[1,p[1]] * ... * A[(n-1),p[n-1]]) за всеки и всеки permutation (p) на номерата {0,1,2,...,(n-1)}.
Терминът означава (sgn(p)) "подписването" (или суап разчита) на permutation (p), където (sgn(p)) е равна на (+1) ако (p) е "дори permutation," и е равна на (-1) ако (p) е "нечетен permutation."
Тъй като съществуват (n!) permutations на номерата {0,1,2,...,(n-1)}, тази формула е (n!) summands.
Други процедури за изчисляване на определящ фактор (примери: тези, които включват Gaussian елиминиране или LU факторинг, или разширяване с малолетни и т.н.) подразбиращо се изчисли йерархия от междинни количества, които да позволяват на тези процедури, за да изчислявате определящ фактор в един ред на (n^3) стъпки.
Ако определящ фактор на матрица е нула, матрицата не са обратни.
Ако определящ фактор на матрицата не е нула, матрицата е една инверсия.
Ако представлява матрицата на коефициентите на линейно уравнение на системата, както и определящ фактор на тази матрица е нула, а след това системата от уравнения не разполага с уникално решение.
Ако представлява матрицата на коефициентите на линейно уравнение на системата, както и определящ фактор на матрицата, че не е нула, а след това системата от уравнения е уникално решение.
Кода по-долу включва процедури за изчисляване на определящ фактор за всеки квадратен матрица.
Детерминанти за 1*1, 2*2, 3*3 и 4*4 вероятностите са изчислени с изричното формули.
Използване на ясни формули за изчисляване на детерминанти за малки матрици е вероятно да изчислявате резултатите бързо, отколкото използването на обща процедура за изчисляване на тези фактори.
Процедурата за общия случай LU ползва факторинг, но има много различни методи за изчисляване детерминанти.
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) тримерно вектори: по продукт
"кръстосано продукт:" функция, свързани с (d) тримерно пространство, което конвертира един поръчан списък на (d-1) (d) тримерно носители, { A(0), A(1), ..., A(d-2) } до (d) тримерно векторно, такава, че е функция на атрибути по-долу:
(1) Ако A(i)=0 (нулев вектор), за всеки 0 <= i <= (d-2), а после Cross( A(0), A(1), ..., A(d-2) )=0 (нулев вектор).
Общите продукт на една поръчани списък на векторите е нулев, ако всеки вектор в този списък е нула;
(2) Ако всяка двойка вектори, A(i) и A(j), за всеки (i!=j) с 0 <= i <= (d-2) и 0 <= j <= (d-2), са паралелни (Abs(Dot(A(i),A(j))) = Length(A(i)) * Length(A(j))), а после Cross( A(0), A(1), ..., A(d-2) )=0 (нулев вектор).
Общите продукт на една поръчани списък на векторите е нулев, ако всеки две вектори в този списък са паралелни;
(3) Ако всеки вектор A(i) не е нула, и ако A(i) не е паралелно с A(j) за всички (i!=j) с 0 <= i <= (d-2) и 0 <= j <= (d-2), а после B = Cross( A(0), A(1), ..., A(d-2) ) е такава, че Dot(A(i),B)=0 за всички 0 <= i <= (d-2).
Ако кръста продукт на една поръчани списък на векторите не е нула, а след това кръста продукт, резултат е перпендикулярна на всеки вектор подредени в списък от вектори;
(4) Ако (c) е цифрова стойност, а след това 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) ).
Като се има предвид, наредени един списък на носители, кръста продукт на поръчаната списък с всички носители на един вектор, умножен по специална цифрова стойност е еквивалентна на умножи по продукт, поръчан от списъка на носители на цифрова стойност;
(5) Като се има предвид, две числа (i) и (j), където (i!=j) и 0 <= i <= (d-2) и 0 <= j <= (d-2), както и списък на подредени (d-1) (d) тримерно носители, { A(0), A(1), ..., A(d-2) }, както и друга подредени списък на (d-1) (d) тримерно носители, { B(0), B(1), ..., B(d-2) }, така че B(k) = A(k) за всички ((k != i) && (k != j)) с 0 <= k <= (d-2), както и такива, че B(i) = A(j) и B(j) = A(i), а после Cross( A(0), A(1), ..., A(d-2) ) = (-1) * Cross( B(0), B(1), ..., B(d-2) ).
Ако кръста продукт на една поръчани списък на векторите не е нула, а след това, че продуктът е по обратния знак на кръста продукт на същата тази поръчка списък на вектори, с изключение на точно две бацилоносители се обменя (замяна).
Като се има предвид, база вектори e(k), за 0 <= k <= (d-1), на (d) тримерно пространство, кръста продукт на една поръчани списък на (d-1) вектори могат да бъдат изчислени чрез изчисляване на определящ фактор на (d*d) матрицата по-долу:
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]
В определящ фактор за тази матрица няма да бъде много, но вместо това ще бъде един вектор, където ще има числен коефициенти за (d) база вектори.
За двумерен пространство и векторно A = ( ax, ay ):
Cross( A )
= Determinant
(
e(x), e(y),
ax, ay
)
= ay * e(x)
- ax * e(y)
= ( ay, -ax )
В резултат на вектор е перпендикулярна на единен вход вектор.
Ако x-y равнина се възприема такъв, че x се увеличава в посока rightward и y увеличения в една възходяща посока, след това кръста продукт Резултатът е една четвърт завой "по часовниковата стрелка," отнасящо се до входа вектор.
Примери за кръстосани продукти в двумерен пространство:
Cross( e(x) ) = (-1) * e(y)
Cross( e(y) ) = e(x)
За тримерно пространство, както и носители A = ( ax, ay, az ) и 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) )
Примери за кръстосани продукти в тримерно пространство:
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)
За срок от четири тримерно пространство, както и носители A = ( ax, ay, az, aw ), B = ( bx, by, bz, bw ) и 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) )
Примери за кръстосани продукти в четири тримерно пространство:
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)
Кода по-долу computes моментното продукт на две или повече вектори.
Изпитът примери показват, че не нулев резултат вектори са перпендикулярни на всички носители на входа.
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 Матрицата: обратния
"обратната матрица:" квадратна матрица, свързани с друга квадратни матрици с еднакъв размер, така че продуктът на двете матрици е "самоличността матрица."
Ако матрицата е асоцииран обратната матрица, има само един такъв обратната матрица.
Ако определящ фактор на матрица е нула това означава, матрицата се нарича "единствено" и матрицата не е свързано обратната матрица.
Кода по-долу computes обърнат на площада, не-единствено матрица.
Обратно за 1*1, 2*2, 3*3 и 4*4 вероятностите са изчислени с изричното формули.
Използване на ясни формули може да бъде бързо, отколкото чрез общия ред.
Изричното формули също са надеждни.
Процедурата за общия случай LU ползва факторинг, но има много различни методи за изчисляване обратно пропорционални на вероятностите.
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 Мултипликативна продукт на матрица и вектор
В Мултипликативна продукт на матрица (M) и (a) вектор е вектор (b); (M)*(a) = (b).
Матрицата (M) "трансформира" вектор (a) за векторно (b).
Кода по-долу computes на Мултипликативна продукт на матрица и вектор.
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 Използването на хомогенна вероятностите да преведат хомогенна вектори
"превод:" операция T(s), определен от специално вектор (s), че всяко конвертира векторен (r) на съответния вектор (r + s); T(s) * (r) = (r + s).
T(-s) е обърнат на T(s); Inverse(T(s)) = T(-s); T(-s) * (r) = (r - s).
Превод операция може да бъде представлявана от една матрица, като се умножи, че преводът матрица с вектор е ефектът на превръщането на този вектор.
Приносът на операция за превод на компонентите на вектора резултат не зависи пряко или косвено, върху компонентите на оригинала вектор.
Така, превод матрицата, трябва да бъдат структурирани по начин, който позволява тя да допринесе превод отношение на резултата, без да са засегнати пряко или косвено, от компонентите на оригинала вектор.
"Хомогенна вектор" е вектор с крайно компонент, равен на 1.
"Хомогенна матрица" е матрица с окончателния ред с всички вписвания, равна на 0 с изключение на окончателен проект на 1.
В Мултипликативна продукт на хомогенна матрица и друга хомогенна матрица е хомогенна матрица.
В Мултипликативна продукт на хомогенна матрица и хомогенна вектор е хомогенна вектор.
Един интересен атрибут на този продукт е, че текстовете в окончателния колона на хомогенна матрица, с изключение на окончателния проект, се добавя директно към компонентите на вектора резултат, без да са засегнати пряко или косвено, от компонентите на първоначалното вектор.
Така, превод операция може да бъде представляван от хомогенна матрица, както и такава матрица може да се използва за да превеждате хомогенна носители.
Ако ние искаме да представляват позиции и направления в (d) тримерно пространство, и ние искаме да използват матрици да преведат позиции, а след това хомогенна матрици и хомогенна вектори могат да бъдат използвани.
Въпреки това, трябва да са хомогенни вектори (d+1) компоненти, както и хомогенна вероятностите трябва да са ((d+1)*(d+1)) вписванията.
Окончателният компонент на хомогенна вектор и окончателния ред на хомогенна матрицата не отговарят на някое от (d) координати, които желаят да представят в (d) тримерно пространство, но вместо да служат само за да направят възможно да получите вероятностите да направите някои операции независимо от (d) координира ценности.
Използването на хомогенна хомогенна вектори и матрици улеснява сложни преобразувания и се прилага трансформации с бацилоносители, но други операции стават все по-Неудобната.
Представляващи позиции или направления в (d) тримерно пространство използване на хомогенни вектори означава, че ще има носители (d+1) компоненти, с крайната компонент, равен на 1.
Компютърни на точкова продукт от две хомогенни носители използват редовно точкова продукт функция би съчетал имплицитно на точкова продукт на първата (d) компоненти с неподходящи стойност на 1 дължи на компонентите на крайното 1.
Следователно, за да получите точкова продукт на съответната (d) компоненти на хомогенна вектор (с (d+1) компоненти), на точкова продукт на хомогенна бацилоносители трябва да бъде направено чрез специален "хомогенна точкова продукт вектор" функция, която просто Пропуснати крайния компонент на вектори.
Освен това, в точка продуктът може да бъде изчислена с редовни точкова продукт функция с разбирането, че размерът на 1 трябва да бъдат извадени от резултата, преди да тълкуване на резултата, както свързани с (d) съответните координати.
Кода по-долу показва как една матрица може да се формира, които могат да бъдат използвани след това да преведат вектор.
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) тримерно пространство: ротации
"ротация в (d) тримерно пространство:" операция (R) определен от двуизмерен ротация равнина (S), един ъгъл на завъртане (t), както и ротация на централната точка (c), че преобразува всяка точка (p) до съответния пункт (q), например, че:
(1) Length( q - c ) = Length( p - c ); разстоянието между нова точка и ротация на централната точка е същата като разстоянието между първоначалното точка и ротация на централната точка.
Така, новата точка е на повърхността на дадено (d) триизмерна сфера;
(2) Да (P) е равнина, която е успоредна на равнината (S) ротация и съдържа точка (p).
Точка (q) ще бъдат в равнината (P).
Така, новата точка е по-специално двумерен равнина, в (d) тримерно пространство;
(3) Dot( (q - c), (p - c) ) = Length( q - c ) * Length( p - c ) * Cos( t ); на точкова продукт на вектора на въртене централната точка на първоначалната точка и вектор от ротация център точка за новия пункт е равен на произведението от тези две дължини на вектори и степенуване на въртене ъгъл.
Така, новата точка е на повърхността на дадено (d) тримерно конус;
Принуда (1) границите на новата точка на повърхността на конкретен (d) триизмерна сфера.
Състояние (2) допълнително constrains новата точка към даден двумерен равнина (d) в тримерно пространство.
Ето защо, новият пункт трябва да бъде по-специално двумерен кръг (d) в тримерно пространство.
Състояние (3) допълнително constrains новата точка за конкретен (d) тримерно конус.
Ето защо, новият пункт трябва да бъде на едно специално точка (за Cos( t ) = 1 или Cos( t ) = (-1)) или на една от двете възможни точки (за (-1) < Cos( t ) < 1).
Едно окончателно ограничение е длъжен да избере една от двете възможни точки, които са резултат от ограниченията, споменати по-горе.
Нека T(s) представлява "превод" операция, когато (s) е произволен вектор, така че T(s) * (r) = (r + s) и Inverse(T(s)) * (r) = (r - s).
Ротационен R(S,t,c), за равнина S, ъгъл t, както и център точка c, могат да бъдат изразени по отношение на превода операции и ротация с център точка по произход: R(S,t,c) = T(c) * R(S,t,0) * T(-c) = T(c) * R(S,t,0) * Inverse(T(c)).
Нека R(S,t) представлява редуване на операцията, с редуване на централната точка на произхода на (d) тримерно пространство.
По този начин, всеки произволен ротация може да се изрази като: R(S,t,c) = T(c) * R(S,t) * T(-c) = T(c) * R(S,t) * Inverse(T(c)).
Нека разгледаме двумерен ротация (S) равнини, които включват произхода на (d) тримерно пространство и са успоредно с две оси, да координира.
Нека (a) и (b) бъдат числа, които представляват координира оси (d) в тримерно пространство, така че 0 <= a <= (d-1) и 0 <= b <= (d-1) и (a != b).
Ние се дефинира R(a,b,t) да бъде ротационен принцип за произхода на (d) тримерно пространство, с редуване на координатната равнина успоредни оси, посочени от координира индекси (a) и (b), както и с ъгъл (t), така че да превърнат R(a,b,(pi/2)) положително координира стойност за координиране ос, посочена от (a) да координира положителна стойност за координатна ос, посочена от (b); Например: R(a,b,(pi/2)) * e(a) = e(b), където e(k) показва единица база вектори, съответстващи на координацията ос (k).
Това определение се дава на крайния многозначност е необходимо изцяло да се определи нова точка, произведено от въртене.
Ротационен, че се движи точки по оста a към оста b според ъгъл t е еквивалентно на въртене, че се движи точки по оста b към оста a с ъгъл (-t); т.е. R(a,b,t) = R(b,a,(-t)).
Ако (d >= 2), а после (d) тримерно пространство е ((d*(d-1))/2) чифта координира оси.
По този начин, за (d >= 2), (d) тримерно пространство е ((d*(d-1))/2) отделна ос-редуване на равнини, успоредни.
(Примери: d=2: () x-y равнина; d=3: (x-y равнина, y-z равнина, z-x равнина); d=4: (x-y равнина, y-z равнина, z-x равнина, x-w равнина, z-w равнина, w-y равнина).)
Кода по-долу показва как една матрица може да се използва за представляване на въртене R(a,b,t).
Един пример е видно, че използва ротация на матрицата в хомогенна начин да се сменят хомогенна вектор на нов хомогенна вектор.
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 Предложено определение на "авъртане"
В двумерните пространство, тримесечието завой ротационен принцип, който може да конвертира векторен e(0) да e(1) се счита авъртане.
Ротацията е R(0,1,(pi/2)) * e(0) = e(1).
Следователно, ротация R(0,1,t) е авъртане.