数值计算day4-求解线性方程组(下)

作者: xkzhai | 来源:发表于2019-07-28 10:34 被阅读0次

数值计算day4-求解线性方程组(下)
[C数值算法].pdf
第二章：Python数据分析简介
Conjugate gradient method
非线性方程组求解
数值计算day3-求解线性方程组(上)
2020-04-18
基本算法思想之概率
数值分析：线性方程组迭代求解
数值分析：非线性方程组求解

上节课主要介绍了线性方程组的几种直接求解法，包括高斯消去法（主元消去）、高斯约当法(可以用来求解矩阵的逆)、LU分解法（高斯消去和Crout 方法），使用高斯消去法实现LU分解的原理可以通过一种拉格朗日形式来理解。本节课主要介绍线性方程组的几种迭代求解法，包括Jacobi法和高斯塞德尔法。

1. 使用MATLAB内置函数实现LU分解以及矩阵求逆

LU分解

MATLAB使用lu命令对矩阵 $A$ 进行LU分解，分解前对矩阵 $A$ 进行行变换，以满足主元消去的条件：
$P A= LU$ $A=(P^{-1}L)U$

>> help lu
lu - lu 矩阵分解

    此 MATLAB 函数 返回矩阵 Y，其中包含严格下三角矩阵 L（即不带其单位对角线）和上三角矩阵 U 作为子矩阵。也即，如果 [L,U,P] =
    lu(A)，Y = U+L-eye(size(A))。不会返回置换矩阵 P。

    Y = lu(A)
    [L,U] = lu(A)
    [L,U,P] = lu(A)
    [L,U,P,Q] = lu(A)
    [L,U,P,Q,R] = lu(A)
    [...] = lu(A,'vector')
    [...] = lu(A,thresh)
    [...] = lu(A,thresh,'vector')

    另请参阅 cond, det, ilu, inv, qr, rref

    lu 的参考页
    名为 lu 的其他函数

>> A = [1 2 3;4 5 6;7 8 9]

A =

     1     2     3
     4     5     6
     7     8     9

>> format long 
>> det(A)

ans =

    -9.516197353929915e-16

>> [l u p] = lu(A)

l =

   1.000000000000000                   0                   0
   0.142857142857143   1.000000000000000                   0
   0.571428571428571   0.500000000000000   1.000000000000000


u =

   7.000000000000000   8.000000000000000   9.000000000000000
                   0   0.857142857142857   1.714285714285714
                   0                   0  -0.000000000000000


p =

     0     0     1
     1     0     0
     0     1     0

>> p*A

ans =

     7     8     9
     1     2     3
     4     5     6

>> l*u

ans =

     7     8     9
     1     2     3
     4     5     6

>> [l_1,u_1] = lu(A);
>> inv(p)*l == l_1

ans =

  3×3 logical 数组

   1   1   1
   1   1   1
   1   1   1

>> [l_2,u_2,p_2] = lu(p*A)

l_2 =

   1.000000000000000                   0                   0
   0.142857142857143   1.000000000000000                   0
   0.571428571428571   0.500000000000000   1.000000000000000


u_2 =

   7.000000000000000   8.000000000000000   9.000000000000000
                   0   0.857142857142857   1.714285714285714
                   0                   0  -0.000000000000000


p_2 =

     1     0     0
     0     1     0
     0     0     1

矩阵求逆

MATLAB使用inv命令求矩阵的逆，但当矩阵的数量级很小时，inv求解会带来较大的舍入误差（Round-off error）:

>> A = 10^(-4)*[1 2 3;4 5 6;7 8 9]

A =

   1.0e-03 *

   0.100000000000000   0.200000000000000   0.300000000000000
   0.400000000000000   0.500000000000000   0.600000000000000
   0.700000000000000   0.800000000000000   0.900000000000000

>> inv(A)
警告: 矩阵接近奇异值，或者缩放错误。结果可能不准确。RCOND =  1.317607e-17。 

