算法的数值稳定性分析实验代码

数值稳定性是算法在处理数据时保持数值解接近真实解的能力。一个数值稳定的算法可以减小舍入误差的影响，使得计算结果更加准确可靠。在进行数值稳定性分析实验时，通常需要对比不同的算法在相同或相似条件下得到的结果，以及这些结果随输入数据的变化情况。

以下是一个简单的数值稳定性分析实验的代码示例，我们以两种不同的方法来计算线性方程组的解：一个是直接法（如高斯消元法），另一个是迭代法（如雅可比迭代法）。我们将分析两种方法在面对不同条件数的系数矩阵时的数值稳定性。

import numpy as np
from scipy.linalg import solve

# 定义计算条件数的函数
def condition_number(A):
    return np.linalg.cond(A)

# 定义高斯消元法求解线性方程组的函数
def gaussian_elimination(A, b):
    return np.linalg.solve(A, b)

# 定义雅可比迭代法求解线性方程组的函数
def jacobi_iteration(A, b, tol=1e-10, max_iter=1000):
    n = len(b)
    x = np.zeros_like(b, dtype=np.double)
    D = np.diag(A)
    R = A - np.diagflat(D)
    for i in range(max_iter):
        x_new = (b - np.dot(R, x)) / D
        if np.linalg.norm(x_new - x, ord=np.inf) < tol:
            break
        x = x_new
    return x

# 生成随机的系数矩阵和常数项
np.random.seed(0)
n = 10  # 方程组的大小
A = np.random.rand(n, n)
b = np.random.rand(n)

# 计算高斯消元法的解和条件数
x_ge = gaussian_elimination(A, b)
cond_ge = condition_number(A)

# 计算雅可比迭代法的解和条件数
x_jacobi = jacobi_iteration(A, b)
cond_jacobi = condition_number(A)

# 打印结果
print(f"高斯消元法解: {x_ge}")
print(f"雅可比迭代法解: {x_jacobi}")
print(f"高斯消元法条件数: {cond_ge}")
print(f"雅可比迭代法条件数: {cond_jacobi}")

# 分析数值稳定性
if np.linalg.norm(x_ge - x_jacobi) < tol:
    print("两种方法的解非常接近，表明高斯消元法数值稳定性较好。")
else:
    print("两种方法的解相差较大，表明雅可比迭代法在本例中数值稳定性较差。")

在这段代码中，我们首先定义了计算矩阵条件数的函数、高斯消元法求解函数和雅可比迭代法求解函数。然后，我们生成了一个随机的系数矩阵`A`和常数项`b`，分别使用高斯消元法和雅可比迭代法求解线性方程组，并计算了它们的解和条件数。最后，我们通过比较两种方法得到的解之间的差异，分析了两种方法的数值稳定性。