Gamer's Show

你知道人生最要紧的事是快乐不停

0%

数值计算方法

Posted on Edited on In

引言

误差

  1. 对于十字相乘法类的二次方程，如果两根相差过大，会发现求根公式算小根是错误的：或许结论推理也不能
  2. 如何运算精确度更高？

误差的来源

知道误差来源属于哪一类

  1. 模型误差：分析方法近似造成的误差eg.物理中的理想化模型
  2. 观测误差：测量过程中不可消除的误差，与设备、方式有问题
  3. 截断误差：数学计算方法中常见，主动的操作，可以保证计算最后的精度；数学中的“无穷”需要转化为计算机中的“有穷”；比如说无穷大化为一个极大数值，无穷小数转为多位的有限小数。表示为：
    $$

$$
* 泰勒展开，级数展开
4. 舍入误差：计算机中由于数据类型的长度有限，在实际位数超过该长度时，
在本课误差估计中，假定数学模型合理，测定值精确，主要讨论截断误差和舍入误差

绝对误差、相对误差与有效数字

绝对误差

$x^*$为准确值x的一个近似值，则称$x-x^*$为近似值$x^*$的绝对误差，简称误差，表示为$\epsilon$:
$$ \epsilon=x-x^*$$
$|\epsilon|$的大小显示出近似值的准确程度，同一量的不同近似中，$|\epsilon|$越小，近似值准确度越高

绝对误差限

实践中无法确定x，但是可以确定误差绝对值的上限，从而确定准确值的范围，定义满足如下不等式的正数称为近似值的误差限：
$$ x^*-\epsilon \le x \le x^*+\epsilon $$
也可以表示为$x=x^*\pm \epsilon$

相对误差$\epsilon^*$

绝对误差与真值之比，称为近似值的相对误差（为了增强可操作性）
$$ \epsilon_r(x)=\frac{\epsilon(x)}{x}=\frac{x-x^*}{x} $$
* 相对误差可以体现误差绝对大小与测量度相对大小关系，也就是说误差的绝对大小并不能决定精度
然而实际问题中x未知，因而常将近似值替换准确值来计算相对误差，衡量准确程度：
$$ \epsilon^*_r(x)=\frac{\epsilon(x)}{x^*}=\frac{x-x^*}{x^*} $$

相对误差限

类似于绝对误差限，满足下列不等式的正数$\epsilon_r$为近似值的相对误差限
$$ \epsilon^*_r(x)=|\frac{\epsilon(x)}{x^*}|=|\frac{x-x^*}{x^*}| \le \epsilon_r $$

有效数字

近似值的误差限是其某一位上的半个单位时，称其“准确”到这意味，且从该位其到前面第一位非零数字为止所有数字都称为有效数字。
一般说，设有一个数x，其近似值的规格化形式为：
$$ x^*=\pm0.\alpha_1\alpha_2···\alpha_n×10^m$$
其中$\alpha_n$是0~9的数字,$\alpha_1\ne0$，n是正整数，m是整数。

绝对误差限

若$x^*$的误差限为：
$$ |\epslion(x)=|x-x^*|\ne\frac{1}{2}×10^{m-n} $$
则称$x^*$具有n位有效数字，或称其精确到$10^{m-n}$

注意：

  1. 绝对误差有单位，相对误差由于做了除法是没有单位的，常用百分号表示
  2. 四舍五入时，绝对误差限不超过末位数的半个单位
  3. 有效数字：
  • 203和0.203有效数字均为n=3，只是m不同
  • 为什么强调采用标准式写法：9000的有效数字位数不一定为4，我们无法确定这是测定准确结果还是仅将0用于占位
  1. 当计算有效数字位数和绝对误差限时，建议写作标准形式，写出m，n，套公式计算。

有效数字与误差的关系

相对误差回推有效数字；
$$ |\epsilon^*(x)|=|frac{\epsilon(x)}{x^*}\lefrac{feac{1}{2}×10^{m-n}}{\alpha_1×10^{m-1}}=frac{1}{2\alpha_1}×10^{-n+1} $$

相对误差限

