数值方法No.4-求解联立线性代数方程

主要内容

• 矩阵定义
• 解的唯一性
• 顺序Gauss消去法
• 全选主元Gauss消去法
• LU分解
• Cholesky分解
• Gauss-Jordan消去法&求矩阵的逆
• 三对角线方程组
• Gauss-Seidel迭代&松弛迭代
• 病态系统

一、引言

• n个方程和n个未知数的方程组

• 可以转化为矩阵的格式

• 定义增广矩阵

• 解线性方程的方法：
  • 直接求解：Gauss消元，LU分解，Cholesky分解，Thomas法，Gauss-Jordan法。
  • 迭代：Gauss-Seidel。

二、定义矩阵

• 
• 对阵矩阵：
• 斜对称矩阵：
• 对角阵：主对角线存在非零元素的方正，其他都为0。
• 上三角U：下半部分全0；下三角L：上半部分全0。
• 稠密矩阵：没有任何0元素；稀疏矩阵：具有个别0元素。
• 迹tr：方阵对角线的元素之和。
• 非奇异矩阵：方阵的行列式不为0，从而有逆矩阵。
• 当行列式很小，矩阵是病态的。
• 秩：最大的非奇异的子方阵的阶数。

三、解的唯一性（以下默认系数A矩阵为n*n的方阵）

• 当且仅当：增广矩阵和系数矩阵 的秩相同，方程组有解；
• 当且仅当：增广矩阵和系数矩阵 的秩相同，且=未知数的个数，方程组有唯一解；
• 当且仅当：增广矩阵和系数矩阵 的秩相同，且<未知数的个数，方程组有无数解；
• 齐次方程组：——总有解为全0向量，这种解称为平凡解。如果系数矩阵奇异，那么齐次方程组存在无数解。

四、顺序Gauss消去法

1. 原理和算法

• 将增广矩阵化简为上三角矩阵，从而向量在此过程中被修改，最后解出向量。

消去。消去和回代步骤的公式如下：

• 消去

• 回代：从计算到

• 缺陷：
  • 如果对角线上的主元素=0的时候，那么消去公式的计算就会出现问题。
  • 当这个元素和其他元素比起来很小的时候，那么结果可能不准确。
  • 解决方法：交换系数矩阵的行和列，一边选择最大的主元素，称为全主元法；不时执行交换行或者列，而不是全部都交换，称为部分主元法。如果对角线选主元找不到非零元素，那么矩阵一定奇异。
• 算法复杂度：
  • 消去和回代需要次乘法和除法运算，以及次加减运算。

2. 例子

• 例子

首先，系数矩阵和增广矩阵分别为
以及

消去：根据公式，n=3, kij依次取数：
k=0, i=1,2, j=0,1,2,3;
k=1, i=2, j=1,2,3；
①k=0, i=1的时候，j=0~3

k=0, i=2的时候，j=0~3

因此，k=0的循环之后，结果为

②k=1, i=2的时候，j=1,2,3

因此，k=1的循环之后，结果为

接着进行回代：
。

3. 代码

#include<iostream>
#include<math.h>
using namespace std;
class Gauss
{
private:
  int i,j,k,n;
  double eps, ratio, sum, *x, **a; //定义x是用来储存求解的结果向量，a是用来储存增广矩阵。
public:
  void gauss_input();
  void gauss_elimination();
  void gauss_output();
  ~Gauss()
  {
    delete [] x; // 释放x向量空间
    // 释放a矩阵空间，需要先释放每个“元素”（每个元素是向量），再最后释放a。
    for(i=0; i<n; i++)
      {
        delete [] a[i];
      }
       delete [] a;
  }
};
void Gauss::gauss_input()
{
  cout << "Enter the number of functions : " ;
  cin >> n;
  // 储存x向量
  x = new double [n]; 
  // 出错点：这里*不可少，储存a矩阵；注意定义矩阵的技巧。
  a = new double*[n];
  for(i=0; i<n; i++){
    a[i] = new double [n+1];
  }
  // 输入a(0:n-1, 0:n-1)的元素即为系数矩阵。
  for(i=0; i<n; i++)
    for(j=0; j<n; j++){
      cout << "Enter a[" << i << "][" << j << "]=";
      cin >> a[i][j];
    }
  // 输入a(0:n-1, n)的元素即为b
  for(i=0; i<n; i++){
    cout << "Enter b[" << i << "]=" ;
    cin >> a[i][n];
  }
  // 输入对角线最小的元素。
  cout << "Min of main element : ";
  cin >> eps;
}
void Gauss::gauss_elimination()
{
  // 消去
  for(k=0; k<n-1; k++){
    int i_start=k+1;
    //如果对角线主元素太小，那么报错。
    if(fabs(a[k][k]) < eps){
      cout << "Main element is too small. Failure..." << endl;
      exit(0);
    }
    for(i=i_start; i<n; i++){
      ratio = a[i][k]/a[k][k];
      for(j=(k+1); j<n+1; j++){
    a[i][j] -= ratio * a[k][j];
      }
    }
  }
  // 回代
  x[n-1] = a[n-1][n] / a[n-1][n-1];
  for(i=n-2; i>=0; i--){
    sum = 0.0;
    for(j=i+1; j<n; j++){
      sum += a[i][j] * x[j];
    }
    x[i] = (a[i][n] - sum) / a[i][i];
  }
}
void Gauss::gauss_output()
{
  cout << "Result is " << endl;
  for(i=0; i<n; i++){
    cout << "x[" << i << "]=" << x[i] << endl;
  }
}
int main()
{
  Gauss g1;
  g1.gauss_input();
  g1.gauss_elimination();
  g1.gauss_output();
  return 0;
}

