2021-02-26-数值分析-Day01-顺序Gauss消元

1.1 - 1.2

数值分析的算法要求

• 结构简单，易于计算机实现
• 理论上保证方法的收敛性和数值稳定性
• 计算效率高：速度快，内存开销小
• 经过数值实验验证

误差的来源和分类

• • 模型误差：将实际问题抽象为数学模型过程中产生的误差，不可避免
  • 观测误差：观测实验得到的参数带来的误差
  • 截断误差：近似方法产生的误差
  • 舍入误差：计算机有限位计算产生的误差

  （后两者为本课程主要研究对象）

• 绝对误差和相对误差

  设 \(x\) 是精确值 \(x^{*}\) 的一个近似值，则 \(e=x^{*}-x\) 为近似值 \(x\) 的绝对误差，简称误差。

  若 \(\varepsilon\) 满足 \(|e|\le\varepsilon\) ，则称 \(\varepsilon\) 为 \(x\) 的绝对误差限，简称误差限，有量纲。则满足 \(x-\varepsilon \le x^{*} \le x+\varepsilon\) ，简写为 \(x^{*} = x\pm\varepsilon\) 。

  \(e_{r} = \displaystyle\frac{e}{x^{*}} = \frac{x^{*}-x}{x^{*}}\) 称相对误差。

  \(\varepsilon_{r} = \displaystyle\frac{\varepsilon}{|x|}\) 称相对误差限， \(|e_{r}| \le \varepsilon_{r}\) ，无量纲。

  （其中误差与相对误差用于理论分析，误差限与相对误差限用于实际应用。）

1.3 有效数字

定义1：设数 \(x\) 是数 \(x^{*}\) 的近似值，如果 \(x\) 的绝对误差限是它的某一数位的半个单位，并且从 \(x\) 左起第一个非零数到该数位共有 \(n\) 位，则称这 \(n\) 个数字为 \(x\) 的有效数字，也称用 \(x\) 近似 \(x*\) 时具有 \(n\) 位有效数字。

有效数字与绝对误差限的关系： \(x\) 作为 \(x^{*}\) 的近似值，具有 \(n\) 位（ \(n \le k\) ）有效数字当且仅当：

\[ \left|x^{*}-x\right| \le \frac{1}{2}\times10 ^{m-n} \]

近似值的有效数字越多，其绝对误差越小。

（PS：精确值的有效数字可认为有无限多位。）

有效数字与相对误差限的关系：若 \(x\) 有 \(n\) 位有效数字，则其相对误差限为： \[ \varepsilon_{r} \le \frac{1}{2 a_{1}}\times10^{-n+1} \] 反之，若 \(x\) 的相对误差限： \[ \varepsilon_{r} \le \frac{1}{2(a_{1}+1)}\times10^{-n+1} \] 则 \(x\) 至少有 \(n\) 位有效数字。

1.4 - 1.6 数值计算中的若干原则

为减少舍入误差

• 避免两个相近的数相减

  利用对数计算法则将减法转换为除法、利用三角函数和差化积公式转换为乘除、倒用分母有理化将减法转换为分母上的乘法等方法。

• 防止大数“吃掉”小数

  指参与计算的数数量级差很大时，则加减运算中，绝对值小的数往往被绝对值大的数“吃掉”（如被计算机存储位数限制舍去）。

  改变计算顺序，在求和差过程中由小到大进行计算。

• 绝对值太小的数不宜作除数

• 注意简化计算程序，减少计算次数

• 选用数值稳定性好的算法

  计算舍入误差积累是可控制的。

2.1 - 2.4 Gauss消去法

