[Serializable]

public class Matrix

{

/// <summary>

/// 提供矩阵求逆、特殊变换和特殊矩阵生成的方法

/// </summary>

int rows;

int columns;

double[,] value;

/// <summary>

/// 获取矩阵行数

/// </summary>

public int Rows

{

get

{

return rows;

}

set

{

rows = value;

}

}

/// <summary>

/// 获取矩阵列数

/// </summary>

public int Columns

{

get

{

return columns;

}

set

{

columns = value;

}

}

/// <summary>

/// 获取与设置矩阵元素值的二维数组

/// </summary>

public double[,] Value

{

get

{

return this.value;

}

set

{

this.value = value;

}

}

/// <summary>

/// 获取与设置矩阵指定行列元素

/// </summary>

/// <param name="i">行号</param>

/// <param name="j">列号</param>

/// <returns></returns>

public double this[int i, int j]//索引器

{

set

{

this.value[i, j] = value;

}

get

{

return this.value[i, j];

}

}

#region 构造函数

private Matrix()//这是一个快速创建什么字段都没有初始化的矩阵对象的方法，只允许内部使用！

{

}

private Matrix(int rows, int cols)

{

this.rows = rows;

this.columns = cols;

this.value = new double[rows, cols];

}

/// <summary>

/// 构造一个一行一列矩阵

/// </summary>

/// <param name="num"></param>

public Matrix(double num)

{

rows = 1;

columns = 1;

value = new double[1, 1] { {num} };

}

/// <summary>

/// 构造一个行向量

/// </summary>

/// <param name="num"></param>

public Matrix(double[] num)

{

rows = 1;

columns = num.GetLength(0);

value = new double[1,columns];

for (int i = 0; i < columns; i++)

{

value[0, i] = num[i];

}

}

/// <summary>

/// 以二维数组初始化一个矩阵

/// </summary>

/// <param name="num"></param>

public Matrix(double[,] num)

{

rows = num.GetLength(0);

columns = num.GetLength(1);

value = new double[rows, columns];

for (int i = 0; i < rows; i++)

for (int j = 0; j < columns; j++)

value[i, j] = num[i, j];

}

/// <summary>

/// 以指定矩阵初始化一个矩阵(产生一个内容相同引用不同的矩阵)

/// </summary>

/// <param name="inMatrix"></param>

public Matrix(Matrix inMatrix)

{

rows = inMatrix.rows;

columns = inMatrix.columns;

value = new double[rows, columns];

for (int i = 0; i < rows; i++)

for (int j = 0; j < columns; j++)

value[i, j] = inMatrix.value[i, j];

}

#endregion

#region 特殊矩阵生成

/// <summary>

/// 生成一个指定阶的元素全为1的方阵

/// </summary>

/// <param name="dimension"></param>

/// <returns></returns>

public static Matrix Ones(int dimension)

{

Matrix result = new Matrix();

result.rows = dimension;

result.columns = dimension;

result.value = new double[dimension, dimension];

for (int i = 0; i < dimension; i++)

for (int j = 0; j < dimension; j++)

result.value[i, j] = 1;

return result;

}

/// <summary>

/// 生成一个指定行列的元素全为1的矩阵

/// </summary>

/// <param name="row"></param>

/// <param name="column"></param>

/// <returns></returns>

public static Matrix Ones(int row, int column)

{

Matrix result = new Matrix();

result.rows = row;

result.columns = column;

result.value = new double[row, column];

for (int i = 0; i < row; i++)

for (int j = 0; j < column; j++)

result.value[i, j] = 1;

return result;

}

/// <summary>

/// 生成一个指定阶数的单位矩阵

/// </summary>

/// <param name="dimension"></param>

/// <returns></returns>

public static Matrix Eye(int dimension)

{

Matrix result = new Matrix();

result.rows = dimension;

result.columns = dimension;

result.value = new double[dimension, dimension];

for (int i = 0; i < dimension; i++)

for (int j = 0; j < dimension; j++)

{

if (i == j)

result.value[i, j] = 1;

else

result.value[i, j] = 0;

}

return result;

}

/// <summary>

/// 生成一个互换初等矩阵

/// </summary>

/// <param name="dimension"></param>

/// <param name="i"></param>

/// <param name="j"></param>

/// <returns></returns>

public static Matrix ElementarySwitch(int dimension, int i, int j)//互换初等矩阵,i、j行(列)互换

{

if (i >= dimension)

throw new Exception("行列互换初等矩阵互换的行列号必须为0到dimension-1");

if (j >= dimension)

throw new Exception("行列互换初等矩阵互换的行列号必须为0到dimension-1");

Matrix result = Matrix.Eye(dimension);

if (i!=j)

{

result.value[i, i] = 0;

result.value[i, j] = 1;

result.value[j, i] = 1;

result.value[j, j] = 0;

}

return result;

}

/// <summary>

/// 生成一个倍乘初等矩阵

/// </summary>

/// <param name="dimension"></param>

/// <param name="i"></param>

/// <param name="k"></param>

/// <returns></returns>

public static Matrix ElementaryMultiple(int dimension, int i, double k)//倍乘初等矩阵，k倍i行(列)

{

if (i >= dimension)

throw new Exception("倍乘初等矩阵的行列号必须为0到dimension-1");

Matrix result = Matrix.Eye(dimension);

result.value[i, i] = k * result.value[i, i];

return result;

}

/// <summary>

/// 生成一个倍加初等矩阵

/// </summary>

/// <param name="dimension"></param>

/// <param name="i"></param>

/// <param name="k"></param>

/// <param name="j"></param>

/// <returns></returns>

public static Matrix ElementaryMulAdd(int dimension, int i, double k, int j)//倍加初等矩阵,i行加k倍的j行(j列加k倍i列)

{

if (i >= dimension)

throw new Exception("行列倍加初等矩阵倍加的行列号必须为0到dimension-1");

if (j >= dimension)

throw new Exception("行列倍加初等矩阵倍加的行列号必须为0到dimension-1");

Matrix result = Matrix.Eye(dimension);

if (i!=j)

{

result.value[i, j] = k;

}

return result;

}

/// <summary>

/// 生成一个指定阶数的零矩阵

/// </summary>

/// <param name="dimension"></param>

/// <returns></returns>

public static Matrix Zeros(int dimension)

{

Matrix result = new Matrix();

result.rows = dimension;

result.columns = dimension;

result.value = new double[dimension, dimension];

for (int i = 0; i < dimension; i++)

for (int j = 0; j < dimension; j++)

result.value[i, j] = 0;

return result;

}

/// <summary>

/// 生成一个指定行列数的零矩阵

/// </summary>

/// <param name="row"></param>

/// <param name="column"></param>

/// <returns></returns>

public static Matrix Zeros(int row, int column)

{

Matrix result = new Matrix();

result.rows = row;

result.columns = column;

result.value = new double[row, column];

for (int i = 0; i < row; i++)

for (int j = 0; j < column; j++)

result.value[i, j] = 0;

return result;

}

/// <summary>

/// 生成一个指定阶数的伪随机数矩阵

/// </summary>

/// <param name="dimension"></param>

/// <returns></returns>

public static Matrix Random(int dimension)

{

Matrix result = new Matrix();

result.rows = dimension;

result.columns = dimension;

result.value = new double[dimension, dimension];

for (int i = 0; i < dimension; i++)

for (int j = 0; j < dimension; j++)

result.value[i, j] = Matrix.GetRandomNum();

return result;

}

/// <summary>

/// 生成一个指定行列数的伪随机数矩阵

/// </summary>

/// <param name="row"></param>

/// <param name="column"></param>

/// <returns></returns>

public static Matrix Random(int row, int column)

{

Matrix result = new Matrix();

result.rows = row;

result.columns = column;

result.value = new double[row, column];

for (int i = 0; i < row; i++)

for (int j = 0; j < column; j++)

result.value[i, j] = Matrix.GetRandomNum();

return result;

}

/// <summary>

/// 生成一个指定阶数和随机范围的伪随机数矩阵

/// </summary>

/// <param name="dimension"></param>

/// <param name="min">伪随机数下限</param>

/// <param name="max">伪随机数上限</param>

/// <returns></returns>

public static Matrix Random(int dimension, double min, double max)

{

if (min > max)

{

min = max + min;

max = min - max;

min = min - max;

}

Matrix result = new Matrix();

result.rows = dimension;

result.columns = dimension;

result.value = new double[dimension, dimension];

for (int i = 0; i < dimension; i++)

for (int j = 0; j < dimension; j++)

result.value[i, j] = Matrix.GetRandomNum() * (max - min) + min;

return result;

}

/// <summary>

/// 生成一个指定行列数和伪随机数范围的矩阵

/// </summary>

/// <param name="row"></param>

/// <param name="column"></param>

/// <param name="min">伪随机数下限</param>

/// <param name="max">伪随机数上限</param>

/// <returns></returns>

public static Matrix Random(int row, int column, double min, double max)

{

if (min > max)

{

min = max + min;

max = min - max;

min = min - max;

}

Matrix result = new Matrix();

result.rows = row;

result.columns = column;

result.value = new double[row, column];

for (int i = 0; i < row; i++)

for (int j = 0; j < column; j++)

{

result.value[i, j] = Matrix.GetRandomNum() * (max - min) + min;

}

return result;

}

/// <summary>

/// 生成一个由指定数组元素填充对角线的对角矩阵

/// </summary>

/// <param name="diag"></param>

/// <returns></returns>

public static Matrix Diagonal(double[] diag)

{

Matrix result = Matrix.Zeros(diag.GetLength(0));

for (int i = 0; i < result.rows; i++)

{

result.value[i, i] = diag[i];

}

return result;

}

/// <summary>

/// 生成一个指定对角元素和行列数的对角矩阵

/// </summary>

/// <param name="diag"></param>

/// <param name="row"></param>

/// <param name="column"></param>

/// <returns></returns>

public static Matrix Diagonal(double[] diag,int row, int column)

{

Matrix result = Matrix.Zeros(row, column);

int num = row > column ? column : row;

num = num > diag.GetLength(0) ? diag.GetLength(0) : num;

for (int i = 0; i < num; i++)

{

result.value[i, i] = diag[i];

}

return result;

}

/// <summary>

/// 生成一个指定对角元素和偏移量的对角矩阵

/// </summary>

/// <param name="diag"></param>

/// <param name="move">偏移量为正表示向右移,为负表示向左移</param>

/// <returns></returns>

public static Matrix Diagonal(double[] diag, int move)

{

Matrix result = Matrix.Zeros(diag.GetLength(0));

for (int i = 0; i < result.rows; i++)

{

int col=(i + move) % result.columns;

col = col >= 0 ? col : (col + result.columns);

result.value[i, col] = diag[i];

}

return result;

}

/// <summary>

/// 生成一个取参数矩阵上三角部分元素矩阵

/// </summary>

/// <param name="x"></param>

/// <param name="fill">默认其他元素填充为0,可设置</param>

/// <returns></returns>

public static Matrix TriUp(Matrix x, double fill = 0)//上三角

{

Matrix result = new Matrix(x);

for (int i = 0; i < result.rows; i++)

{

for (int j = 0; j < result.columns; j++)

{

if (i > j)

{

result.value[i, j] = fill;

}

}

}

return result;

}

/// <summary>

/// 生成一个取参数矩阵下三角部分元素矩阵

/// </summary>

/// <param name="x"></param>

/// <param name="fill">默认其他元素填充为0,可设置</param>

/// <returns></returns>

public static Matrix TriLow(Matrix x, double fill = 0)//下三角

{

Matrix result = new Matrix(x);

for (int i = 0; i < result.rows; i++)

{

for (int j = 0; j < result.columns; j++)

{

if (i < j)

{

result.value[i, j] = fill;

}

}

}

return result;

}

/// <summary>

/// 生成一个对称矩阵

/// </summary>

/// <param name="x"></param>

/// <param name="reverse">默认为false,为true表示保留下三角部分元素</param>

/// <returns></returns>

public static Matrix Symmetry(Matrix x, bool reverse=false)//对称,inverse参数表示将上三角(下三角)对称到下方(上方)

{

Matrix result = new Matrix(x);

int num = x.rows < x.columns ? x.rows : x.columns;

result = Matrix.SubMatrix(x, num, num);

if (!reverse)

{

for (int i = 0; i < x.rows; i++)

{

for (int j = 0; j < x.columns; j++)

{

if (i > j)

{

result.value[i, j] = result.value[j, i];

}

}

}

}

else

{

for (int i = 0; i < x.rows; i++)

{

for (int j = 0; j < x.columns; j++)

{

if (i < j)

{

result.value[i, j] = result.value[j, i];

}

}

}

}

return result;

}

/// <summary>

/// 生成一个反对称矩阵

/// </summary>

/// <param name="x"></param>

/// <param name="reverse">默认为false,为true表示保留下三角部分元素</param>

/// <returns></returns>

public static Matrix Antisymmetry(Matrix x, bool reverse = false)//反对称,inverse参数表示将上三角(下三角)反对称到下方(上方)

{

Matrix result = new Matrix(x);

int num = x.rows < x.columns ? x.rows : x.columns;

result = Matrix.SubMatrix(x, num, num);

if (!reverse)

{

for (int i = 0; i < x.rows; i++)

{

for (int j = 0; j < x.columns; j++)

{

if (i > j)

{

result.value[i, j] = -result.value[j, i];

}

}

result.value[i, i] = 0;

}

}

else

{

for (int i = 0; i < x.rows; i++)

{

for (int j = 0; j < x.columns; j++)

{

if (i < j)

{

result.value[i, j] = -result.value[j, i];

}

}

result.value[i, i] = 0;

}

}

return result;

}

/// <summary>

/// 生成一个平面坐标旋转变换矩阵(结构力学中由局部坐标系到整体坐标系转换)

/// </summary>

/// <param name="angle">顺时针旋转角度值(不是弧度)</param>

/// <returns></returns>

public static Matrix TransMatrix(double angle)//转换矩阵,结构力学局部坐标系到整体坐标系转换

{

double cosVal = Math.Cos(angle * Math.PI / 180);

double sinVal = Math.Sin(angle * Math.PI / 180);

return new Matrix(new double[6, 6]{{cosVal,sinVal,0,0,0,0},{-sinVal,cosVal,0,0,0,0},{0,0,1,0,0,0},

{0,0,0,cosVal,sinVal,0},{0,0,0,-sinVal,cosVal,0},{0,0,0,0,0,1}});

}

/// <summary>

/// 生成一个平面刚架计算的刚度矩阵(结构力学)

/// </summary>

/// <param name="EI">弯曲刚度</param>

/// <param name="EA">拉压刚度</param>

/// <param name="L">刚架长度</param>

/// <returns></returns>

public static Matrix StiffnessMatrix(double EI, double EA,double L)//刚度矩阵//结构力学

{

return new Matrix(new double[6, 6]{{EA/L,0,0,-EA/L,0,0},{0,12*EI/(L*L*L),6*EI/(L*L),0,-12*EI/(L*L*L),6*EI/(L*L)},

{0,6*EI/(L*L),4*EI/L,0,-6*EI/(L*L),2*EI/L},{-EA/L,0,0,EA/L,0,0},

{0,-12*EI/(L*L*L),-6*EI/(L*L),0,12*EI/(L*L*L),-6*EI/(L*L)},

{0,6*EI/(L*L),2*EI/L,0,-6*EI/(L*L),4*EI/L}});

}

/// <summary>

/// 生成一个指定范围行向量(默认步距为1)

/// </summary>

/// <param name="begin">起始值</param>

/// <param name="end">终止值</param>

/// <returns></returns>

public static Matrix RangeVector(double begin, double end)

{

Matrix result = Matrix.Zeros(1, (int)Math.Floor(Math.Abs(end - begin))+1);

double incre = end >= begin ? 1 : -1;

int i = 0;

for (double x = begin; i<result.columns; x += incre)

{

result.value[0, i++] = x;

}

return result;

}

/// <summary>

/// 生成一个指定范围和步距的行向量

/// </summary>

/// <param name="begin">起始值</param>

/// <param name="incre">步距</param>

/// <param name="end">终止值</param>

/// <returns></returns>

public static Matrix RangeVector(double begin, double incre, double end)

{

if (end >= begin && incre < 0)

throw new Exception("IncreaseVector函数使用时试图由小数递减到大数");

if (end <= begin && incre > 0)

throw new Exception("IncreaseVector函数使用时试图由大数递增到小数");

if (incre == 0)

throw new Exception("IncreaseVector函数使用时递增量为0错误");

Matrix result = Matrix.Ones(1, (int)Math.Floor((end - begin) / incre) + 1);

decimal val = (decimal)begin;

for (int i = 0; i < result.columns; i++)

{

result.value[0, i] = result.value[0, i] * (double)val;

val += (decimal)incre;

}

return result;

}

/// <summary>

/// 生成一个线性空间行向量

/// </summary>

/// <param name="begin">起始值</param>

/// <param name="end">终止值</param>

/// <param name="num">元素个数</param>

/// <returns></returns>

public static Matrix LinspaceVector(double begin, double end, int num)

{

Matrix result = new Matrix(1, num + 1);

for (int i = 0; i <= num; i++)

{

result.value[0, i] = i * (end - begin) / num + begin;

}

return result;

}

/// <summary>

/// 生成一个对数空间行向量

/// </summary>

/// <param name="begin">起始值</param>

/// <param name="end">终止值</param>

/// <param name="num">元素个数</param>

/// <returns></returns>

public static Matrix LogspaceVector(double begin, double end, int num)

{

Matrix result = new Matrix(1, num + 1);

for (int i = 0; i <= num; i++)

{

result.value[0, i] = Math.Pow(10, i * (end - begin) / num + begin);

}

return result;

}

#endregion

#region 矩阵的控制台显示

/// <summary>

/// 在控制台显示矩阵

/// </summary>

public void DisplayInConsole()

{

for(int i=0;i<rows;i++)

for (int j = 0; j < columns; j++)

{

if (j == 0)

Console.Write("[");

Console.Write("{0}\t", value[i, j]);

if (j == columns - 1)

Console.WriteLine("]");

}

}

/// <summary>

/// 在控制台显示矩阵(元素取四位小数)

/// </summary>

public void DisplayLimited()

{

for (int i = 0; i < rows; i++)

for (int j = 0; j < columns; j++)

{

if (j == 0)

Console.Write("[");

Console.Write("{0:f4}\t", value[i, j]);

if (j == columns - 1)

Console.WriteLine("]");

}

}

#endregion

#region 矩阵的运算处理

/// <summary>

/// 返回参数矩阵的对角元素

/// </summary>

/// <param name="x"></param>

/// <returns></returns>

public static double[] Diagonal(Matrix x)//其实是生成对角矩阵方法的反方法,取矩阵的对角元素

{

int num = x.rows > x.columns ? x.columns : x.rows;

double[] result = new double[num];

for (int i = 0; i < num; i++)

{

result[i] = x.value[i, i];

}

return result;

}

/// <summary>

/// 返回参数矩阵的逆矩阵

/// </summary>

/// <param name="x"></param>

/// <returns></returns>

public static Matrix Reverse(Matrix x)

{

if (x.rows != x.columns)

throw new Exception("Reverse函数使用时,求逆的矩阵必须为方阵,行列数必须相同");

if (x.rows == 1)

if (x.value[0, 0] != 0)

{

return new Matrix(1 / x.value[0, 0]);

}

else

throw new Exception("Reverse函数使用时,对数字0求逆错误");

double detVal=Matrix.Det(x);

if (detVal == 0)

throw new Exception("Reverse函数使用时,方阵的行列式为0不可逆");

detVal = 1 / detVal;

return Matrix.CompanionMatrix(x) * detVal;

}

/// <summary>

/// 返回参数矩阵i行j列的余子式

/// </summary>

/// <param name="x">注意必须为方阵</param>

/// <param name="i">i从0起到行数减1</param>

/// <param name="j">j从0起到列数减1</param>

/// <returns></returns>

public static double Cofactor(Matrix x, int i, int j)//求逆的辅助函数，求余子式，注意判断行列数相同

{

double[,] array = new double[x.rows - 1, x.columns - 1];

for (int m = 0; m < x.rows - 1; m++)

{

for (int n = 0; n < x.columns - 1; n++)

{

int s = m, t = n;

if (s >= i)

s++;

if (t >= j)

t++;

array[m, n] = x.value[s, t];

}

}

Matrix cofactor = new Matrix(array);

return Matrix.Det(cofactor);

}

/// <summary>

/// 返回参数矩阵的代数余子式

/// </summary>

/// <param name="x">注意必须为方阵</param>

/// <param name="i">i从0起到行数减1</param>

/// <param name="j">j从0起到列数减1</param>

/// <returns></returns>

public static double AlgeCofactor(Matrix x, int i, int j)//代数余子式

{

return Math.Pow(-1, i + j) * Matrix.Cofactor(x, i, j);

}

/// <summary>

/// 返回参数矩阵的伴随矩阵

/// </summary>

/// <param name="x">注意必须为方阵</param>

/// <returns></returns>

public static Matrix CompanionMatrix(Matrix x)//伴随矩阵

{

double[,] array = new double[x.rows, x.columns];

for (int i = 0; i < x.rows; i++)

{

for (int j = 0; j < x.columns; j++)

{

array[i, j] = Matrix.AlgeCofactor(x, i, j);

}

}

return Matrix.Transfer(new Matrix(array));

}

/// <summary>

/// 返回参数矩阵的转置矩阵

/// </summary>

/// <param name="x"></param>

/// <returns></returns>

public static Matrix Transfer(Matrix x)

{

Matrix result = new Matrix();

result.rows = x.columns;

result.columns = x.rows;

result.value = new double[result.rows, result.columns];

for (int i = 0; i < result.rows; i++)

for (int j = 0; j < result.columns; j++)

result.value[i, j] = x.value[j, i];

return result;

}

/// <summary>

/// 返回参数矩阵的行列式值

/// </summary>

/// <param name="x">注意必须为方阵</param>

/// <returns></returns>

public static double Det(Matrix x)

{

if(x.rows!=x.columns)