ans =

   1.0e+19 *

  -0.527049830677416   1.054099661354831  -0.527049830677416
   1.054099661354831  -2.108199322709663   1.054099661354832
  -0.527049830677415   1.054099661354831  -0.527049830677416

>> A*inv(A)
警告: 矩阵接近奇异值，或者缩放错误。结果可能不准确。RCOND =  1.317607e-17。 

ans =

   1.250000000000000                   0                   0
   0.500000000000000   1.000000000000000  -0.500000000000000
   1.000000000000000                   0                   0

2. Jacobi法求解线性方程组

Jacobi法是一种同步迭代求解法，考虑如下线性方程构成的系统： $\begin{cases}a_{11}x_1+\cdots+a_{1n}x_n = b_1\\a_{21}x_1+\cdots+a_{2n}x_n = b_2\\ \vdots \\ a_{n1}x_1+\cdots+a_{nn}x_n = b_n\end{cases}$ 简单运算，可推导出：

$x_1= \frac{b_1-a_{12}x_2-...-a_{1n}x_n}{a_{11}}$ $x_2 = \frac{b_2-a_{21}-a_{23}x_3-...-a_{2n}x_n}{a_{22}}$ $\vdots$ $x_n = \frac{b_n-a_{n1}x_1-...-a_{n,n-1}x_{n-1}}{a_{n,n}}$ 这一转换过程可以描述为如下矩阵形式： $Ax=b\rightarrow x = Cx+D$ 其中矩阵 $C$ 和 $D$ 的元素为 $C_{ij} =\begin{cases} \frac{a_{ij}}{a_{ii}},\ i\neq j \\ 0,\ i=j\end{cases},D_{i} = \frac{b_i}{a_{ii}}$ Jacobi迭代法就是给出一个初始估计值 $x^{(1)}$ ，使用 $x^{(k+1)}=Cx^{(k)}+D$ 进行迭代求解，特点是每次迭代，同时更新 $x^{(k)}$ 的所有元素（同线性回归与逻辑回归中的梯度下降法）： $x^{(k+1)}_{i}=\frac{1}{a_{ii}}(b_i-\sum^n_{j=1,\ j\neq i}a_{ij}x^{(k)}_j)$ 通常初始值可以选择所有值为 $0$ 。
接下来，讨论该迭代法的收敛性条件。Jacobi迭代法收敛的一个充分条件为： $|a_{ii}|>\sum^n_{j=1,\ j\neq i}|a_{ij}|,i=1,2,...,n$ 满足这一条件的矩阵称为对角占优的(diagonally dominant)。但是，不满足这一条件时，Jacobi迭代法也有可能会收敛，这只是一个充分条件。
Jacobi迭代法的停止准则采取估计值的相对误差来定义 $|\frac{x^{(k+1)}_{i}-x^{(k)}_i}{x^{(k)}_i}|<\epsilon,i=1,2,...,n$ 其中 $\epsilon>0$ 是预先给定的一个精度。

3. 高斯-塞德尔法求解线性方程组

高斯塞德尔法与Jacobi法的推导过程是一致的，不同的是高斯塞德尔法使用异步迭代的方式进行迭代求解： $x^{(k+1)}_{1}=\frac{1}{a_{11}}[b_1-\sum^{n}_{j=2}a_{1j}x^{(k)}_j]$ $x^{(k+1)}_{i}=\frac{1}{a_{ii}}[b_i-\sum^{i-1}_{j=1}a_{ij}x^{(k+1)}_{j}-\sum^n_{j=i+1}a_{ij}x^{(k)}_j],i=2,3,....,n-1$ $x^{(k+1)}_n=\frac{1}{a_{nn}}[b_n-\sum^{n-1}_{j=1}a_{nj}x^{(k+1)}_{j}]$ 可以看到，第 $k+1$ 步的 $x^{(k+1)}_{i}$ 的值会使用 $k+1$ 步计算出来的 $x^{(k+1)}_{j}(j<i)$ 的值以及 $k$ 步计算的 $x^{(k)}_{j}(j>i)$ 值进行计算。高斯塞德尔的停止准则与Jacobi方法一致，一般情况下高斯塞德尔方法的收敛速度比Jacobi法快，且占用内存小。
例： $9x_1-2x_2+3x_3+2x_4 = 54.5$ $2x_1+8x_2-2x_3+3x_4 = -14$ $-3x_1+2x_2+11x_3-4x_4 = 12.5$ $-2x_1+3x_2+2x_3+10x_4 = -21$ MATLAB实现：

