类文件
一共写了三个类
Matrix_Subtract.cs：矩阵减法
Matrix_Multiply.cs：矩阵乘法
QR.cs：QR算法
创建三个同名的类，依次把下面的代码Copy进去就可以直接用
QR.cs
/// <summary>
/// QR分解
/// </summary>
/// <param name="matrix"></param>
/// <returns></returns>
private static double[,] QRDecpt(ref double[,] matrix)
{
int rowNumber = matrix.GetLength(0);
int columnNumber = matrix.GetLength(1);
double x_norm;
double u_norm;
int count;
double[,] Q = new double[rowNumber, columnNumber];
for (int i = 0; i < rowNumber; i++)
{
Q[i, i] = 1;//先让Q矩阵为单位阵
}
for (int j = 0; j < columnNumber - 1; j++)
{
//先判断该列的主对角元下的元素是否为0，如果都为0就不需要进行化0步骤
int judgment = 0;
for (int i = j + 1; i < rowNumber; i++)
{
if (matrix[i, j] == 0)
{
judgment++;
}
}
if (judgment == rowNumber - j - 1)
{
continue;
}
else
{
double[,] x = new double[rowNumber - j, 1];
double[,] y = new double[rowNumber - j, 1];
double[,] u = new double[rowNumber - j, 1];
double[] uTrans = new double[rowNumber - j];
double[,] UnitMatrix = new double[rowNumber - j, rowNumber - j];
double[,] H = new double[rowNumber, columnNumber];
x_norm = 0;
u_norm = 0;
count = 0;
for (int i = j; i < rowNumber; i++)
{
x[count, 0] = matrix[i, j];
count++;
}
for (int i = 0; i < x.GetLength(0); i++)
{
x_norm += Math.Pow(x[i, 0], 2);//x的二范数
UnitMatrix[i, i] = 1;//单位阵
}
y[0, 0] = -Math.Pow(x_norm, 0.5);
u = Matrix_Subtract.Subtract(x, y);
for (int i = 0; i < u.GetLength(0); i++)
{
uTrans[i] = u[i, 0];
}
for (int i = 0; i < u.GetLength(0); i++)
{
u_norm += Math.Pow(u[i, 0], 2);//u的二范数
}
double[,] h = Matrix_Multiply.Multiply(u, uTrans);
for (int m = 0; m < h.GetLength(0); m++)
{
for (int n = 0; n < h.GetLength(1); n++)
{
h[m, n] = 2 * h[m, n] / u_norm;
}
}
h = Matrix_Subtract.Subtract(UnitMatrix, h);
if (j == 0)
{
H = h;
matrix = Matrix_Multiply.Multiply(H, matrix);
Q = Matrix_Multiply.Multiply(Q, H);
}
else
{
for (int i = 0; i < j; i++)
{
H[i, i] = 1;
}
for (int m = j; m < rowNumber; m++)
{
for (int n = j; n < columnNumber; n++)
{
H[m, n] = h[m - j, n - j];
}
}
matrix = Matrix_Multiply.Multiply(H, matrix);
Q = Matrix_Multiply.Multiply(Q, H);
}
}
}
return Q;
}
/// <summary>
/// 特征值计算
/// </summary>
/// <param name="lamda">A阵</param>
/// <returns></returns>
public static double[,] Eig(double[,] lamda)
{
double[,] Q;
int count = 0;
do
{
//开始迭代后lamda就是R矩阵
Q = QRDecpt(ref lamda);
lamda = Matrix_Multiply.Multiply(lamda, Q);
count++;
}
while (count < 50);//设定迭代50次
return lamda;
}
Matrix_Subtract.cs
public static class Matrix_Subtract
{
/// <summary>
/// 矩阵的减法
/// </summary>
/// <param name="matrix1">二维矩阵</param>
/// <param name="matrix2">二维矩阵</param>
/// <returns></returns>
public static double[,] Subtract(double[,] matrix1, double[,] matrix2)
{
int number = matrix1.GetLength(0);
//列向量相减
if (matrix1.GetLength(1) == 1 && matrix2.GetLength(1) == 1)
{
double[,] MatrixSubtract = new double[number, 1];
for (int i = 0; i < number; i++)
{
MatrixSubtract[i, 0] = matrix1[i, 0] - matrix2[i, 0];
}
return MatrixSubtract;
}
//方阵相减
else
{
double[,] MatrixSubtract = new double[number, number];
for (int i = 0; i < number; i++)
{
for (int j = 0; j < number; j++)
{
MatrixSubtract[i, j] = matrix1[i, j] - matrix2[i, j];
}
}
return MatrixSubtract;
}
}
}
Matrix_Multiply.cs
public static class Matrix_Multiply
{
/// <summary>
/// 矩阵乘法的运算
/// </summary>
/// <param name="matrix1">左矩阵</param>
/// <param name="matrix2">右矩阵</param>
/// <returns></returns>
public static double[,] Multiply(double[,] matrix1, double[,] matrix2)
{
int matrix1_row = matrix1.GetLength(0);
// int matrix1_column = matrix1.GetLength(1);
int matrix2_row = matrix2.GetLength(0);
int matrix2_column = matrix2.GetLength(1);
double[,] Matrix_Multiply = new double[matrix1_row, matrix2_column];
for (int i = 0; i < matrix1_row; i++)
{
for (int j = 0; j < matrix2_column; j++)
{
for (int m = 0; m < matrix2_row; m++)
{
Matrix_Multiply[i, j] += matrix1[i, m] * matrix2[m, j];
}
}
}
return Matrix_Multiply;
}
/// <summary>
/// 向量乘法的运算
/// </summary>
/// <param name="matrix1">列向量</param>
/// <param name="matrix2">行向量</param>
/// <returns></returns>
public static double[,] Multiply(double[,] matrix1, double[] matrix2)
{
int number = matrix2.Length;
double[,] matrix = new double[number, number];
for (int i = 0; i < number; i++)
{
for (int j = 0; j < number; j++)
{
matrix[i, j] = matrix1[i, 0] * matrix2[j];
}
}
return matrix;
}
}
算例验证
求以下矩阵的特征值
[ 1 2 3 2 1 3 3 3 6 ] \left[ \begin{matrix} 1 & 2 & 3 \\ 2 & 1 & 3 \\ 3 & 3 & 6 \end{matrix} \right]⎣⎡123213336⎦⎤
理论解为-1,0,9
Main函数里面写
class Program
{
static void Main(string[] args)
{
double[,] matrix1 = { { 1, 2, 3 }, { 2, 1, 3 }, { 3, 3, 6 } };
double[,] matrix;
matrix = QR.Eig(matrix1);
for (int i = 0; i < matrix.GetLength(0); i++)
{
Console.Write(matrix[i, i] + " ");
}
Console.ReadKey();
}
}
输出
广义特征值
表达式：A x = λ B x Ax=λBxAx=λBx，此时λ λλ即为广义特征值。可以把式子写成：B − 1 A x = λ x B^{-1}Ax=λxB−1Ax=λx，即转化为求B − 1 A B^{-1}AB−1A的特征值。
这样处理在有些场合不适用，比如计算力学里面求解临界力，固有频率等，有可能会出现B BB矩阵奇异的情况，一般可以这样处理：1 λ x = A − 1 B x \frac{1}{λ}x=A^{-1}Bxλ1x=A−1Bx，转化为求A − 1 B A^{-1}BA−1B的特征值再求倒数即可。
版权声明：本文为wangbo8366534原创文章，遵循CC 4.0 BY-SA版权协议，转载请附上原文出处链接和本声明。