目录
考虑方程组=b的解,其中系数矩阵H为Hilbert矩阵:
这是一个著名的病态问题。通过首先给定解（例如取为各个分量均为1）再计算出右端的办法给出确定的问题。
(1）选择问题的维数为6，分别用 Gauss消去法、J迭代法、GS迭代法和SOR迭代法求解方程组,其各自的结果如何?将计算结果与问题的解比较,结论如何。
(2）逐步增大问题的维数，仍用上述的方法来解它们，计算的结果如何?计算的结果说明的什么?
(3）讨论病态问题求解的算法。
实验目的:通过本实验加深对于高斯消去法，Jacobi迭代法，G-S迭代法和SOR迭代法的理解和运用，同时理解不同阶数对于方程的影响，同时加深对于病态方程的理解。
假设给定一个线性方程组 AX=b，其中 A 是一个 n×n 的系数矩阵，X 和 b 分别是未知数向量和常数向量。高斯消去法的主要步骤如下：
1.构建增广矩阵：将系数矩阵 A 和常数向量 b 合并成一个增广矩阵 [A | b]。
2.消元过程：通过一系列的初等行变换将增广矩阵转化为上三角形矩阵。这些初等行变换包括以下操作：
3.交换两行的位置；
4.用一个非零常数乘以某一行的所有元素；
5.用一行的倍数加到另一行上。
在消元过程中，每一步都以一个主元为基准，将主元以下的列中的元素通过消元操作变为零，以实现逐步消去下方的元素。
回代过程：从上三角形矩阵的最后一行开始，通过回代求解未知数的值。回代的步骤如下：
求解最后一行的未知数；继续向上一行回代，依次求解每个未知数，直到求解出所有未知数。回代过程中，每个未知数的解是通过将已知的未知数代入方程进行求解得到的。
得到方程组的解：回代过程完成后，我们得到了线性方程组的解向量 X。
高斯消去法的关键思想是通过消元操作将线性方程组转化为上三角形矩阵，从而使回代过程更加简单和可行。它的优点是简单易实现，但它的缺点是可能会面临主元为零或接近零的情况，导致算法的稳定性下降，需要进行一些修正，如选取主元时使用行部分选主元（部分主元消去法）或使用列部分选主元（完全主元消去法）。
这里还有LU分解的过程：
对于
（A非奇异）求解时，可以先将A分解成一个下三角矩阵L和一个上三角矩阵U的乘积，即
，就可以通过
求解出
的值。
欲把一个给定的矩阵A分解为一个下三角阵L与一个上三角阵U的乘积，最自然的做法便是通过一系列的初等变换，逐步将A约化为一个上三角阵，而又能保证这些变换的乘积是一个下三角阵。这可归结为：对于一个任意给定的向量
找一个尽可能简单的下三角阵，使
经这一矩阵作用之后的第
至第
个分量均为零。能够完成这一任务的最简单的下三角阵便是如下形式的初等下三角阵：
其中
即
这种类型的初等下三角阵称作Gauss变换，而称向量
为Gauss向量。
对于一个给定的向量
我们有：
由此立即可知，只要取
便有
当然，这里我们要求
而后经过多次变换可以得到
从而求出上三角阵U，而后通过
求得下三角阵
将结果带入到
式中求出
的值即可。
Jacobi迭代法是一种用于解线性方程组的迭代方法，其原理基于将线性方程组的解逐步逼近的思想。
假设给定一个线性方程组 AX=b，其中 A 是一个 n×n 的系数矩阵，X 和 b 分别是未知数向量和常数向量。Jacobi迭代法的主要步骤如下：
1.将线性方程组转化为等价的形式：将 A 分解为 A = D - L - U，其中 D 是 A 的对角线矩阵，L 是 A 的严格下三角矩阵（即除去对角线的元素），U 是 A 的严格上三角矩阵（即除去对角线的元素）。
2.初始化：将解向量 X 的初始值设置为一个任意的向量，通常选择一个全零向量。
迭代过程：根据迭代公式 X(k+1) = D^(-1) * (b - (L + U) * X(k))，其中 X(k) 表示第 k 次迭代的解向量，X(k+1) 表示第 k+1 次迭代的解向量。
3.在每一次迭代中，根据上述迭代公式，将 X(k) 代入计算 X(k+1)。
4.迭代过程将 X 的值逐步逼近最终的解向量。
5.判断终止条件：根据需要设置一个终止条件，如达到一定的迭代次数或解向量的变化量小于某个阈值。
6.输出解向量：当满足终止条件时，迭代停止，输出最终的解向量 X。
以下是Jacobi迭代法的公式推导:
考虑非奇异线性代数方程组:
令
其中
那么合并之后可以写成一下形式:
其中
若给定初始向量
并代入(1.4)式右边，又可得到一个向量
；一次类推，有
这就是所谓的Jacobi迭代法，其中
叫做Jacobi迭代法的迭代矩阵，
叫做Jacobi迭代法的常数项。
Gauss-Seidel（G-S）迭代法是一种用于解线性方程组的迭代方法，其原理基于将线性方程组的解逐步逼近的思想，并结合了前向替代和后向替代。
假设给定一个线性方程组 AX=b，其中 A 是一个 n×n 的系数矩阵，X 和 b 分别是未知数向量和常数向量。G-S迭代法的主要步骤如下：
1. 将线性方程组转化为等价的形式：将 A 分解为 A = D - L - U，其中 D 是 A 的对角线矩阵，L 是 A 的严格下三角矩阵（即除去对角线的元素），U 是 A 的严格上三角矩阵（即除去对角线的元素）。
2. 初始化：将解向量 X 的初始值设置为一个任意的向量，通常选择一个全零向量。
3. 迭代过程：根据迭代公式 X(k+1) = D^(-1) * (b - L * X(k+1) - U * X(k))，其中 X(k) 表示第 k 次迭代的解向量，X(k+1) 表示第 k+1 次迭代的解向量。
在每一次迭代中，根据上述迭代公式，先使用已知的 X(k) 来计算新的 X(k+1) 的值，然后将新的 X(k+1) 代入计算下一次迭代的 X(k+1)。
在计算 X(k+1) 的过程中，使用的是当前已经计算出来的新值，这就是 G-S 迭代法与 Jacobi 迭代法的区别。
4. 判断终止条件：根据需要设置一个终止条件，如达到一定的迭代次数或解向量的变化量小于某个阈值。
5. 输出解向量：当满足终止条件时，迭代停止，输出最终的解向量 X。
G-S迭代法的关键思想是通过迭代过程，利用当前解向量的近似值来计算下一次迭代的解向量，逐步逼近方程组的解。相比于 Jacobi 迭代法，G-S 迭代法在每次迭代中直接使用新的解向量，因此收敛速度相对较快。然而，对于某些特殊的方程组，G-S 迭代法可能不收敛或收敛较慢。在实际应用中，可以通过调整迭代次数或使用其他更高效的迭代方法来提高求解效率。
以下是G-S迭代法的公式推导：
注意到Jacobi迭代法中各分量的计算顺序是没有关系的，先算那个分量都一样。现在，假设不按Jacobi迭代格式，而是在计算
的第一个分量用
的各个分量计算，但当计算
的第二个分量
时，因
已经算出，用它代替
，其他分量仍用
。类似地，计算
时，因
都已算出，用它们代替
其他分量仍用
的分量，于是有
我们称这种迭代格式为Gauss-Seidel迭代法，简称为G-S迭代法。它的一个明显的好处是在编写程序是存储量减少了。如果
存在，G-S迭代法可以改写成
我们把
叫做G-S迭代法的迭代矩阵，而把
叫做G-S迭代法的常数项。
Successive Over-Relaxation（SOR）迭代法是一种用于解线性方程组的迭代方法，其原理基于将线性方程组的解逐步逼近的思想，并结合了Gauss-Seidel迭代法和松弛因子的概念。
假设给定一个线性方程组 AX=b，其中 A 是一个 n×n 的系数矩阵，X 和 b 分别是未知数向量和常数向量。SOR迭代法的主要步骤如下：
1. 将线性方程组转化为等价的形式：将 A 分解为 A = D - L - U，其中 D 是 A 的对角线矩阵，L 是 A 的严格下三角矩阵（即除去对角线的元素），U 是 A 的严格上三角矩阵（即除去对角线的元素）。
2. 初始化：将解向量 X 的初始值设置为一个任意的向量，通常选择一个全零向量。
3. 迭代过程：根据迭代公式 X(k+1) = (1 - w) * X(k) + w * D^(-1) * (b - L * X(k+1) - U * X(k))，其中 X(k) 表示第 k 次迭代的解向量，X(k+1) 表示第 k+1 次迭代的解向量，w 是松弛因子（0 < w < 2）。
在每一次迭代中，根据上述迭代公式，首先使用已知的 X(k) 来计算新的 X(k+1) 的值，然后将新的 X(k+1) 代入计算下一次迭代的 X(k+1)。
松弛因子 w 控制着迭代的加速程度。当 w=1 时，SOR迭代法就等价于Gauss-Seidel迭代法；当 w<1 时，迭代过程为欠松弛；当 w>1 时，迭代过程为过松弛。
4.判断终止条件：根据需要设置一个终止条件，如达到一定的迭代次数或解向量的变化量小于某个阈值。
5. 输出解向量：当满足终止条件时，迭代停止，输出最终的解向量 X。
SOR迭代法的关键思想是通过迭代过程，利用当前解向量的近似值来计算下一次迭代的解向量，逐步逼近方程组的解。通过调整松弛因子 w，可以控制迭代的加速程度。在实际应用中，选择合适的松弛因子 w 对于迭代的收敛性和效率非常重要。如果选择合适的 w，SOR迭代法可以加快收敛速度。
以下是SOR迭代法的推导公式:
Q1：选择问题的维数为 6，分别用 Gauss 消去法、J迭代法、GS 选代法和SOR迭代法求解方程组,其各自的结果如何?将计算结果与问题的解比较,结论如何
（1）首先将实验的数据导入并且做成numpy矩阵：
# 构建 Hilbert 矩阵 H
H = np.zeros((n, n))
for i in range(1, n + 1):
for j in range(1, n + 1):
H[i-1, j-1] = 1 / (i + j - 1)
print(H)
# 构建解向量 x
x = np.ones(n)
# 计算右端向量 b
b = np.dot(H, x)
#高斯消元法
def gaussian_elimination(A, b):
n = len(A)
# 构建增广矩阵
aug_matrix = np.concatenate((A, b.reshape(-1, 1)), axis=1)
# 消元过程
for i in range(n):
pivot = aug_matrix[i, i] # 主元
# 消除下方的元素
for j in range(i + 1, n):
factor = aug_matrix[j, i]/pivot
aug_matrix[j, :] -= factor * aug_matrix[i, :]
return aug_matrix //返回消去最后的矩阵
# 调用高斯消去法求解方程组
aug_matrix = gaussian_elimination(H, b)
# 输出最后的 A 矩阵
UA= aug_matrix[:, :-1]
print("最后的 A 矩阵:")
for i in range(n):
for j in range(n):
if(i>j):
UA[i][j]=0
print(UA)
rUA = np.linalg.inv(UA)
LA=np.dot(H,rUA)
print("LA")
for i in range(n):
for j in range(n):
if(i<j):
LA[i][j]=0
print(LA)
rLA = np.linalg.inv(LA)
y=np.dot(rLA,b)
rUA = np.linalg.inv(UA)
x=np.dot(rUA,y)
print(x)
（程序运行结果1）高斯消去法:
计算返回的A矩阵的最后形式就是LU分解所得的U矩阵
再通过
得出矩阵L的值为:
最后通过Ly=b y=Ux两个公式的代换求出要求x的值为：
故该方程的解为:
首先计算L,U,D,B矩阵:
D = np.diag(np.diag(H))
D_inv = np.linalg.inv(D)
L = -np.tril(H, k=-1)
U = -np.triu(H, k=1)
B=np.dot(D_inv,(L+U))
print(D)
print(L)
print(U)
print(B)
D矩阵如下：
L矩阵
U矩阵:
计算B矩阵的值: 
计算B矩阵的谱半径判断其是否收敛:
得出
所以矩阵式发散的,Jacobi迭代法不能求解x的值。
首先先求
和 
两个矩阵的值:
C = np.linalg.inv(D-L)
L1=np.dot(C,U)
G1=np.dot(C,b)
print(L1)
print(G1)
的矩阵如下:
的矩阵如下:
接下来求谱半径的值:
(eva1, evt1) = np.linalg.eig(L1)
print("eigenvalue of matrix L1:\n{}".format(eva1))
max_abs_value1 = max(map(abs, eva1))
print("lamda1中绝对值的最大值：", max_abs_value1)
if max_abs_value1>1:
print("矩阵的谱半径大于1,该矩阵发散,无法求解x的值")
if max_abs_value1<1:
print("矩阵的谱半径小于1,该矩阵收敛")
得出的结果如下:
此时谱半径小于1,可以进行迭代求解x的值:
以下是求解的迭代的过程:设定eps=1e-6
x_test1=np.zeros(n)
x_result=G1+np.dot(L1,x_test1)
k=0
print(x_result)
eps=1e-6
while max(map(abs, x_test1-x_result))>eps:
x_test1=x_result
x_result=G1+np.dot(L1,x_test1)
k=k+1
print(x_test1)
print(k)
结果显示再经过14521次迭代之后得到这样的结果：
故此时x的值为：
首先判断松弛因子是否存在:
#SOR迭代法
rob=max_abs_value#max_abs_value=p(B) 谱半径
flag=1-rob*rob
if flag<0:
print("松弛因子不存在,取w=1.5!")
w=1.5
if flag>=0:
w=2/(1+cmath.sqrt(flag))
由于
,所以flag<0,松弛因子不存在，故我们自行设定松弛因子为1.5
再求Lw矩阵的值:
Tw=np.linalg.inv(D-w*L)
Gw=(1-w)*D+w*U
Lw=np.dot(Tw,Gw)
print(Lw)
再求Lw的谱半径的值：
(eva2, evt2) = np.linalg.eig(Lw)
print("eigenvalue of matrix B:\n{}".format(eva2))
max_abs_value2 = max(map(abs, eva2))
print("lamda中绝对值的最大值：", max_abs_value2)
lamda2=max_abs_value2
if lamda2>1:
print("矩阵的谱半径大于1,该矩阵发散,无法求解x的值")
if lamda2<=1:
print("SOR迭代法收敛!")
结果如下:
那么迭代法收敛就进行迭代：
x_test2=np.zeros(n)
x_result2=np.dot(Lw,x_test2)+w*np.dot(Tw,b)
k1=0
eps1=1e-5
print(x_result2)
while max(map(abs, x_test2-x_result2))>eps1:
x_test2=x_result2
x_result2=np.dot(Lw,x_test2)+w*np.dot(Tw,b)
k1=k1+1
print(x_test2)
print(k1)
设定eps1=1e-6,求出的结果如下所示:
经过15533次的迭代之后，得到x的值为:
在n=6的时候用高斯消去法直接求解出了x的精确解，而利用SOR迭代法或者G-S迭代法则是通过不断相当大次数的迭代来求解出x的近似解，在eps=1e-6的情况下还是比较精确的，然而在这道题中我们没办法利用Jacobi迭代法来进行x的求解。
---------------------------ENDQ1-------------------------------
Q2:逐步增大问题的维数，仍用上述的方法来解它们，计算的结果如何? 计算的结果说明的什么?
由于上述证明了Jacobi迭代法不能进行迭代求解，以下就展示剩下三个方法的结果。
当n=20时：
高斯消去法绘制的结果图:
我们可以发现n=20的时候，x的值大幅度溢出了，这就是病态方程的特征，其在计算过程中存在一些特别小的数，使得经过相除使得得到的结果远远超出了我们预计，比如在后面的值都超过了10000，最高的达到了140000.
相比之下G-S迭代法的效果表现的更加优秀，每一个值都没有溢出，而且都更加接近于精确值，此时迭代了26121次
相比之下SOR迭代法的效果表现的比较优秀，每一个值都没有溢出，而且都更加接近于精确值，此时迭代了34712次
当n=100时:
高斯消去法绘制的效果图如下:
我们可以看到受到随着n的变大，其值的溢出也越来越明显，溢出的指数都达到了e24
所以在规模大的方程组中，高斯消去法的稳定性还是不够的，受到极小数的影响比较明显，溢出程度比较明显。
相比之下G-S迭代法的效果表现的更加优秀，每一个值都没有溢出，而且都更加接近于精确值，此时迭代了30993次
相比之下SOR迭代法的效果表现的比较优秀，每一个值都没有溢出，而且都更加接近于精确值，此时迭代了98084次
但是此时迭代次数已经时G-S迭代法的3倍了，消耗还是比较大的。
---------------------------ENDQ2-------------------------------
(3)讨论病态问题求解的算法
显然，对于病态方程来说，高斯消去法一般是不推荐使用的，或者在n<10的范围内可以适当使用，解决病态问题还是要考虑迭代法来求解，基于本题，我们来讨论G-S迭代法和SOR法求解病态方程的算法都比较合理。但是还是比较推荐G-S迭代法来求解病态方程,两种算法对于X的值的逼近是差不多的，但是G-S迭代法的运算效率要高于SOR方法，特别是在n的规模比较大的时候，两者的运算的差距就显得比较大了，所G-S迭代法还是比较合适的求解病态方程的算法。
当然求解病态方程还可以有以下的细节处理：
1.正则化（Regularization）：正则化是一种常见的处理病态问题的方法。它通过在目标函数中引入额外的约束或惩罚项，来限制问题的解空间，使得问题更加稳定和可靠。例如，岭回归和Lasso回归是常见的线性回归问题的正则化方法。
2.截断（Truncation）和截取（Clipping）：对于一些问题，可能存在极大或极小的输入数据值，导致问题的不稳定性。截断和截取是一种方法，通过限制输入数据的范围来缓解这个问题。例如，在信号处理中，可以对信号进行截断或截取操作，以移除极端的数据值。
3.数据平滑和滤波（Data Smoothing and Filtering）：数据平滑和滤波技术可以用来减少噪声和异常值对病态问题的影响。这些技术包括移动平均、卡尔曼滤波、高斯滤波等，可以通过对数据进行平滑处理来减少噪声的影响。
4.参数调整和模型选择：在病态问题中，模型参数的选择和调整非常重要。不恰当的参数选择可能导致问题的不稳定性。通过调整参数或选择合适的模型，可以改善问题的求解稳定性。
5.数值稳定性技术：在数值计算中，病态问题经常涉及到数值不稳定性，例如数值溢出、舍入误差等。数值稳定性技术，如数值精度控制、矩阵条件数估计、迭代算法的收敛性分析等，可以帮助处理这些数值稳定性问题。
总之解决病态问题的算法还是比较多的,当我们通过合理的数据处理和选择合适的模型等等来降低病态方程极小数据对于求解的影响是比较好的，选择比较合适的方法和技术可以大大降低误差和求解的难度。
---------------------------ENDQ3-------------------------------
最后附上整体的源程序代码：
import numpy as np
from scipy.linalg import solve, norm
import cmath
from scipy.linalg import lu
import matplotlib.pyplot as plt
def draw_bar_chart(x):
# 创建画布和子图
fig, ax = plt.subplots()
# 绘制条形图
ax.bar(range(1, len(x)+1), abs(x))
# 添加标签和标题
ax.set_xlabel('X')
ax.set_ylabel('Y')
ax.set_title('n=100 Gauss-Seidel method')
# 显示图形
plt.show()
# 设置问题的维数为 6
n = 6
# 构建 Hilbert 矩阵 H
H = np.zeros((n, n))
for i in range(1, n + 1):
for j in range(1, n + 1):
H[i-1, j-1] = 1 / (i + j - 1)
print(H)
# 构建解向量 x
x = np.ones(n)
# 计算右端向量 b
b = np.dot(H, x)
#Guess消去法
#高斯消元法
def gaussian_elimination(A, b):
n = len(A)
# 构建增广矩阵
aug_matrix = np.concatenate((A, b.reshape(-1, 1)), axis=1)
# 消元过程
for i in range(n):
pivot = aug_matrix[i, i] # 主元
# 消除下方的元素
for j in range(i + 1, n):
factor = aug_matrix[j, i]/pivot
aug_matrix[j, :] -= factor * aug_matrix[i, :]
return aug_matrix
# 调用高斯消去法求解方程组
aug_matrix = gaussian_elimination(H, b)
# 输出最后的 A 矩阵
UA= aug_matrix[:, :-1]
print("最后的 A 矩阵:")
for i in range(n):
for j in range(n):
if(i>j):
UA[i][j]=0
print(UA)
rUA = np.linalg.inv(UA)
LA=np.dot(H,rUA)
print("LA")
for i in range(n):
for j in range(n):
if(i<j):
LA[i][j]=0
print(LA)
rLA = np.linalg.inv(LA)
y=np.dot(rLA,b)
rUA = np.linalg.inv(UA)
x=np.dot(rUA,y)
print(x)
#Jacobi迭代法
D = np.diag(np.diag(H))
D_inv = np.linalg.inv(D)
L = -np.tril(H, k=-1)
U = -np.triu(H, k=1)
B=np.dot(D_inv,(L+U))
print(D)
print(L)
print(U)
print(B)
(eva, evt) = np.linalg.eig(B)
print("eigenvalue of matrix B:\n{}".format(eva))
max_abs_value = max(map(abs, eva))
print("lamda中绝对值的最大值：", max_abs_value)
if max_abs_value>1:
print("矩阵的谱半径大于1,该矩阵发散,无法求解x的值")
#G-S迭代法
C = np.linalg.inv(D-L)
L1=np.dot(C,U)
G1=np.dot(C,b)
print(L1)
print(G1)
(eva1, evt1) = np.linalg.eig(L1)
print("eigenvalue of matrix L1:\n{}".format(eva1))
max_abs_value1 = max(map(abs, eva1))
print("lamda1中绝对值的最大值：", max_abs_value1)
if max_abs_value1>1:
print("矩阵的谱半径大于1,该矩阵发散,无法求解x的值")
if max_abs_value1<1:
print("矩阵的谱半径小于1,该矩阵收敛")
x_test1=np.zeros(n)
x_result=G1+np.dot(L1,x_test1)
k=0
print(x_result)
eps=1e-6
while max(map(abs, x_test1-x_result))>eps:
x_test1=x_result
x_result=G1+np.dot(L1,x_test1)
k=k+1
print(x_test1)
print(k)
draw_bar_chart(x_test1)
#SOR迭代法
rob=max_abs_value#max_abs_value=p(B) 谱半径
flag=1-rob*rob
if flag<0:
print("松弛因子不存在,取w=1.5!")
w=1.5
if flag>=0:
w=2/(1+cmath.sqrt(flag))
print(w)
Tw=np.linalg.inv(D-w*L)
Gw=(1-w)*D+w*U
Lw=np.dot(Tw,Gw)
print(Lw)
(eva2, evt2) = np.linalg.eig(Lw)
print("eigenvalue of matrix B:\n{}".format(eva2))
max_abs_value2 = max(map(abs, eva2))
print("lamda中绝对值的最大值：", max_abs_value2)
lamda2=max_abs_value2
if lamda2>1:
print("矩阵的谱半径大于1,该矩阵发散,无法求解x的值")
if lamda2<=1:
print("SOR迭代法收敛!")
x_test2=np.zeros(n)
x_result2=np.dot(Lw,x_test2)+w*np.dot(Tw,b)
k1=0
eps1=1e-6
print(x_result2)
while max(map(abs, x_test2-x_result2))>eps1:
x_test2=x_result2
x_result2=np.dot(Lw,x_test2)+w*np.dot(Tw,b)
k1=k1+1
print(x_test2)
print(k1)