k = 1; x1 = 0; x2 = 0; x3 = 0; x4 = 0;
disp('  k         x1        x2        x3         x4')
fprintf(' %2.0f       %-8.5f  %-8.5f  %-8.5f  %-8.5f \n', k, x1, x2, x3, x4)
for k=2 : 8
    x1 = (54.5 - (-2*x2 + 3*x3 + 2*x4))/9;
    x2 = (-14 - (2*x1 - 2*x3 + 3*x4))/8;
    x3 = (12.5 - (-3*x1 + 2*x2 - 4*x4))/11;
    x4 = (-21 - (-2*x1 + 3*x2 + 2*x3))/10;
    fprintf(' %2.0f       %-8.5f  %-8.5f  %-8.5f  %-8.5f \n', k, x1, x2, x3, x4)
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
>> Program4_8
  k         x1        x2        x3         x4
  1       0.00000   0.00000   0.00000   0.00000  
  2       6.05556   -3.26389  3.38131   -0.58598 
  3       4.33336   -1.76827  2.42661   -1.18817 
  4       5.11778   -1.97723  2.45956   -0.97519 
  5       5.01303   -2.02267  2.51670   -0.99393 
  6       4.98805   -1.99511  2.49806   -1.00347 
  7       5.00250   -1.99981  2.49939   -0.99943 
  8       5.00012   -2.00040  2.50031   -0.99992

4. 三对角系统(Tridiagonal systems)与Thomas算法