求解线性方程组 \[ \left\{ \begin{aligned} a_{11}x_{1} + a_{12}x_{2} + ... + a_{1n}x_{n} = b_{1}\\ a_{21}x_{1} + a_{22}x_{2} + ... + a_{2n}x_{n} = b_{2}\\ \dots\\ a_{n1}x_{1} + a_{n2}x_{2} + ... + a_{nn}x_{n} = b_{n}\\ \end{aligned} \right. (1) \] 矩阵形式 \(Ax = b\) \[ A= \begin{pmatrix} a_{11}& a_{12}& \cdots & a_{1n} \\ a_{21}& a_{22}& \cdots & a_{2n} \\ \vdots & \vdots & \ddots & \vdots \\ a_{n1}& a_{n2}& \cdots & a_{nn} \end{pmatrix}, x= \begin{pmatrix} x_{1} \\ x_{2} \\ \cdots \\ x_{n} \\ \end{pmatrix}, b= \begin{pmatrix} b_{1} \\ b_{2} \\ \cdots \\ b_{n} \\ \end{pmatrix} \\ 矩阵A奇异，det(A)\neq0 \]

Cramer法则：直接方法，计算量过大，不适合实际应用

Gauss消去法：直接方法，逐次消元，转换为上三角方程组。

顺序Gauss消去法

顺序Gauss消去法求解线性方程组 \((1)\)：

记 \(A^{(1)}=A, b^{(1)}=b, a_{ij}^{(1)}=a_{ij}, b_{i}*{(1)}=b_{i}\)

则，线性方程组 \((1)\) 的增广矩阵为： \[ (A^{(1)},b^{(1)})= \begin{pmatrix} a_{11}^{(1)}& a_{12}^{(1)}& a_{13}^{(1)}& \dots & a_{1n}^{(1)}& b_{1}^{(1)} \\ a_{21}^{(1)}& a_{22}^{(1)}& a_{23}^{(1)}& \dots & a_{2n}^{(1)}& b_{2}^{(1)} \\ a_{31}^{(1)}& a_{32}^{(1)}& a_{33}^{(1)}& \dots & a_{3n}^{(1)}& b_{3}^{(1)} \\ \dots& \dots& \dots& \dots& \dots& \dots \\ a_{n1}^{(1)}& a_{n2}^{(1)}& a_{n3}^{(1)}& \dots & a_{nn}^{(1)}& b_{n}^{(1)} \\ \end{pmatrix}, \]

• 第 \(1\) 步：设 \(a_{11}^{(1)}\neq0\) ，依次用 \(-l_{i1}=-\displaystyle\frac{a_{i1}^{(1)}}{a_{11}^{(1)}}, (i=2,3,\dots,n)\) 乘以矩阵的第 \(1\) 行加到第 \(i\) 行，得到矩阵： \[ \begin{array}{l} (A^{(2)},b^{(2)})= \begin{pmatrix} a_{11}^{(1)}& a_{12}^{(1)}& a_{13}^{(1)}& \dots & a_{1n}^{(1)}& b_{1}^{(1)} \\ 0& a_{22}^{(2)}& a_{23}^{(2)}& \dots & a_{2n}^{(2)}& b_{2}^{(2)} \\ 0& a_{32}^{(2)}& a_{33}^{(2)}& \dots & a_{3n}^{(2)}& b_{3}^{(2)} \\ \dots& \dots& \dots& \dots& \dots& \dots \\ 0& a_{n2}^{(2)}& a_{n3}^{(2)}& \dots & a_{nn}^{(2)}& b_{n}^{(2)} \\ \end{pmatrix},\\ 其中：\begin{aligned} a_{ij}^{(2)} = a_{ij}^{(1)} - l_{i1}a_{1j}^{(1)},\quad i,j=2,3,\dots,n \\ b_{i}^{(2)} = b_{i}^{(1)} - l_{i1}b_{1}^{(1)},\qquad i=2,3,\dots,n \end{aligned} \end{array} \]

• 第 \(2\) 步：设 \(a_{22}^{(2)}\neq0\) ，依次用 \(-l_{i2}=-\displaystyle\frac{a_{i2}^{(2)}}{a_{22}^{(2)}}, (i=3,4,\dots,n)\) 乘以矩阵的第 \(2\) 行加到第 \(i\) 行，得到矩阵： \[ \begin{array}{l} (A^{(2)},b^{(2)})= \begin{pmatrix} a_{11}^{(1)}& a_{12}^{(1)}& a_{13}^{(1)}& \dots & a_{1n}^{(1)}& b_{1}^{(1)} \\ 0& a_{22}^{(2)}& a_{23}^{(2)}& \dots & a_{2n}^{(2)}& b_{2}^{(2)} \\ 0& 0& a_{33}^{(3)}& \dots & a_{3n}^{(3)}& b_{3}^{(3)} \\ \dots& \dots& \dots& \dots& \dots& \dots \\ 0& 0& a_{n3}^{(3)}& \dots & a_{nn}^{(3)}& b_{n}^{(3)} \\ \end{pmatrix},\\ 其中：\begin{aligned} a_{ij}^{(3)} = a_{ij}^{(2)} - l_{i2}a_{2j}^{(2)},\quad i,j=3,4,\dots,n \\ b_{i}^{(3)} = b_{i}^{(2)} - l_{i2}b_{2}^{(2)},\qquad i=3,4,\dots,n \end{aligned} \end{array} \]

• 第 \(3\) 步，第 \(4\) 步， \(\dots\) ，第 \((n-1)\) 步后得到矩阵： \[ (A^{(n)},b^{(n)})= \begin{pmatrix} a_{11}^{(1)}& a_{12}^{(1)}& a_{13}^{(1)}& \dots & a_{1n}^{(1)}& b_{1}^{(1)} \\ 0& a_{22}^{(2)}& a_{23}^{(2)}& \dots & a_{2n}^{(2)}& b_{2}^{(2)} \\ 0& 0& a_{33}^{(3)}& \dots & a_{3n}^{(3)}& b_{3}^{(3)} \\ \dots& \dots& \dots& \dots& \dots& \dots \\ 0& 0& 0& \dots & a_{nn}^{(n)}& b_{n}^{(n)} \\ \end{pmatrix} \]

• 对应的方程组则变成： \[ \left\{ \begin{aligned} a_{11}^{(1)}x_{1} + a_{12}^{(1)}x_{2} + ... + a_{1n}^{(1)}x_{n} = b_{1}^{(1)}\\ a_{22}^{(2)}x_{2} + ... + a_{2n}^{(2)}x_{n} = b_{2}^{(2)}\\ \dots\\ a_{nn}^{(n)}x_{n} = b_{n}^{(n)}\\ \end{aligned} \right. \]

• 对此方程组进行回代，就可求出方程组的解

顺序Gauss消去法的乘除计算量：\(\displaystyle\sum_{k=1}^{n}(k^{2}-1)\ +\ \displaystyle\sum_{k=1}^{n}k\ =\ \displaystyle\frac{1}{3}(n^{3}+3n^{2}-n)\)

顺序Gauss消去法也简称为Gauss消去法，其中 \(a_{kk}^{(k)}(k=1,2,\dots,n)\) 也称为主元素，且主元素不为零（使用顺序Gauss消去法的充要条件）。

列主元Gauss消去法

顺序Gauss消去法的局限性：主元素过小时出现“大数吃小数”情况，舍入误差较大。常采用列主元Gauss消去法和全主元Gauss消去法解决。

线性方程组 \(Ax=b\)，记 \(A^{(1)}=A, b^{(1)}=b,\) 列主元Gauss消去法过程如下： + 在增广矩阵 \(B^{(1)}=(A^{(1)},b^{(1)})\) 的第一列元素中，取： \(\left|a_{k1}^{(1)}\right| = \max\limits_{1\le i\le n}{\left|a_{i1}^{(1)}\right|}\) 为主元素， \(r_{k}\leftrightarrow r_{1}\) 。 然后进行第 \(1\) 步消元得到增广矩阵 \(B^{2}=(A^{(2), b^{(2)}})\) 。 + 再在矩阵 \(B^{(2)}=(A^{(2)},b^{(2)})\) 的第二列元素中，取： \(\left|a_{k2}^{(2)}\right| = \max\limits_{2\le i\le n}{\left|a_{i2}^{(2)}\right|}\) 为主元素， \(r_{k}\leftrightarrow r_{2}\) 。 然后进行第二步消元得到增广矩阵 \(B^{3}=(A^{(3), b^{(3)}})\) 。 + 依此类推，经过 \((n-1)\) 步选主元和消元运算，得到增广矩阵 \(B^{n}=(A^{(n), b^{(n)}})\) 。则方程组 \(A^{(n)}x=b^{(n)}\) 是与原方程组等价的上三角方程组，回代求解。 + 易证：只要 \(|A|\neq 0\) 列主元Gauss消去法就可顺利进行。

全主元Gauss消去法

在每一步消元前，在所有元素中选取绝对值最大的元素作为主元素。但由于运算量大增，实际应用中不经常使用。

代码：顺序Gauss消去法

在每一步消元前，在所有元素中选取绝对值最大的元素作为主元素。但由于运算量大增，实际应用中不经常使用。

//
// Created by xa on 2021-02-26.
//
#include <iostream>
#include <vector>

using std::vector;

vector<vector<double>> matrix;
int size = 0;

vector<double> solve(vector<vector<double>> input);
void upperTriangular();
void swapRow(int index);
void gauss(int index);

template<typename T>
void swap(T* a, T* b);

void printMatrix();

int main() {
    vector<vector<double>> test {{8, 1, 7, 3}, {3, 7, 9, 8}, {9, 1, 5, 6}};
    vector<double> ans = solve(test);
    for (int i = 0; i < ans.size(); ++i)
        std::cout << ans[i] << " ";
}

vector<double> solve(vector<vector<double>> input) {
    matrix.assign(input.begin(), input.end());
    size = matrix.size();
    vector<double> solution(size);
    upperTriangular();

    if (matrix[size-1][size-1] == 0) std::cout << "No Solution.";

    solution[size-1] = matrix[size-1][size]/matrix[size-1][size-1];
    for (int i = size-2; i >= 0; --i) {
        double tmp = matrix[i][size];
        for (int j = i+1; j < size; ++j)
            tmp -= matrix[i][j] * solution[j];
        solution[i] = tmp/matrix[i][i];
    }

    return solution;
}

void upperTriangular() {
    for (int i = 0; i < size; ++i) {
        swapRow(i);
        gauss(i);
    }
    printMatrix();
}

void swapRow(int index) {
    int max = index;
    if (index < size - 1) {
        for (int i = index; i < size; ++i)
            if (matrix[i][index] > matrix[max][index]) max = i;
        swap(matrix[index], matrix[max]);
    }
    std::cout << "swap" << std::endl;
    printMatrix();
}

void gauss(int index) {
    for (int i = index + 1; i < size; ++i) {
        double l = matrix[i][index] / matrix[index][index];
        for (int j = 0; j < size + 1; ++j) {
            matrix[i][j] -= l * matrix[index][j];
            if (matrix[i][j] < 1e-16 && matrix[i][j] > -1e-10) matrix[i][j] = 0;
        }
    }
    std::cout << "gauss" << std::endl;
    printMatrix();
}

template<typename T>
void swap(T* a, T* b) {
    double tmp = a;
    a = b;
    b = tmp;
}

void printMatrix() {
    for (int i = 0; i < size; ++i) {
        for (int j = 0; j < size + 1; ++j)
            std::cout << matrix[i][j] << ' ';
        std::cout << std::endl;
    }
    std::cout << std::endl;
}