ubuntu中运行

danielshen@DESKTOP-VOSNAQJ:~/nr/chp2$ g++ c5.cpp -o c5 -std=c++11
danielshen@DESKTOP-VOSNAQJ:~/nr/chp2$ ./c5
Enter the number of functions : 3
Enter a[0][0]=0.8
Enter a[0][1]=0.02
Enter a[0][2]=0.06
Enter a[1][0]=0.1
Enter a[1][1]=0.83
Enter a[1][2]=0.12
Enter a[2][0]=0.1
Enter a[2][1]=0.15
Enter a[2][2]=0.82
Enter b[0]=50
Enter b[1]=30
Enter b[2]=20
Min of main element : 0.001
Result is
x[0]=60.9226
x[1]=27.0682
x[2]=12.0091

4. 其他补充

• 高斯消元法不适合奇异矩阵。
• 行阶梯形矩阵（Row echelon form）：每行的第一个非零元素（称为leading coefficient）为1；第i行的leading coefficient 位于第i+1行的leading coefficient的左上；全0的行位于具有非零行的下方。
• 主要是进行行操作，从而将矩阵转化为上三角矩阵，变成行阶梯形矩阵；操作包括了“交换行”、“将行乘以非零常数”、“将一行加到另一行。”。

• 算法：

  • 部分旋转：通过交换行，从而将具有最大绝对值的输入移动到轴的位置，保持计算的稳定性；

• 代码

• 
  

五、全选主元Gauss消去法

https://www.geeksforgeeks.org/program-for-gauss-jordan-elimination-method/

1. 介绍

• 全选主元Gauss消去法更加robust。
• 新增了两个数据成员pivrow和pivcol，用于储藏在执行过程中的主元元素的行和列，储存需要的信息。因为在“行列”交换的时候，x也需要进行交换行的位置，因此要用到pivcol中的信息。

2. 例子

例如，方程组，直接求解的话，，。原因是太小了，导致在求解的时候变得巨大。

那么这里考虑将上下进行交换位置，就能得到正确的解。
，

因此，在进行高斯消去的时候，避免采用绝对值很小的元素，保持。

3. 算法

• 第一步
  • 选主元：在系数矩阵A中选出绝对值最大的元素作为主元素，确定坐标，使得
  • 交换行列式：
    • 当，交换增广矩阵的第0行和第行元素；
    • 当，交换增广矩阵的第0列列元素；
    • 注意调换和，并记录。
  • 消元
    • 
    • 
• 第k步：重复进行，假设已经完成了第1步~第k-1步，
• 消去：对于，从到进行以下运算。
  • 全选主元：
  • 通过行列交换，把绝对值最大的元素交换到主元位置上。
  • 系数矩阵归一化：

待完善

六、带有正向和反向代入的LU分解

1. 介绍

• 如果，系数矩阵元素不变，但是右边的b向量变了，那么对系数矩阵的操作不必重复进行，可以节省大量的计算工作。
• 因此，采用LU分解，L是下三角矩阵，U是上三角矩阵。从而系数矩阵能够表示为两个矩阵的乘积

• 注意，下三角矩阵L对角线元素是1，从而使用前向和后向扫描得到了向量，定义新的向量，那么有

利用进行前向扫描，得到，利用进行后向扫描得到。

2. 算法

• 步骤1：设
• 步骤2：对于，执行
  • 
  • 
• 系数矩阵A中的每个元素只使用1次，因此可以储存在A的对应位置上，从而分解为

• 前向代入

• 后向带入

元素不需要储存。

1. 例子