三对角系统： $\begin{bmatrix}A_{11} & A_{12} & 0 & 0 & 0 & 0 & 0 & 0\\ A_{21} & A_{22} & A_{23} & 0 & 0 & 0 & 0 & 0\\ 0 & A_{32} & A_{33} & A_{34} & 0 &0 & 0 & 0\\ & & \cdots & \cdots & \cdots & & &\\ & & & \cdots & \cdots & \cdots & & \\0 & 0 & 0 & 0 A_{n-2,n-3} & A_{n-2,n-2} & A_{n-2,n-1} & 0\\ 0 & 0& 0& 0 & 0 & A_{n-1,n-2} & A_{n-1,n-1} & A_{n-1,n}\\0 & 0 & 0 & 0& 0 & 0 & A_{n,n-1} & A_{n,n} \end{bmatrix}\begin{bmatrix}x_1\\x_2\\x_3\\ \vdots \\ \vdots \\ x_{n-2} \\ x_{n-1} \\ x_{n}\end{bmatrix}=\begin{bmatrix}B_1\\B_2\\B_3\\ \vdots \\ \vdots \\ B_{n-2} \\ B_{n-1} \\ B_{n}\end{bmatrix}$ 这样的系统可以通过高斯法，高斯约当法以及LU分解法求解，但基于系统矩阵的特殊结构（Sparse matrix, 稀疏矩阵），可以使用另一种更节省内存和时间的方法——Thomas算法。
Thomas算法将三角矩阵的对角线元素存储到一个向量 $d$ 中： $d=\begin{bmatrix}A_{11} & A_{22} & \cdots & A_{nn}\end{bmatrix}^T$ ，上对角元素存储到向量 $a$ 中： $a=\begin{bmatrix}A_{12} & A_{23} & \cdots & A_{n-1,n} & 0\end{bmatrix}^T$ ，下对角元素存储到向量 $b$ 中： $b=\begin{bmatrix}0 & A_{21} & A_{32} & \cdots & A_{n,n-1}\end{bmatrix}^T$ ，上述三对角系统可以转换为： $\begin{bmatrix}d_{1} & a_{1} & 0 & 0 & 0 & 0 & 0 & 0\\ b_{2} & d_{2} & a_{2} & 0 & 0 & 0 & 0 & 0\\ 0 & b_{3} & d_{3} & a_{3} & 0 &0 & 0 & 0\\ & & \cdots & \cdots & \cdots & & &\\ & & & \cdots & \cdots & \cdots & & \\0 & 0 & 0 & 0 & b_{n-2} & d_{n-2} & a_{n-2} & 0\\ 0 & 0& 0& 0 & 0 & b_{n-1} & d_{n-1} & a_{n-1}\\0 & 0 & 0 & 0& 0 & 0 & b_{n} & d_{n} \end{bmatrix}\begin{bmatrix}x_1\\x_2\\x_3\\ \vdots \\ \vdots \\ x_{n-2} \\ x_{n-1} \\ x_{n}\end{bmatrix}=\begin{bmatrix}B_1\\B_2\\B_3\\ \vdots \\ \vdots \\ B_{n-2} \\ B_{n-1} \\B_{n}\end{bmatrix}$ Thomas算法的第一步是将第一行元素规范化，得到： $\begin{bmatrix}1 & a'_{1} & 0 & 0 & 0 & 0 & 0 & 0\\ b_{2} & d_{2} & a_{2} & 0 & 0 & 0 & 0 & 0\\ 0 & b_{3} & d_{3} & a_{3} & 0 &0 & 0 & 0\\ & & \cdots & \cdots & \cdots & & &\\ & & & \cdots & \cdots & \cdots & & \\0 & 0 & 0 & 0 & b_{n-2} & d_{n-2} & a_{n-2} & 0\\ 0 & 0& 0& 0 & 0 & b_{n-1} & d_{n-1} & a_{n-1}\\0 & 0 & 0 & 0& 0 & 0 & b_{n} & d_{n} \end{bmatrix}\begin{bmatrix}x_1\\x_2\\x_3\\ \vdots \\ \vdots \\ x_{n-2} \\ x_{n-1} \\ x_{n}\end{bmatrix}=\begin{bmatrix}B'_1\\B_2\\B_3\\ \vdots \\ \vdots \\ B_{n-2} \\ B_{n-1} \\B_{n}\end{bmatrix}$ 其中 $a'_{1}=a_1/d_1,B'_1=B_1/d_1$ ，消去 $b_2$ ，得 $\begin{bmatrix}1 & a'_{1} & 0 & 0 & 0 & 0 & 0 & 0\\0 & d'_{2} & a_{2} & 0 & 0 & 0 & 0 & 0\\ 0 & b_{3} & d_{3} & a_{3} & 0 &0 & 0 & 0\\ & & \cdots & \cdots & \cdots & & &\\ & & & \cdots & \cdots & \cdots & & \\0 & 0 & 0 & 0 & b_{n-2} & d_{n-2} & a_{n-2} & 0\\ 0 & 0& 0& 0 & 0 & b_{n-1} & d_{n-1} & a_{n-1}\\0 & 0 & 0 & 0& 0 & 0 & b_{n} & d_{n} \end{bmatrix}\begin{bmatrix}x_1\\x_2\\x_3\\ \vdots \\ \vdots \\ x_{n-2} \\ x_{n-1} \\ x_{n}\end{bmatrix}=\begin{bmatrix}B'_1\\B'_2\\B_3\\ \vdots \\ \vdots \\ B_{n-2} \\ B_{n-1} \\B_{n}\end{bmatrix}$ 其中 $d'_2=d_2-b_2a'_1,B'_2=B_2-B'_1b_2$ ，重复这一过程直到有 $\begin{bmatrix}1 & a'_{1} & 0 & 0 & 0 & 0 & 0 & 0\\0 & 1& a'_{2} & 0 & 0 & 0 & 0 & 0\\ 0 & 0 & 1 & a'_{3} & 0 &0 & 0 & 0\\ & & \cdots & \cdots & \cdots & & &\\ & & & \cdots & \cdots & \cdots & & \\0 & 0 & 0 & 0 & 0 & 1 & a'_{n-2} & 0\\ 0 & 0& 0& 0 & 0 & 0 & 1 & a'_{n-1}\\0 & 0 & 0 & 0& 0 & 0 & 0 & 1 \end{bmatrix}\begin{bmatrix}x_1\\x_2\\x_3\\ \vdots \\ \vdots \\ x_{n-2} \\ x_{n-1} \\ x_{n}\end{bmatrix}=\begin{bmatrix}B'_1\\B'_2\\B'_3\\ \vdots \\ \vdots \\ B'_{n-2} \\ B'_{n-1} \\B'_{n}\end{bmatrix}$ 这样一个系统很容易可以通过回代法求解。

