线性代数

Posted 2021-02-20 bzdyl

线性代数

矩阵消元

总的来说就是第n行消去从n+1行开始所有行的第n个元素

行列式

定义：

[对于n阶矩阵A=left[egin{matrix}a_{1,1} & …… &a_{1,n} \ & …… & \a_{n,1} & …… & a_{n,n}end{matrix} ight] ]

记(A)的行列式为(|A|)或(det(A))

[egin{align} & det(A)=sum_{j_1j_2……j_n}(-1)^{r(j_1,……j_n)}cdot a_{1,j_1}a_{2,j_2}……a_{n,j_n}& 其中{j_1,……,j_n}是1-n的一个排列，r(j_1,……,j_n)表示这个排列的逆序对数量 end{align} ]

上三角矩阵

形如

[A= left[ egin{matrix} a_{1,1} & a_{1,2} & a_{1,3} & a_{2,2} & a_{2,3} & 0 & a_{3,3} end{matrix} ight] ]

这样主对角线下方全为0的矩阵

对于一个上三角矩阵(A)，它的行列式(det(A))就等于主对角线乘积

行列式计算

性质：对于一个矩阵，在对它进行矩阵消元之后，行列式的值不改变

根据这一性质，我们可以先对矩阵进行消元，将它化为上三角矩阵的形式，再求解

代码：

int a[N][N],ans=1;

void Matrix_elimination()
{
	for(int i=2;i<=n;i++)
		for(int j=1;j<i;j++)
		{
			double tem=a[i][j]/a[j][j];
			for(int k=1;k<=n;k++) a[i][k]-=tem*a[j][j];
		}
	return;
}

for(int i=1;i<=n;i++) ans*=a[i][i];

线性方程组

形如

[egin{cases} a_1x_1+b_1x_2+cdots+m_1x_n=k_1a_2x_1+b_2x_2+cdots+m_2x_n=k_2\cdotsa_nx_n+b_nx_2+cdots+m_nx_n=k_n end{cases} ]

对于这样的方程，我们可以将他们变成一个形如

[left[ egin{matrix} a_1,cdots,m_1,k_1\cdotsa_n,cdots,m_n,k_n end{matrix} ight] ]

的增广矩阵，然后运用上面所说的矩阵消元，再依次回带就可以得出方程的解

高斯消元法：

1.把方程组转换成增广矩阵。

2.利用初等行变换来把增广矩阵转换成行阶梯阵。
枚举k从0到equ – 1，当前处理的列为col(初始为0) ，每次找第k行以下(包括第k行)，col列中元素绝对值最大的列与第k行交换。如果col列中的元素全为0，那么则处理col + 1列，k不变。

3.转换为行阶梯阵，判断解的情况。

① 无解
当方程中出现(0, 0, …, 0, a)的形式，且a != 0时，说明是无解的。

② 唯一解
条件是k = equ，即行阶梯阵形成了严格的上三角阵。利用回代逐一求出解集。

③ 无穷解。
条件是k < equ，即不能形成严格的上三角形，自由变元的个数即为equ – k，但有些题目要求判断哪些变元是不缺定的。
这里单独介绍下这种解法：
首先，自由变元有var - k个，即不确定的变元至少有var - k个。我们先把所有的变元视为不确定的。在每个方程中判断不确定变元的个数，如果大于1个，则该方程无法求解。如果只有1个变元，那么该变元即可求出，即为确定变元。

代码：

#include <iostream>
#include <cmath>
#include <algorithm>

using namespace std;

const double eps = 1e-6;

bool is_zero( const double num )//用于判断有无逆矩阵
{
    return fabs(num) < eps;
}

void create( double ** & matrix, const int n )
{
    matrix = new double* [n];
    for ( int i = 0; i < n; ++i )
        matrix[i] = new double[n];
}

void input ( double ** matrix, const int n )
{
    for ( int i = 0; i < n; ++i )
    {
        for ( int j  = 0; j < n; ++ j )
            cin >> matrix[i][j];
    }
}

bool inverse ( double ** matrix1, double ** matrix2, const int n )
{
    int i, j;
    for ( i = 0; i < n; ++ i )//初始化一个单位矩阵
    {
        for ( j = 0; j < n; ++ j )
        {
            if ( i == j )
                matrix2[i][j] = 1;
            else
                matrix2[i][j] = 0;
        }
    }
    for ( i = 0; i < n; ++i )
    {
        int rowmaxpos = i;
        for ( j = i + 1; j < n; ++j )
        {
            if ( matrix1[i][j] > matrix1[i][rowmaxpos] )
                rowmaxpos = j;
        }
        for ( j = i; j < n; ++ j )//按从大到小的顺序排列矩阵
        {
            swap( matrix1[j][rowmaxpos], matrix1[j][i]);
            swap( matrix2[j][rowmaxpos], matrix2[j][i]);
        }
        if ( !is_zero(matrix1[i][i]) )
        {
            int divisor = matrix1[i][i];
            for ( j = i; j < n; ++ j )//归一化矩阵
            {
                matrix1[i][j] /= divisor;
                matrix2[i][j] /= divisor;
            }
            for ( j = i + 1; j < n; ++ j )//高斯消元法处理行列式
            {
                int multiple = matrix1[j][i];
                for ( int k = i; k < n; ++ k )
                {
                    matrix1[i][j] -= matrix1[i][k] * multiple;
                    matrix2[i][j] -= matrix2[i][k] * multiple;
                }
            }
        }
        else
            return false;
    }
    return true;
}

void output( double ** matrix, const int n )
{
    for ( int i = 0; i < n; ++i )
    {
        for ( int j = 0; j < n; ++ j )
            cout << matrix[i][j] << ‘ ‘;
        cout<<endl;
    }
}

void destroy( double ** matrix, const int n )
{
    for ( int i = 0; i < n; ++ i )
        delete [] matrix[i];
    delete [] matrix;
}

int main()
{
    int n;
    double ** matrix1;
    double ** matrix2;
    while ( cin >> n )
    {
        create( matrix1, n );
        create( matrix2, n );
        input( matrix1, n);
        if ( inverse(matrix1, matrix2, n) )
            output( matrix2, n );
        else
            cout << "No inverse matrix" << endl;
        destroy( matrix1, n );
        destroy( matrix2, n );
    }
    return 0;
}