从相对误差限的定义式可知，有效数字的位数反映了近似值的相对精确度：
$$ \alpha_r=frac{1}{2\alpha_1}×10^{-n+1} $$
当用相对误差表示时，如果相对误差能满足$$|\epsilon_r(x)| \le frac{1}{2(\alpha_1+1)}×10^{-n+1} $$
则$x^*$至少具有n位有效数字

数值运算中应注意的几个原则

  1. 选用数值稳定性好的算法
  • 什么是稳定性好
  • 如何判断一种算法具有良好的稳定性？
    假定数值具有初始误差$\epsilon_0$，中间不再产生新误差，考察由$\epsilon_0$引起的误差累计是否会增长，如果不增长则认为是稳定的，如果严重增长则认为不稳定
  • 如果改变为一种正确的算法，或许相同的计算类型稳定性可以增强（比如说改变递推方式，可以得到越来越精确的值）
  • 重要的是误差传递的方式
  • 可以通过误差传递方式由最后的误差倒退回误差初始值$\epsilon_0$
  1. 避免两个（有效位数）相近的数相减（会严重消损有效数字）
  • 在数值计算中两个相近的数相减会造成有效数字的严重损失，有一些好的解决方案：
    • 分母“无理化”
    • 降幂公式逆运算$1-cosx=s(sin(x/2))^2$
    • $arctan\theta_1-arctan\theta_2 = arctan(\frac{1}{\theta_1 \theta_2})$
    • x很小的时候泰勒展开计算$e^x$
  • 总结：无法改变算式时，可采用增加有效数字位数进行运算（计算机上采用双精度计算，但是会增大时空开销）

有效数字的运算规则

  1. 加减运算
    • 取舍以小数点后位数最少的数值为准
    • 先按照小数点后位数最少的数据保留其他各数的位数，再进行加减计算，结果也应当使小数点后保留相同的位数
  2. 乘除运算
    • 以有效数字最少的为准
    • 先按有效数字最少的数据保留其他各数，再进行乘除运算，计算结果仍保留相同有效数字
  3. 绝对值太小的不易做除数
    • 容易数据溢出
    • 小数点后高位略有误差会有较大影响
  4. 防止大数吃掉小数（Kahan算法，大数吃小数）
    • 字长有限，绝对值大小差距悬殊运算时容易出现
    • 小数进行一定累计之后再参与运算

范数

是一种映射，并且有相类似的规范

向量范数

三角不等式中是矢量相加
Img
Img
L2范数实际上是向量空间中的距离，无穷范数取最大值,p范数其实是其他特殊范数的一般形式；

矩阵范数

相对于向量范数，多了相容性；
求范数之后成为标量；
Img
根据向量的常用范数得到常用的矩阵范数：
Img

  • 1范数：每列和的绝对值取最大值
  • 无穷范数（行范数）：每行和的绝对值取最大值
  • 2范数实际上是欧氏空间中对距离的测定，最常用；ATA的正负不需要考虑，一定为正；
    Img
    2范数实际上是ATA与谱半径积的平方根
    单位矩阵任何范数都是1

扰动分析

病态方程与良态方程

对于线性方程组Ax=b，若系数矩阵A或常数项b的元素

讨论：常数项与矩阵的扰动对方程组解的影响

进行数据存在扰动时的误差分析，这里主要使用了矩阵范数的相容性性质，对于任何范数进行证明中的运算均成立

  1. 常数项的扰动
    Img

A二范数的条件数是A最大谱半径与最小谱半径商的平方根
根据A矩阵的条件数，可以确定解扰动的范围，体现了方程数据扰动对解相对误差的影响，初步判断是否病态。
其中A是扰动误差产生的核心，也是计算软件进行运算与否的重要依据之一

线性方程直接求解

两类解法：直接法与迭代法

直接法

高斯消去法

先消元，即按一定的规律逐步消去未知量，将方程组化为等价的上三角形方程组；然后进行回代，即由上三角形方程组逐个求出每个未知数的解

顺序高斯消去法