Thomas算法的步骤可以归结为：

定义 $d=\begin{bmatrix}d_{1} & d_{2} & \cdots & d_{n}\end{bmatrix}^T$ , $b=\begin{bmatrix}0 & b_{2} & b_{3} & \cdots & b_{n}\end{bmatrix}^T$ , $a=\begin{bmatrix}a_{1} & a_{2} & \cdots & a_{n-1} & 0\end{bmatrix}^T$ , $B=\begin{bmatrix}B_{1} & B_{2} & \cdots & B_{n}\end{bmatrix}^T$ ,
更新 $a_1$ , $B_1$ ： $a_1=\frac{a_1}{d_1}$ , $B_1=\frac{B_1}{d_1}$
对 $i=2,3,...,n-1$ ： $a_i=\frac{a_i}{d_i-b_ia_{i-1}},\ B_i=\frac{B_i-b_iB_{i-1}}{d_i-b_ia_{i-1}}$
$B_n=\frac{B_{n}-b_nB_{n-1}}{d_n-b_na_{n-1}}$
使用回代法求解： $x_n=B_n$ , $x_i=B_i-a_ix_{i+1},i=n-1,n-2,...,1$

MATLAB实现：

function x = Tridiagonal(A,B)
% The function solve a tridiagonal system of linear equations [a][x]=[b]
% using Thomas algorithm.
% Input variables:
% A  The matrix of coefficients.
% B  A column vector of constants.
% Output variable:
% x  A colum vector with the solution.

[nR, nC] = size(A);
for i = 1:nR
    d(i) = A(i,i);
end
for i = 1:nR-1
    ad(i) = A(i,i+1);
end
for i = 2:nR
    bd(i) = A(i,i-1);
end
ad(1) = ad(1)/d(1);
B(1) = B(1)/d(1);
for i = 2:nR-1
    ad(i) = ad(i)/(d(i)-bd(i)*ad(i-1));
    B(i)=(B(i)-bd(i)*B(i-1))/(d(i)-bd(i)*ad(i-1));
end
B(nR)=(B(nR)-bd(nR)*B(nR-1))/(d(nR)-bd(nR)*ad(nR-1));
x(nR,1) = B(nR);
for i = nR-1:-1:1
    x(i,1) = B(i)-ad(i)*x(i+1);
end

5. 范数与条件数

5.1 向量范数

一般来说，向量范数需要满足一下三个条件：

正定性 $\|x\|\geq 0,\forall x\in R^{n}, \|x\|=0\text{当且仅当}x=0$
齐次性 $\|kx\|=k\|x\|,\forall k\in R, x\in R^{n}$
三角不等式 $\|x+y\|\leq \|x\|+\|y\|,\forall x,y\in R^{n}$

几种常见的向量范数：

1范数： $\|x\|_1=\sum^n_{i=1}|x_i|$
2范数： $\|x\|_2=\sqrt{x^2_1+...+x^2_n}$
$\infty$ 范数： $\|x\|_{\infty}=\max_{1\leq i\leq n}|x_i|$

