矩阵和多项式代数
4.84/5 (33投票s)
本文展示了矩阵和多项式代数基本算法的实现
引言
对许多人来说，甚至对开发人员来说，线性代数相关的问题似乎是解决其他各种问题的障碍。为了应对这种情况，我们寻找现成的工具，这些工具我们不得不将其视为黑箱。本文的目的是展示最常见的矩阵和多项式代数案例的理论和代码实现。
矩阵
矩阵最常见的定义与其外观有关。我们将矩阵描述为存储在矩形表中的数字集合
$\begin{aligned} \begin{bmatrix} a_{1,1} & a_{1,2} & \cdots & a_{1,n} \\ a_{2,1} & a_{2,2} & \cdots & a_{2,n} \\ \vdots & \vdots & \ddots & \vdots \\ a_{m,1} & a_{m,2} & \cdots & a_{m,n} \end{bmatrix} \end{aligned}$
事实上，矩阵是定义在两个向量的笛卡尔积上的一个二元函数
$\begin{aligned} \begin{bmatrix} 1, ..., m \end{bmatrix} \times \begin{bmatrix} 1, ..., n \end{bmatrix} \end{aligned}$
因此，矩阵描述了单个向量元素之间的关系。
描述矩阵性质的要素有：迹、秩、行列式和特征多项式。矩阵的典型运算包括：转置、标量乘法、加法、减法、乘法和求逆。
Trace
矩阵的迹只对方阵有定义。它是矩阵对角线上所有元素的总和。
$\begin{aligned} \mathrm{trace}(A_{n \times n}) = \sum^{n}_{i=1} a_{i,i} \end{aligned}$
/// <summary>
/// Matrix trace
/// </summary>
/// <returns>Matrix trace</returns>
public Double Trace()
{
// int Rows, Columns - fileds of class Matrix, store number of rows and columns
// Double[,] coefs - private field of class Matrix, stores values of matrix elements
if (Rows == Columns)
{
Double Trace = 0;
for (int m = 0; m < Rows; m++)
{
Trace += coefs[m, m];
}
return Trace;
}
else
{
throw new Exception("Matrix must be square.");
}
}
此外，矩阵的迹将用于计算矩阵的特征多项式和逆矩阵。
秩
秩是矩阵最重要的特征之一。如果我们想象矩阵是由水平（行优先方法）或垂直（列优先方法）向量组成的，那么秩就回答了这些向量中有多少是线性无关的。
矩阵的秩与我们采取的方法无关，即矩阵A的行秩等于该矩阵的列秩。
秩可以分配给任何矩形矩阵。一般来说，对于维度为n x m的非零矩阵A，其秩的范围是
$\begin{aligned} 1 \le \mathrm{rank}(\textbf{A}) \le \mathrm{MIN}(m, n) \end{aligned}$
确定矩阵的秩等同于找到维度最大的子式（其维度对应于原始矩阵的秩）。因此，一个非零方阵A是满秩的，当且仅当其行列式不为零。
然而，子式法虽然在视觉上很优雅，但在计算大型矩阵的秩时效率不高。在这些情况下，我们可以使用高斯消元法。
这种方法的核心是利用对行（列）的线性变换，从矩阵A创建矩阵B。为了使矩阵A满秩，矩阵B应该是一个对角线以下的元素为零的三角矩阵，其中每行（或每列）所有元素的和都非零，因为A的秩是非零和的行（列）的数量。
设矩阵A的元素为
$\begin{aligned} \begin{bmatrix} a_{1,1} & a_{1,2} & \cdots & a_{1,n} \\ a_{2,1} & a_{2,2} & \cdots & a_{2,n} \\ \vdots & \vdots & \ddots & \vdots \\ a_{n,1} & a_{n,2} & \cdots & a_{n,n} \end{bmatrix} \end{aligned}$
第一步，我们将第一行的每个元素乘以
$\begin{aligned} c_{i,1} = \frac{a_{i,1}}{a_{1,1}} & & \mathrm{其中} & & i = 2, 3, ..., n \end{aligned}$
然后将乘法结果从第i行中减去。我们得到一个新的矩阵 A(1)
$\begin{aligned} \begin{bmatrix} a_{1,1} & a_{1,2} & \cdots & a_{1,n} \\ 0 & a_{2,2}^{(1)} & \cdots & a_{2,n}^{(1)} \\ \vdots & \vdots & \ddots & \vdots \\ 0 & a_{n,2}^{(1)} & \cdots & a_{n,n}^{(1)} \end{bmatrix} \end{aligned}$
实际上所做的是消除了ai,1元素（除了a1,1）。
第二步，我们将新第二行的每个元素乘以
$\begin{aligned} c_{i,2} = \frac{a_{i,2}}{a_{2,2}} & & \mathrm{其中} & & i = 3, 4, ..., n \end{aligned}$
然后将乘法结果从第i行中减去。我们得到新的矩阵 A(2)，消除了a(2)i,2元素（除了a(2)2,2）。
// coefs is a Double[,] structure that holds our matrix's coeficients
public int Rank()
{
int rank = 0;
Double[,] mat = coefs; // coefficients of our matrix
for (int k = 0; k < Rows; k++ )
{
for (int i = k + 1; i < Rows; i++ )
{
Double c = mat[i, k] / mat[k, k];
for (int k1 = 0; k1 < Columns; k1++ )
{
mat[i, k1] = mat[i, k1] - mat[k, k1] * c;
}
}
// Check if created row's elements sum is non-zero
Double sum = 0;
for (int i = 0; i < Columns; i++)
{
sum += mat[k, i];
}
if (sum != 0) { rank++; } // Increase rank if sum of new row is non-zero.
}
return rank;
}
请注意，在此算法的每一步k中，我们需要假设元素a(k)k,k非零（以便进行除法）。但是，如果我们找到元素a(k)k,k= 0，我们可以应用称为行主元法（或列主元法）的方法。我们选择一个第一行（从剩余的行中）具有非零a(k)l,l（其中l = k+1, k+2, ..., n）的行，并替换k行和l行（或列）的位置。
public int Rank()
{
int rank = 0;
Double[,] mat = coefs; // coefficients of our matrix
for (int k = 0; k < Rows; k++)
{
for (int i = k + 1; i < Rows; i++)
{
// Check if divider is equal to zero
if (mat[k,k] == 0)
{
// It is - pivot row
mat = RowPivot(mat, k);
}
Double c = mat[i, k] / mat[k, k];
for (int k1 = 0; k1 < Columns; k1++)
{
mat[i, k1] = mat[i, k1] - mat[k, k1] * c;
}
}
// Check if created row's elements sum is non-zero
Double sum = 0;
for (int i = 0; i < Columns; i++)
{
sum += mat[k, i];
}
if (sum != 0) { rank++; } // Increase rank if sum of new row is non-zero.
}
return rank;
}
private Double[,] RowPivot(Double[,] matrix, int k)
{
// k - starting row to search for non-zero element
for (int i = k + 1; i < matrix.GetLength(0); i++ )
{
if (matrix[i, i] != 0)
{
Double[] x = new Double[matrix.GetLength(1)];
for (int j = 0; j < matrix.GetLength(1); j++)
{
x[j] = matrix[k, j];
matrix[k, j] = matrix[i, j];
matrix[i, j] = x[j];
}
break;
}
}
return matrix;
}
为了快速检查我们的秩确定算法，让我们考虑以下示例，其中矩阵A由线性无关的行构成，矩阵B的第1行和第3行是线性相关的，而矩阵C的第一行以0开头
$\begin{aligned} (\mathrm{Ex. 1}) & A = \begin{bmatrix} 1 & 2 & 7 \\ 3 & 4 & 1 \\ -1 & 3 & 1 \end{bmatrix} & (\mathrm{Ex. 2}) & B = \begin{bmatrix} 1 & 2 & 7 \\ 3 & 4 & 1 \\ 2 & 4 & 14 \end{bmatrix} & (\mathrm{Ex. 3}) & C = \begin{bmatrix} 0 & 2 & 7 \\ 3 & 4 & 1 \\ -1 & 3 & 1 \end{bmatrix} \end{aligned}$
根据我们的算法计算，示例1中的矩阵是
$\begin{aligned} A^{(0)} = \begin{bmatrix} 1 & 2 & 7 \\ 3 & 4 & 1 \\ -1 & 3 & 1 \end{bmatrix} \end{aligned}$
$\begin{aligned} A^{(1)} = \begin{bmatrix} 1 & 2 & 7 \\ 0 & -2 & -20 \\ 0 & 5 & 8 \end{bmatrix} \end{aligned}$
$\begin{aligned} A^{(2)} = \begin{bmatrix} 1 & 2 & 7 \\ 0 & -2 & -20 \\ 0 & 0 & -42 \end{bmatrix} \end{aligned}$
每行的所有元素之和非零，因此矩阵A的秩为3（满秩）。
根据我们的算法计算，示例2中的矩阵是
$\begin{aligned} B^{(0)} = \begin{bmatrix} 1 & 2 & 7 \\ 3 & 4 & 1 \\ 2 & 4 & 14 \end{bmatrix} \end{aligned}$
$\begin{aligned} B^{(1)} = \begin{bmatrix} 1 & 2 & 7 \\ 0 & -2 & -20 \\ 0 & 0 & 0 \end{bmatrix} \end{aligned}$
所有元素之和非零的行数为2。矩阵B的秩为2。
根据我们的算法计算，示例3中的矩阵是
步骤0：算法发现元素c1,1为零，因此应用了行主元法
$\begin{aligned} C^{(0)} = \begin{bmatrix} 3 & 4 & 1 \\ 0 & 2 & 7 \\ -1 & 3 & 1 \end{bmatrix} \end{aligned}$
$\begin{aligned} C^{(1)} = \begin{bmatrix} 3 & 4 & 1 \\ 0 & 2 & 7 \\ 0 & 4.33 & 1.33 \end{bmatrix} \end{aligned}$
$\begin{aligned} C^{(2)} = \begin{bmatrix} 3 & 4 & 1 \\ 0 & -2 & -7 \\ 0 & 0 & -13.83 \end{bmatrix} \end{aligned}$
每行的所有元素之和非零，因此矩阵C的秩为3（满秩）。
行列式
根据形式定义，矩阵行列式是一个函数，它将每个方阵An x n{R}（其元素来自交换环R）映射到该交换环R中的一个元素。
为了我们的目的，我们将矩阵行列式定义为一个描述矩阵（实际上是描述矩阵元素）的另一个函数。
在符号数学中，行列式使用Leibniz公式（也称为行列式的排列定义）计算
$\begin{aligned} \mathrm{det}A = \sum_{s \in S_n} -1^{In(s)} \prod_{i=1}^{n} a_{i, s_{i}} \end{aligned}$
其中Sn是集合{1, 2, ..., n}的所有可能排列的集合，In(s)是排列s的可能反转次数。
虽然上述定义看起来很复杂，但我们习惯用它来计算2x2矩阵的行列式。我们有
$\begin{aligned} \mathrm{det} A = \mathrm{det}\begin{bmatrix} a_{1,1} & a_{1,2}\\ a_{2,1} & a_{2,2} \end{bmatrix} = a_{1,1} a_{2,2} - a_{1,2}a_{2,1} \end{aligned}$
在计算3x3矩阵的行列式时，所有简化都已由Pierre Sarrus完成
$\begin{aligned} \begin{matrix} \mathrm{det} A = \mathrm{det}\begin{bmatrix} a_{1,1} & a_{1,2} & a_{1,3} \\ a_{2,1} & a_{2,2} & a_{2,3} \\ a_{3,1} & a_{3,2} & a_{3,3} \end{bmatrix} = \\ = (a_{1,1}a_{2,2}a_{3,3} + a_{1,2}a_{2,3}a_{3,1} + a_{1,3}a_{2,1}a_{3,2}) - (a_{1,3}a_{2,2}a_{3,1} + a_{1,1}a_{2,3}a_{3,2} + a_{1,2}a_{2,1}a_{3,3}) \end{matrix} \end{aligned}$
在离散数学中，我们通常使用拉普拉斯展开递归算法，其中
$\begin{aligned} \mathrm{det} A_{n \times n} = \sum_{j=1}^n a_{i,j} A_{i,j} \end{aligned}$
其中i是我们进行展开的固定行号，ai,j是坐标为i,j的矩阵元素，而Ai,j是元素ai,j的代数余子式。
上述公式递归使用，直到代数余子式维度为2x2（此时可以直接通过定义计算行列式）。元素ai,j的代数余子式定义为
$\begin{aligned} A_{i,j} = -1^{i+j} M_{i,j} \end{aligned}$
其中Mi,j是矩阵A的子式，通过删除第i行和第j列创建。
下面的代码与附加库中的代码不同，以便更好地展示算法的工作原理。
public Double MatrixDeterminant(Double[,] Matrix)
{
int M = Matrix.GetLength(0);
int N = Matrix.GetLength(1);
if (M == N)
{
Double Det = 0;
if (M == 1)
{
Det = Matrix[0, 0];
}
else if (M == 2)
{
Det = Matrix[0, 0] * Matrix[1, 1] - Matrix[0, 1] * Matrix[1, 0]; // by definition
}
else
{
Double[,] AlgebraicComplement = new Double[M - 1, M - 1];
for (int m = 0; m < M; m++)
{
int a = 0;
for (int k = 1; k < M; k++)
{
int b = 0;
for (int l = 0; l < M; l++)
{
if (l != m)
{
AlgebraicComplement[a, b] = Matrix[k, l];
b++;
}
}
a++;
}
Det += Math.Pow(-1, m) * Matrix[0, m] * MatrixDeterminant(AlgebraicComplement);
}
}
return Det;
}
else
{
throw new Exception("Matrix must be square.");
}
}
考虑矩阵A（上一部分的示例1）
$\begin{aligned} A = \begin{bmatrix} 1 & 2 & 7 \\ 3 & 4 & 1 \\ -1 & 3 & 1 \end{bmatrix} \end{aligned}$
计算出的行列式为84。但让我们回顾一下高斯消元法之后该矩阵的结构
$\begin{aligned} A^{(2)} = \begin{bmatrix} 1 & 2 & 7 \\ 0 & -2 & -20 \\ 0 & 0 & -42 \end{bmatrix} \end{aligned}$
如果我们将其对角线上的所有元素相乘，也会得到84。这是三角矩阵的一个性质——行列式等于对角线元素之积。因此，我们也可以使用高斯消元法来计算矩阵的行列式。
特征多项式
方阵An x n的特征多项式定义为
$\begin{aligned} p(x) = \mathrm{det} (A - x \textbf{1}) \end{aligned}$
其中1是单位矩阵（对角线上为1）。
如果我们使用多项式记法，我们可以将p(x)重写为
$\begin{aligned} p(x) = x^n + p_1 x^{n-1} + ... + p_{n-1} x + p_n \end{aligned}$
为了计算p的系数，我们可以使用Faddeev - Leverrier算法，在该算法中，我们每一步都“减去”上一步计算出的系数
$\begin{aligned} \begin{matrix} A^{(1)} = A & \to & p_1 = \frac{1}{1} \mathrm{trace} A^{(1)} \\ A^{(2)} = A (A^{(1)} - p_1 \textbf{1}) & \to & p_2 = \frac{1}{2} \mathrm{trace} A^{(2)} \\ \vdots & &\vdots \\ A^{(n)} = A (A^{(n-1)} - p_{n-1} \textbf{1}) & \to & p_n = \frac{1}{n} \mathrm{trace} A^{(n)} \end{matrix} \end{aligned}$
下面的实现与附加库中的代码不同。
/// <summary>
/// The method calculates characteristic polynominal coefficients of given matrix
/// </summary>
/// <param name="Matrix">Real numbers matrix</param>
/// <returns>Characteristic polynominal coefficients (vector starting with free term)</returns>
public Double[] CharacteristicPolynominalCoefficients(Double[,] Matrix)
{
int M = Matrix.GetLength(0);
int N = Matrix.GetLength(1);
if (M == N)
{
Double[] Coeffs = new Double[M + 1];
Double[] CoeffsSorted = new Double[M + 1];
Double[,] B;
Double LastCoeff;
Coeffs[0] = 1;
B = Matrix;
LastCoeff = MatrixTrace(B);
Coeffs[1] = -LastCoeff;
for (int m = 2; m < M + 1; m++)
{
B = MatrixMultiplication(Matrix, MatrixSubtraction(B, MatrixScaling(LastCoeff, IdentityMatrix(B.GetLength(0)))));
LastCoeff = MatrixTrace(B) / m;
Coeffs[m] = -LastCoeff;
}
for (int m = 0; m < M + 1; m++)
{
CoeffsSorted[m] = Coeffs[M - m];
}
return CoeffsSorted;
}
else
{
throw new Exception("Input data matrix must be square.");
}
}
函数p(x)的根称为矩阵A的特征值，在下一部分将展示计算多项式根的算法。
转置
矩阵转置运算会创建一个新矩阵，该矩阵是原始矩阵在其主对角线上的镜像。换句话说，转置函数交换行和列。
$\begin{aligned} B = A^{T} & \to & \bigwedge_{i} \bigwedge_{j} b_{i,j} = a_{j,i} \end{aligned}$
/// <summary>
/// Transpose matrix
/// </summary>
/// <returns>Transposed matrix</returns>
public Matrix Transpose()
{
// int Rows, Columns - fileds of class Matrix, store number of rows and columns
// Double[,] coefs - private field of class Matrix, stores values of matrix elements
Double[,] _MatrixT = new Double[Columns, Rows];
for (int m = 0; m < Rows; m++)
{
for (int n = 0; n < Columns; n++)
{
_MatrixT[n, m] = coefs[m, n];
}
}
return new Matrix(_MatrixT);
}
扩展
矩阵的标量乘法将每个矩阵元素乘以给定的标量。
$\begin{aligned} B = s \cdot A & \to & \bigwedge_{i} \bigwedge_{j} b_{i,j} = s \cdot a_{i,j} \end{aligned}$
/// <summary>
/// Multiply each matrix element by scalar
/// </summary>
/// <param name="Coefficient">Multiplier scalar</param>
/// <returns>Scalerd matrix</returns>
public Matrix Scale(Double Coefficient)
{
// int Rows, Columns - fileds of class Matrix, store number of rows and columns
// Constructor Matrix(int r, int c, Double v), creates new matrix in dimension r,c filled with
// values v
// method Element(int r, int c) retuirn value to element in position r,c
// method SetElement(int r, int c, Double v) set element r,c value to v
Matrix result = new Matrix(Rows, Columns, 0);
for (int m = 0; m < Rows; m++)
{
for (int n = 0; n < Columns; n++)
{
result.SetElement(m, n, this.Element(m, n) * Coefficient);
}
}
return result;
}
加法和减法
矩阵的加法和减法运算仅对相同维度的矩阵定义，并按如下方式计算：
$\begin{aligned} C_{m \times n} = A_{m \times n} \pm B_{m \times n} & \to & \bigwedge_{i} \bigwedge_{j} c_{i,j} = a_{i,j} \pm b_{i,j} \end{aligned}$
/// <summary>
/// Adds new matrix to existing one
/// </summary>
/// <param name="MatrixB">Real numbers matrix</param>
/// <returns>Sum</returns>
public Matrix Add(Matrix MatrixB)
{
// int Rows, Columns - fileds of class Matrix, store number of rows and columns
// method Element(int r, int c) retuirn value to element in position r,c
// method SetElement(int r, int c, Double v) set element r,c value to v
if (Rows == MatrixB.Rows && Columns == MatrixB.Columns)
{
Matrix result = new Matrix(Rows, Columns, 0);
for (int n = 0; n < Rows; n++)
{
for (int m = 0; m < Columns; m++)
{
result.SetElement(n, m, (this.Element(n, m) + MatrixB.Element(n, m)));
}
}
return result;
}
else
{
throw new Exception("Matrices dimensions must be equal.");
}
}
/// <summary>
/// Subtracts new matrix from existing one
/// </summary>
/// <param name="MatrixB">Real numbers matrix</param>
/// <returns>Result of subtraction</returns>
public Matrix Sub(Matrix MatrixB)
{
if (Rows == MatrixB.Rows && Columns == MatrixB.Columns)
{
Matrix result = new Matrix(Rows, Columns, 0);
for (int n = 0; n < Rows; n++)
{
for (int m = 0; m < Columns; m++)
{
result.SetElement(n, m, (this.Element(n, m) - MatrixB.Element(n, m)));
}
}
return result;
}
else
{
throw new Exception("Matrices dimensions must be equal.");
}
}
乘法
矩阵乘法运算使用柯西乘法公式定义如下：
$\begin{aligned} C_{m \times p} = A_{m \times n} \cdot B_{n \times p} & \to & \bigwedge_{i: 1 \le i \le m} \bigwedge_{j: 1 \le j \le p} c_{i,j} = \sum_{k=1}^n a_{i,k} b_{k,j} \end{aligned}$
只有当矩阵A的列数等于矩阵B的行数时，才能应用该运算。
/// <summary>
/// Multiplies (right-side) existing matrix (multiplicand) by given one (multiplier)
/// </summary>
/// <param name="MatrixB">Multiplier matrix</param>
/// <returns>Result of multiplification</returns>
public Matrix Multiply(Matrix MatrixB)
{
if (Columns != MatrixB.Rows)
{
throw new Exception("Number of columns in A matrix must be equal to number of rows in B matrix.");
}
else
{
Double[,] _MatrixAB = new Double[Rows, MatrixB.Columns];
for (int m_a = 0; m_a < Rows; m_a++)
{
for (int n_b = 0; n_b < MatrixB.Columns; n_b++)
{
_MatrixAB[m_a, n_b] = 0;
for (int n_a = 0; n_a < MatrixB.Rows; n_a++)
{
_MatrixAB[m_a, n_b] = _MatrixAB[m_a, n_b] + this.Element(m_a, n_a) * MatrixB.Element(n_a, n_b);
}
}
}
return new Matrix(_MatrixAB);
}
}
求逆
根据定义，矩阵An x n是可逆的，当且仅当存在矩阵B满足方程
$\begin{aligned} A B = B A = \textbf{1}_n \end{aligned}$
矩阵B称为矩阵A的逆，通常表示为A-1。并非所有矩阵都可逆，为了使矩阵A可逆，它必须是非奇异的——其行列式必须非零。
我们可以找到许多计算逆矩阵的算法，但一旦我们计算了Faddeev - Leverrier方法，我们将再次使用它，因为矩阵A的逆可以使用以下公式确定：
$\begin{aligned} A^{-1} = \frac{1}{p_n} (A^{(n-1)} - p_{n-1} \textbf{1}) \end{aligned}$
其中A(n-1)是在算法的n-1步中确定的矩阵，pn和pn-1是在n-1和n步中确定的特征多项式系数。
/// <summary>
/// Calculates the inverse matrix using Faddeev-Leverrier algorithm
/// </summary>
/// <returns>Inverse matrix</returns>
public Matrix Inverse()
{
if (Rows == Columns)
{
if (this.Determinant() != 0)
{
Matrix LastB;
Matrix NextB;
Matrix Inv = new Matrix(Rows, Columns, 0);
Double LastCoeff;
Double NextCoeff;
LastB = this;
LastCoeff = LastB.Trace();
for (int m = 2; m < Rows; m++)
{
LastB = this.Multiply(LastB.Sub(new Matrix(Rows, 1).Scale(LastCoeff)));
LastCoeff = LastB.Trace() / m;
}
NextB = this.Multiply(LastB.Sub(new Matrix(Rows, 1).Scale(LastCoeff)));
NextCoeff = NextB.Trace() / Rows;
Inv = LastB.Sub(new Matrix(Rows, 1).Scale(LastCoeff)).Scale(1 / NextCoeff);
return Inv;
}
else
{
throw new Exception("The matrix is not invertible over commutative ring.");
}
}
else
{
throw new Exception("Input data matrix must be square.");
}
}
例如，我们将再次使用矩阵
$\begin{aligned} A = \begin{bmatrix} 1 & 2 & 7 \\ 3 & 4 & 1 \\ -1 & 3 & 1 \end{bmatrix} \end{aligned}$
计算出的逆矩阵是
$\begin{aligned} A^{-1} = \begin{bmatrix} 0.0119047619047619 & 0.226190476190476 & -0.30952380952381 \\ -0.0476190476190476 & 0.0952380952380952 & 0.238095238095238 \\ 0.154761904761905 & -0.0595238095238095 & -0.0238095238095238 \end{bmatrix} \end{aligned}$
根据定义，A A-1或A-1 A的乘积应该得到单位矩阵。我们实际得到的是
$\begin{aligned} A A^{-1} = \begin{bmatrix} 1 & 0 & -2.77555756156289E-17 \\ 0 & 1 & -8.32667268468867E-17 \\ 0 & 1.38777878078145E-17 & 1 \end{bmatrix} \end{aligned}$
多项式
单变量多项式是一种数学表示法（也称为代数和），表示单项式的和。
$\begin{aligned} \sum_{i=0}^n a_i x^i \end{aligned}$
其中i, n是整数，ai和x可以是复数。上述可以写成
$\begin{aligned} a_n x^n + a_{n-1} x^{n-1} + ... + a_1 x + a_0 \end{aligned}$
上述多项式的次数自然是n。当x是一组值时，我们处理的是多项式函数
$\begin{aligned} f(x) = a_n x^n + a_{n-1} x^{n-1} + ... + a_1 x + a_0 \end{aligned}$
评估版
多项式求值（或更精确地说，多项式函数求值）是计算给定x值的f(x)值。
/// <summary>
/// The method calculates polynominal evaluation with given value
/// </summary>
/// <param name="Value">Real value</param>
/// <returns>Value</returns>
public dynamic Evaluate(Complex Value)
{
Complex sum = 0;
// Complex[] coefs filed stores polynominal coefficients
for (int i = 0; i < coefs.Length; i++)
{
sum += coefs[i] * Complex.Pow(Value, i);
}
if (sum.Imaginary == 0)
{
// if result is real we return Double
return sum.Real;
}
else
{
// otherwise, we return Complex
return sum;
}
}
加法和减法
将两个（具有相同变量的）多项式相加（相减）会产生一个新多项式，其次数是相加（相减）分量中最大次数的次数。对应于x相同次数的系数可以相加（相减）。
/// <summary>
/// The method returns product of polynominal addition
/// </summary>
/// <param name="poly">Polynominal</param>
/// <returns>Polynominal</returns>
public Polynominal Add(Polynominal poly)
{
int M = this.Degree();
int N = poly.Degree();
int K = Math.Max(M, N);
Polynominal Poly1Ext = new Polynominal(K);
Polynominal Poly2Ext = new Polynominal(K);
Polynominal Add = new Polynominal(K);
for (int m = 0; m < M + 1; m++)
{
Poly1Ext.SetElement(m, this.Element(m));
}
for (int n = 0; n < N + 1; n++)
{
Poly2Ext.SetElement(n, poly.Element(n));
}
for (int k = 0; k < K + 1; k++)
{
Add.SetElement(k, Poly1Ext.Element(k) + Poly2Ext.Element(k));
}
return Add;
}
扩展
多项式标量乘法是定义的操作
$\begin{aligned} c \cdot (a_n x^n + a_{n-1} x^{n-1} + ... + a_1 x + a_0) = c \cdot a_n x^n + c \cdot a_{n-1} x^{n-1} + ... + c \cdot a_1 x + c \cdot a_0 \end{aligned}$
/// <summary>
/// The method returns product of polynominal multiplication by scalar
/// </summary>
/// <param name="Value">Complex value</param>
/// <returns>Polynominal</returns>
public Polynominal Scale(Complex Value)
{
Polynominal scale = new Polynominal(this.Degree());
for (int i = 0; i < this.Degree() + 1; i++)
{
scale.SetElement(i, this.Element(i) * Value);
}
return scale;
}
多项式除以二项式
用二项式除以多项式是通过消除多项式的一个已知根来降低次数的常用方法——我们将使用它来描述确定多项式函数根的算法。
为了用二项式除以多项式，我们将使用Horner方法。
设p(x)为n次多项式
$\begin{aligned} p(x) = a_n x^n + a_{n-1} x^{n-1} + ... + a_1 x + a_0 \end{aligned}$
设b(x)为二项式
$\begin{aligned} d(x) = x - c \end{aligned}$
我们可以写成
$\begin{aligned} p(x) = d(x) b(x) + r \end{aligned}$
其中r是除法的余数，而b(x)是n-1次多项式。
我们有
$\begin{aligned} a_n x^n + a_{n-1} x^{n-1} + ... + a_1 x + a_0 = (x - c)(b_{n-1} x^{n-1} + ... + b_1 x + b_0) + r \end{aligned}$
通过将方程右侧相乘并比较左侧的系数，我们可以看到
$\begin{aligned} \begin{matrix} b_{n-1} = a_n \\ b_{n-2} = a_{n-1} + c b_{n-2} \\ \vdots \\ b_{n-k} = a_{n-k+1} + c b_{n-k} \\ \vdots \\ b_0 = a_1 + c b_1 \\ r = a_0 + c b_0 \end{matrix} \end{aligned}$
/// <summary>
/// The method returns product of polynominal by binominal division
/// </summary>
/// <param name="BinominalRoot">Binominal complex root</param>
/// <returns>Polynominal</returns>
public Polynominal ByBinominalDivision(Complex BinominalRoot)
{
int M = this.Degree() + 1;
if (M > 1)
{
int N = M - 1;
Complex[] Quotient = new Complex[N];
Complex[] QuotientSorted = new Complex[N];
Complex[] CoeffsSorted = new Complex[M];
Complex Last;
for (int m = 0; m < M; m++)
{
CoeffsSorted[m] = this.Element(M - m - 1);
}
Last = CoeffsSorted[0];
Quotient[0] = Last;
for (int m = 1; m < N; m++)
{
Last = CoeffsSorted[m] + BinominalRoot * Last;
Quotient[m] = Last;
}
Complex Remainder = CoeffsSorted[M - 1] + BinominalRoot * Last; // the method is not returning Reminder!
for (int n = 0; n < N; n++)
{
QuotientSorted[n] = Quotient[N - n - 1];
}
return new Polynominal(QuotientSorted);
}
else
{
throw new Exception("Given coefficients number is not corresponding to correct polynominal structure.");
}
}
导数
多项式的导数可以计算为单项式导数的和，这是一个基本运算。
$\begin{aligned} \frac{d}{dx} a_n x^n = a_n \frac{d}{dx} x^n = a_n n x^{n-1} \end{aligned}$
/// <summary>
/// The method calculates polynominal derivative
/// </summary>
/// <returns>Polynominal</returns>
public Polynominal Derivative()
{
int Degree = this.Degree(); // degree of our polynominal
int NumCoefs = Degree + 1;
Polynominal Derivative = new Polynominal(this.Degree()-1); // new polynominal
if (Degree > 0)
{
for (int i = 1; i < NumCoefs; i++)
{
Derivative.SetElement(i - 1, i * this.Element(i));
}
}
else
{
return new Polynominal(new Double[1] { 0 });
}
return Derivative;
}
根
为了找到多项式函数的根，我们将使用Laguerre方法，该方法对重根的排列不敏感（与Horner方法相反）。Laguerre方法基于每次找到一个根时（使用二项式除法）来降低多项式次数。找到的根也会使用原始多项式进行平滑。
递归地，我们使用以下公式计算每个根的下一个近似值：
$\begin{aligned} r_{i+1} = r_i - \frac{n p(r_i)}{ \frac{d}{dx}p(r_i) \pm \sqrt{(n-1)( (n-1)(\frac{d}{dx}p(r_i))^2 - n p(r_i) \frac{d^2}{dx^2}p(r_i) ) } } \end{aligned}$
其中z0是根的第一个近似值——我们可以使用任何随机数。
分母中的符号应选择为最大化其模。
/// <summary>
/// The method calculates polynominal roots
/// </summary>
/// <param name="Seed">Seed root (if unknown any complex number can be given)</param>
/// <param name="Accuracy">Result accuracy</param>
/// <returns>Polynominal roots</returns>
public Complex[] Roots(Complex Seed, Double Accuracy)
{
int N = this.Degree();
Complex[] Roots = new Complex[N];
int degree = N;
Polynominal tmpoly = this;
for (int n = 0; n < N; n++)
{
while (true)
{
Complex _tmp0 = tmpoly.Derivative().Evaluate(Seed);
Complex _tmp1 = (degree - 1) * Complex.Pow(tmpoly.Derivative().Evaluate(Seed), 2);
Complex _tmp2 = degree * tmpoly.Evaluate(Seed) * tmpoly.Derivative().Derivative().Evaluate(Seed);
Complex _tmp3 = Complex.Pow((degree - 1) * (_tmp1 - _tmp2), 0.5);
Complex _tmp4 = MaximizeResultAbsolute(_tmp0, _tmp3);
Seed = Seed - (degree * tmpoly.Evaluate(Seed)) / _tmp4;
Complex result = tmpoly.Evaluate(Seed);
if (Complex.Abs(result) < Math.Abs(Accuracy))
{
Roots[n] = Seed;
tmpoly = tmpoly.ByBinominalDivision(Seed);
degree--;
Random _Random = new Random();
Seed = -10 + (_Random.NextDouble() * 15);
break;
}
}
}
return Roots;
}
private Complex MaximizeResultAbsolute(Complex Value1, Complex Value2)
{
if (Complex.Abs(Value1 + Value2) > Complex.Abs(Value1 - Value2))
{
return Value1 + Value2;
}
else
{
return Value1 - Value2;
}
}
例如，设p(x)为
$\begin{aligned} p(x) = x^3 + 3x^2 -3x -3 \end{aligned}$
我们可以计算出p(x)有三个实根：1.2618022452599722、-3.601679131883154和-0.660123113376818。
使用我们的算法，以种子2和精度0.0001，我们得到以下根（复数）：
x1 = 1.26180224809477 + 0 i
x2 = -3.60167913003098 - 1.46390615026496E-15 i
x3 = -0.660123118063785 + 5.44361994589775E-16 i
使用代码
本文中介绍的所有算法都已在附加的源代码中实现。要开始使用该库，只需声明其命名空间。
using MatrixPolynominal;
创建新矩阵和多项式
Matrix c = new Matrix(new Double[3, 3] { { 1, 2, 7 }, { 3, 4, 1 }, { -1, 3, 1 } });
Polynominal p = new Polynominal(new Double[4] {-3, -3, 3, 1});
每次我们想在控制台打印矩阵或多项式时，只需调用Print方法。
c.Print(); p.Print();
矩阵运算
// new matrix from array
Matrix c = new Matrix(new Double[3, 3] { { 1, 2, 7 }, { 3, 4, 1 }, { -1, 3, 1 } });
// new matrix with 3 rows, 3 columns, all elements equal 10
Matrix c1 = new Matrix(3, 3, 10);
// new matrix 3 x 3, all diagonal values equal 1
Matrix c2 = new Matrix(3, 1);
Matrix result;
// matrix transposition
result = c.Transpose();
// matrix scaling
result = c.Scale(10);
// multiply by matrix
result = c.Multiply(new Matrix(new Double[3, 3] { { 12, 12, 17 }, { -3, -4, -1 }, { 0, 0, 1 } }));
// inverse matrix
result = c.Inverse();
// add, sub matrix
result = c.Add(new Matrix(new Double[3, 3] { { 12, 12, 17 }, { -3, -4, -1 }, { 0, 0, 1 } }));
result = c.Sub(new Matrix(new Double[3, 3] { { 12, 12, 17 }, { -3, -4, -1 }, { 0, 0, 1 } }));
// get matrix rank
int rank = c.Rank();
// get determinant
Double det = c.Determinant();
// characteristic polynominal
Polynominal cp = c.CharacteristicPolynominal();
// Get matrix element
Double el = c.Element(1, 1);
// Set matrix element
c.SetElement(1,1, el);
多项式运算
// new polynominal from array (Double or Complex)
Polynominal p = new Polynominal(new Double[4] {-3, -3, 3, 1});
Polynominal p1 = new Polynominal(new Double[4] { 1, 2, 3, 4 });
Polynominal result;
// add, sub polynominals
result = p.Add(p1);
result = p.Sub(p1);
// polynominal derivative
result = p.Derivative();
// calculate division by binominal (x - 10)
result = p.ByBinominalDivision(10);
// evaluate with value 10
Complex value = p.Evaluate(10);
// get degree
int degree = p.Degree();
// scale with value 4
result = p.Scale(4);
// get element on position 2
Complex el = p.Element(2);
// set element on position 2
p.SetElement(2, el);
// calculate roots with seed 2 and accuracy 0.0001
Complex[] x = p.Roots(2, 0.0001);