消元过程按方程和未知量的自然顺序进行：方程排列顺序，未知量位置顺序，类似于上三角形矩阵的变换；

  • 顺序高斯消元法的规则：
    • 依从左到右、自上而下的次序将主对角元下方的元素化为零。
    • 不作行交换，也不用非零数乘某行。
    • 第k列消元（将主对角元下方化零）时，将下面各行分别减去第k行的适当倍数，不作其他变换。
  • 即便是非奇异矩阵，也有可能是一整行都是0，或者出现一个较其他行相对较小系数，会造成严重的精度问题（数值稳定性差）
  • 使用上角标$^($)表示迭代次数
  • 行变换不会改变矩阵的非奇异性，理论上不会有某一行为0
  • 消去过程的时间复杂度更高
    Img
  • 消元过程可以按顺序进行的充要条件？
    • $a_{kk}^{(k-1)}\ne0,k=1,2···$
  • 递推Img
    Img

定理要会用：对于方程组：

  • 顺序消元过程能进行到底的充要条件是系数矩阵A的顺序主子矩阵Ak (k=1,2,…,n-1)非奇异。
  • 用顺序高斯消去法求解的充要条件是系数矩阵的一切顺序主子矩阵均非奇异。
  • 一些小补充

    顺序高斯消去法求解和顺序消去可以进行到底有什么区别?
    顺序高斯消去法（Sequential Gaussian Elimination）和顺序消元法（Sequential Elimination）都是线性代数中用于解决线性方程组的方法，但它们有一些区别。

  1. 顺序高斯消去法（Sequential Gaussian Elimination）：
    • 这是高斯消元法的一种变体，通过一系列的行变换将方程组化为上三角形式，然后通过回代求解方程组。在顺序高斯消去法中，通常会按顺序处理每一行，逐步进行消元操作。
    • 顺序高斯消去法通常用于解决稠密的线性方程组，因为它需要进行大量的矩阵操作。
    • 它的计算复杂度为 O(n^3)，其中 n 是方程组的大小。
  2. 顺序消元法（Sequential Elimination）：
    • 顺序消元法也是一种线性方程组求解方法，但不同于高斯消元法，它更侧重于直接消去变量而不是通过矩阵操作。
    • 在顺序消元法中，可以按照任意顺序依次消去方程组中的变量，并解出剩余的变量。通常采用的顺序是从第一个变量开始，逐步向后推进。
    • 顺序消元法通常适用于稀疏的线性方程组，因为它的操作更侧重于变量之间的消去，而不是整个矩阵的操作。
    • 它的计算复杂度取决于消元的顺序，通常为 O(n^2) 到 O(n^3)。

因此，两种方法的主要区别在于其操作的对象和计算复杂度。顺序高斯消去法更侧重于矩阵操作，适用于稠密的线性方程组，而顺序消元法更侧重于变量之间的消去，适用于稀疏的线性方程组。

选主元高斯消去法

Img

  • 基本思想：每次消元前按一定的范围选取绝对值最大的元素作为主元素。以便减少舍入误差的影响。

    • 主要有列主元高斯消去法和全主元高斯消去法.一般只使用行交换进行消去，可以降低复杂度
    • 只要detA≠０，就必能找到k行以下该列某不为0的元素
    • 若r>k 则交换第k行和第r行，然后进行消元。
      • 意思就是谁大把谁放入当前操作行作为保留元素
    • 不影响x，只是对排序调整
      Img
      最后得方程组，有-1，1，5
      Img
      Img

  • 一般使用列主元就已经足够

  • 进阶的选主元消去

    • 对算法的说明
    • 运算量：$\frac{n^3}{3}$(计算机中做乘除运算的时间远远超过做加减运算时间，故我们只估计 乘除运算 的次数)

  • Gauss-Jordan消去法
    Img
    Img

    用途
    高斯约当消去法（Gauss-Jordan elimination）是高斯消去法的一种改进形式，用于求解线性方程组并求解矩阵的逆。与高斯消去法不同，高斯约当消去法一次性将矩阵化为阶梯形式，并最终化为单位矩阵。因此，它能够同时解决线性方程组和求解逆矩阵的问题。
    高斯约当消去法的主要用途包括：

    1. 求解线性方程组：通过高斯约当消去法，可以将线性方程组的系数矩阵化为阶梯形式，并最终化为单位矩阵。这使得方程组的解可以直接从化简后的矩阵中读取。
    2. 求解逆矩阵：在高斯约当消去法中，将待求的矩阵与单位矩阵连接在一起，通过一系列的行变换将待求矩阵化为单位矩阵，同时单位矩阵也会变为原矩阵的逆。因此，高斯约当消去法也可用于求解矩阵的逆。
    3. 求解行列式：通过高斯约当消去法，可以化简原矩阵为上三角形式，从而可以直接计算矩阵的行列式。
      总之，高斯约当消去法是一种功能强大的方法，可以同时解决线性方程组、求解逆矩阵以及计算行列式等问题。