5.2 矩阵 $p$ 范数

矩阵范数要比向量范数多一个条件： $\|AB\|\leq \|A\|\|B\|$
矩阵范数中的 $p$ 范数均是由向量范数诱导而来： $\|A\|_p=max\frac{\|Ax\|_p}{\|x\|_p}=\max_{\|x\|_p=1}\|Ax\|_p$

1范数 $\|A\|_1=\max_{\|x\|_1=1}\|Ax\|_1=\max_{1\leq i\leq n}\|Ae_i\|_1\text{（列绝对值和中的最大值）}$ 证明：对任意 $\|x\|_1=1$ ，注意 $x=k_1e_1+...+k_ne_n$ 因此 $|k_1|+...|k_n|=1$ 根据向量范数的三角不等式，可得 $\|Ax\|_1=\|k_1Ae_1+...+k_nAe_n\|_1\leq |k_1|\|Ae_1\|_1+...|k_n|\|Ae_n\|_1\leq (|k_1|+...+|k_n|)\max_{1\leq i\leq n}\|Ae_i\|_1=\max_{1\leq i\leq n}\|Ae_i\|_1$ 即 $\|A\|_1\leq \max_{1\leq i\leq n}\|Ae_i\|_1$ 当 $x=e_i$ 时，能够取等号，因此 $\|A\|_1= max_{1\leq i\leq n}\|Ae_i\|_1$
2范数 $\|A\|_2=\max_{\|x\|_2=1}\|Ax\|_2=\sqrt{\lambda_{max}(A^TA)}$ 证明：假设 $A^TA$ 的特征值为 $\lambda_1,...,\lambda_n\geq 0$ , 对应的单位特征向量为 $x^((1)),...,x^{(n)}\text{（该组特征向量为$ R^{n} $空间的一组单位正交基有待证明）}$ 对 $\|x\|_2=1,x=k_1x^{(1)}+...+k_{n}x^{(n)}$ ，有 $\sqrt{k^2_1+...+k^2_n}=1$ 因此 $\|Ax\|_2=\sqrt{x^TA^TAx}=\sqrt{(k_1x^{(1)}+...+k_{n}x^{(n)})^T(\lambda_1k_1x^{(1)}+...+\lambda_nk_{n}x^{(n)})}=\sqrt{\lambda_1k^2_1+...+\lambda_nk^2_n}\leq \sqrt{\lambda_{max}(A^TA)}$ 当 $x$ 取最大特征值对应的单位特征向量时，等号成立，因此 $\|A\|_2=\sqrt{\lambda_{max}(A^TA)}$
$\infty$ 范数 $\|A\|_{\infty}=\max_{\|x\|_{\infty}=1}\|Ax\|_{\infty}=\max_{1\leq i\leq n}\|A^Te_i\|_1\text{（行绝对值和中的最大值）}$ 证明： $x=k_1e_1+...k_ne_n$ 当 $\|x\|_{\infty}=1$ 时，有 $\max_{1\leq i\leq n}|k_i|=1$ ，因此 $\|Ax\|_{\infty}=\|k_1Ae_1+...+k_nAe_n\|_{\infty}\leq \|Ae_1+...+Ae_n\|_{\infty}=\|A^Te_i\|_1\text{（行绝对值和中的最大值）}$ 当 $x^T=[1 ... 1]$ 时，等号成立

5.3 条件数与相对误差

对线性系统 $Ax=b$ ，假设其真实解为 $x_E$ ，数值解为 $x_{NS}$ ，则真实误差为 $TrueError = x_E-x_{NS}$ ，一般真实解无法给出，因此，真实误差难以计算，但系统残差（residual）则可以通过下式给出： $R=Ax_E-Ax_{NS}=b-Ax_{NS}$
例： $\begin{cases}1.02x_1+0.98x_2=2\\0.98x_1+1.02x_2\end{cases}$ 该系统的真实解为 $x_E=\begin{bmatrix}1\\1\end{bmatrix}$

若数值解为 $x_{NS}=\begin{bmatrix}1.02\\1.02\end{bmatrix}$ 真实误差为 $TrueError=-\begin{bmatrix}0.02\\0.02\end{bmatrix}$ 系统残差为 $R=\begin{bmatrix}1.02 & 0.98\\0.98 & 1.02\end{bmatrix}\begin{bmatrix}-0.02\\-0.02\end{bmatrix}=\begin{bmatrix}-0.04\\-0.04\end{bmatrix}$ 系统残差与真实误差的数量级时一致的
若数值解为 $x_{NS}=\begin{bmatrix}2\\0\end{bmatrix}$ 真实误差为 $TrueError=\begin{bmatrix}-1\\1\end{bmatrix}$ 系统残差为 $R=\begin{bmatrix}1.02 & 0.98\\0.98 & 1.02\end{bmatrix}\begin{bmatrix}-1\\1\end{bmatrix}=\begin{bmatrix}-0.04\\0.04\end{bmatrix}$ 此时，系统残差的数量级与真实误差的数量级差别很大，因此无法总是用残差来估计真实误差

一个好的想法是用范数来衡量估计误差， $TrueError = x_E-x_{NS}$ ， $R=A(x_E-x_{NS})=A*TrueError$ ，注意 $TrueError=A^{-1}R$ ，根据范数的三角不等式，有 $\|TrueError\|\leq \|A^{-1}\|\|R\|$ 进一步考虑相对误差和相对残差 $\frac{\|TrueError\|}{\|x_E\|},\ \frac{\|R\|}{\|b\|}$ 有 $\frac{\|TrueError\|}{\|x_E\|}\leq \|A^{-1}\|*\frac{\|b\|}{\|x_E\|}*\frac{\|R\|}{\|b\|}$ 注意 $Ax_E=b,\ \|b\|\leq \|A\|\|x_E\|\rightarrow \|x_E\|\geq \frac{\|b\|}{\|A\|}\rightarrow \frac{1}{\|x_E\|}\leq \frac{\|A\|}{\|b\|}$ 因此有 $\frac{\|TrueError\|}{\|x_E\|}\leq \|A^{-1}\|\|A\|\frac{\|R\|}{\|b\|}$ 另一方面， $\|R\|\leq \|A\|\|TrueError\|\rightarrow \|TrueError\|\geq \frac{\|R\|}{\|A\|}$ 因此 $\frac{\|TrueError\|}{\|x_E\|}\geq \frac{1}{\|x_E\|}\frac{\|R\|}{\|A\|}\rightarrow \frac{\|TrueError\|}{\|x_E\|}\geq \frac{1}{\|A\|}\frac{\|b\|}{\|x_E\|}\frac{\|R\|}{\|b\|}$ 同时 $x_E=A^{-1}b\rightarrow \|x_E\|\leq \|A^{-1}\|\|b\|\rightarrow \frac{\|b\|}{\|x_E\|}\geq \frac{1}{\|A^{-1}\|}$ 因此有 $\frac{\|TrueError\|}{\|x_E\|}\geq \frac{1}{\|A\|}\frac{\|b\|}{\|x_E\|}\frac{\|R\|}{\|b\|}\geq \frac{1}{\|A^{-1}\|\|A\|}\frac{\|R\|}{\|b\|}$ 综上，有 $\frac{1}{\|A^{-1}\|\|A\|}\frac{\|R\|}{\|b\|}\leq \frac{\|TrueError\|}{\|x_E\|}\leq \|A^{-1}\|\|A\|\frac{\|R\|}{\|b\|}$ $\|A\|\|A^{-1}\|$ 就称作条件数(Condition number)。