求方程解的例子

考虑以下线性方程组：
[
\begin{cases}
2x + 3y - z = 7 \
3x + y + 2z = 5 \
x - 2y + z = -1
\end{cases}
]
可以将这个线性方程组表示为增广矩阵形式：
[
\left[\begin{array}{ccc|c}
2 & 3 & -1 & 7 \
3 & 1 & 2 & 5 \
1 & -2 & 1 & -1
\end{array}\right]
]
现在，我们将通过高斯约当消元法将这个增广矩阵化简为阶梯形式。然后，从阶梯形式中读取方程组的解。
1. 第一步：将第一行第一个元素（主元素）变为 1，同时使用第一行的主元素消去其他行的对应元素。
[
\left[\begin{array}{ccc|c}
1 & \frac{3}{2} & -\frac{1}{2} & \frac{7}{2} \
3 & 1 & 2 & 5 \
1 & -2 & 1 & -1
\end{array}\right]
]
2. 第二步：使用第一行的第一个元素将第二行的第一个元素变为 0，并将第三行的第一个元素变为 0。
[
\left[\begin{array}{ccc|c}
1 & \frac{3}{2} & -\frac{1}{2} & \frac{7}{2} \
0 & -\frac{7}{2} & \frac{7}{2} & -\frac{11}{2} \
0 & -\frac{7}{2} & \frac{3}{2} & -\frac{9}{2}
\end{array}\right]
]
3. 第三步：将第二行第二个元素（主元素）变为 1，并使用第二行的主元素消去第一行和第三行的对应元素。
[
\left[\begin{array}{ccc|c}
1 & \frac{3}{2} & -\frac{1}{2} & \frac{7}{2} \
0 & 1 & -1 & 1 \
0 & -\frac{7}{2} & \frac{3}{2} & -\frac{9}{2}
\end{array}\right]
]
4. 第四步：使用第二行的第二个元素将第一行的第二个元素变为 0，并使用第二行的第二个元素将第三行的第二个元素变为 0。
[
\left[\begin{array}{ccc|c}
1 & 0 & 1 & 3 \
0 & 1 & -1 & 1 \
0 & 0 & 2 & 4
\end{array}\right]
]
5. 第五步：将第三行的第三个元素（主元素）变为 1，并使用第三行的主元素消去第一行和第二行的对应元素。
[
\left[\begin{array}{ccc|c}
1 & 0 & 0 & 1 \
0 & 1 & 0 & 2 \
0 & 0 & 1 & 2
\end{array}\right]
]
现在，我们得到了阶梯形式的增广矩阵。根据阶梯形式，我们可以直接读出方程组的解：
[
\begin{cases}
x = 1 \
y = 2 \
z = 2
\end{cases}
]
这就是线性方程组的解，通过高斯约当消元法得到。

重点&小结

Img

  1. 解线性方程组Gauss消去法的计算程序如何实现？可以进行Gauss消去法的条件是什么？

    1
    2
    3
    4
    5
    6
    7
    8
    9
    10
    11
    12
    13
    14
    15
    16
    17
    18
    19
    20
    21
    22
    23
    24
    25
    26
    27
    28
    29
    30
    31
    32
    33
    34
    35
    36
    37
    38
    39
    40
    41
    def gauss_elimination(A, b):
        n = len(A)
        
        # 合并增广矩阵
        for i in range(n):
            A[i].append(b[i])
        
        # 消元过程
        for i in range(n):
            # 首先找到主元素所在行
            max_index = i
            for j in range(i + 1, n):
                if abs(A[j][i]) > abs(A[max_index][i]):
                    max_index = j
            
            # 将主元素所在行交换至当前行
            A[i], A[max_index] = A[max_index], A[i]
            
            # 将当前主元素所在列的其他元素消为0
            for j in range(i + 1, n):
                factor = A[j][i] / A[i][i]
                for k in range(i, n + 1):
                    A[j][k] -= factor * A[i][k]
        
        # 回代求解
        x = [0] * n
        for i in range(n - 1, -1, -1):
            x[i] = A[i][n] / A[i][i]
            for j in range(i):
                A[j][n] -= A[j][i] * x[i]
        
        return x

    # 例子
    A = [[2, 1, -1],
         [-3, -1, 2],
         [-2, 1, 2]]
    b = [8, -11, -3]

    solution = gauss_elimination(A, b)
    print("Solution:", solution)

    Gauss 消去法可以应用于求解线性方程组的条件包括：

  2. 方程组的系数矩阵是一个非奇异矩阵：这意味着方程组的系数矩阵的行列式不为零。非奇异矩阵保证了方程组有唯一解。

  3. 方程组的系数矩阵是一个方阵：Gauss 消去法通常应用于求解方阵形式的线性方程组。如果方程组的系数矩阵不是方阵，那么方程组可能有无穷解或者无解。

  4. 可以按顺序进行消去的充要条件：主对角线元素不为零：Gauss 消去法在进行消元的过程中会涉及到主对角线元素的除法操作，因此要求主对角线元素不为零，避免除以零的情况。

  5. 方程组的条件数较小：条件数反映了方程组的稳定性，较小的条件数意味着方程组求解的稳定性较高，更容易得到准确的解。高条件数可能导致数值不稳定性和舍入误差的累积。

  6. 解线性方程组列主元Gauss消去法的计算程序如何实现？可以进行列主元Gauss消去法的条件是什么？
    条件同上；顺序消元过程进行到底&高斯消去法求解<->系数矩阵一切顺序主子矩阵均非奇异

    1
    2
    3
    4
    5
    6
    7
    8
    9
    10
    11
    12
    13
    14
    15
    16
    17
    18
    19
    20
    21
    22
    23
    24
    25
    26
    27
    28
    29
    30
    31
    32
    33
    34
    35
    36
    37
    38
    39
    40
    41
    42
    43
    44
    45
    def swap_rows(A, i, j):
        """交换矩阵的两行"""
        A[i], A[j] = A[j], A[i]

    def gauss_elimination(A, b):
        n = len(A)
        
        # 合并增广矩阵
        for i in range(n):
            A[i].append(b[i])
        
        # 消元过程
        for i in range(n):
            # 在当前列中找到绝对值最大的元素所在的行
            max_index = i
            for j in range(i + 1, n):
                if abs(A[j][i]) > abs(A[max_index][i]):
                    max_index = j
            
            # 将找到的主元素所在行交换至当前行
            swap_rows(A, i, max_index)
            
            # 将当前主元素所在列的其他元素消为0
            for j in range(i + 1, n):
                factor = A[j][i] / A[i][i]
                for k in range(i, n + 1):
                    A[j][k] -= factor * A[i][k]
        
        # 回代求解
        x = [0] * n
        for i in range(n - 1, -1, -1):
            x[i] = A[i][n] / A[i][i]
            for j in range(i):
                A[j][n] -= A[j][i] * x[i]
        
        return x

    # 例子
    A = [[2, 1, -1],
         [-3, -1, 2],
         [-2, 1, 2]]
    b = [8, -11, -3]

    solution = gauss_elimination(A, b)
    print("Solution:", solution)

  7. Gauss消去法的计算复杂度是？
    $O(n^3)$，其中 n 是方程组的未知数个数。

直接三角分解法

对于方程组AX=b，只要对系数矩阵作了三角分解A=LU，就可以很容易地求出线性方程的解。
Img

  • A的行列式就是U的行列式:
    $$ A_k=L_kU_k,detA_k=detU_k=\prod^k_{i=1} a^{(i-1)}_{ii} $$

会计算Dolittle（4，5）

核心：使用尽可能少的计算量写出非零数：通过矩阵乘法规则可以推出来A在LU上的“映射”，每行依次进行，先行后列，可以直接把LU中的非0数写出来
步骤如下：

  • L矩阵对角线元素全为1，U矩阵第一行元素与A第一行相同
  • 取k=1,2,…,n.依次计算k行U和k列L
    • $a_ij$ = Li行*Uj列,按这个来列方程解值
      Img
  • 求解方程组Ax＝b 的Doolittle 三角分解方法,包括分解和回代两步。
    Img
    Img