例：考虑第3节的线性方程组 $9x_1-2x_2+3x_3+2x_4 = 54.5$ $2x_1+8x_2-2x_3+3x_4 = -14$ $-3x_1+2x_2+11x_3-4x_4 = 12.5$ $-2x_1+3x_2+2x_3+10x_4 = -21$ 其真实解为 $x_E=\begin{bmatrix}5 \\ -2 \\2.5\\-1\end{bmatrix}$ 而使用第3节的高斯赛德尔法迭代6次后，得到的数值解为 $x_{NS}=\begin{bmatrix}4.98805\\ -1.99511\\ 2.49806 \\ -1.00347 \end{bmatrix}$ 因此真实误差为 $TrueError=x_E-x_{NS}=\begin{bmatrix}0.0119\\-0.0049\\0.0019\\0.0035\end{bmatrix}$ 系统残差为 $R=\begin{bmatrix}0.13009\\ -0.00869\\-0.03817\\0.00001\end{bmatrix}$ 现在以向量 $\infty$ 为例，检查一下上述不等式的真实性： $\|x_E\|_{\infty}=5,\ \|TrueError\|_{\infty}=0.0119,\ \|R\|_{\infty} = 0.13009,\ \|b\|=54.5, \|A\|_{\infty}=20,\ \|A^{-1}\|_{\infty}=0.19019$ 注意 $\|A\|\|A^{-1}\|=3.8038\rightarrow \frac{1}{\|A^{-1}\|\|A\|}\frac{\|R\|}{\|b\|}=0.00062749 \leq \frac{\|TrueError\|}{\|x_E\|}=0.00238 \leq \|A^{-1}\|\|A\|\frac{\|R\|}{\|b\|}=0.009079$

病态系统（ill-conditioned）：病态系统是指系统参数有很小的变化时都会引起解很大的变化的一类系统，且可以证明当条件数远大于 $1$ 时，系统是一个病态的。

6. 总结

本节课主要介绍了求解线性方程组的两种迭代数值求解法，包括Jacobi法和高斯塞德尔法，其中Jacobi法是同步更新，而高斯赛德尔法是异步更新。对于特殊的三对角矩阵，可以采用Thomas算法求解，相比高斯消去、LU分解法更节省内存和时间。最后，介绍了向量范数和矩阵范数的概念，并引入条件数来确定线性系统数值解的相对误差范围。

数值计算day4-求解线性方程组(下)
上节课主要介绍了线性方程组的几种直接求解法，包括高斯消去法（主元消去）、高斯约当法(可以用来求解矩阵的逆)、LU分...
[C数值算法].pdf
【下载地址】本书编写了300多个实用而有效的数值算法C语言程序。其内容包括：线性方程组的求解，逆矩阵和行列式计算...
第二章：Python数据分析简介
2.1Scipy求解非线性方程组和数值积分 2.2Matplotlib作图基本代码如果使用的是中文标签，就会发现...
Conjugate gradient method
共轭梯度法。一种求解数学特定线性方程组Ax=b的数值解的迭代方法。要求矩阵A对称（symmetric）且正定（po...
非线性方程组求解
非线性方程组求解 fsolve() 可以对非线性方程组进行求解，它的基本调用形式为fsolve(func，x0)。...
数值计算day3-求解线性方程组(上)
上节课主要介绍了非线性方程的几种数值解法，其中包括交叉法（二分法、线性插值法）和开放法（牛顿法、割线法、固定点法）...
2020-04-18
有会用matlab求解线性方程组的吗，想请教一下，急求，真诚感谢。
基本算法思想之概率
因为很多数学问题，往往没有或者很难计算解析，这时候便需要通过数值计算来求解近似值.概率算法依照概率统计的思路来求解...
数值分析：线性方程组迭代求解
前言常用的线性方程组的直接法求解方法有LU分解、高斯消去、列主元消去法等，但是这些直接法最大的缺点就是对系数矩阵...
数值分析：非线性方程组求解
前言非线性方程组，顾名思义就是未知数的幂除了不是1，其他都有可能！线性方程组其实只是非常小的一类，非线性方程组才...