使用python求Jacobi迭代与SOR迭代求解希尔伯特矩阵

希尔伯特矩阵是一个非常特殊的矩阵，由于其行列式的值非常小，导致其求解困难。在这里，我们将使用Jacobi迭代和SOR迭代来求解希尔伯特矩阵。

首先，我们需要生成一个 $n \times n$ 的希尔伯特矩阵 $H_n$，其定义如下：

$$
H_{i,j} = \frac{1}{i+j-1}, \ \ \ \ 1 \leq i,j \leq n
$$

下面是Python代码：

import numpy as np

def hilbert(n):
    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)
    return H

接下来，我们将使用Jacobi迭代和SOR迭代来求解希尔伯特矩阵。对于线性方程组 $Ax=b$，Jacobi迭代的迭代公式为：

$$
x^{(k+1)}_i = \frac{1}{a_{ii}} \left(b_i - \sum_{j \neq i} a_{ij} x^{(k)}_j\right), \ \ \ \ 1 \leq i \leq n
$$

SOR迭代的迭代公式为：

$$
x^{(k+1)}_i = (1-\omega) x^{(k)}_i + \frac{\omega}{a_{ii}} \left(b_i - \sum_{j<i} a_{ij} x^{(k+1)}_j - \sum_{j>i} a_{ij} x^{(k)}_j\right), \ \ \ \ 1 \leq i \leq n
$$

其中，$\omega$ 是松弛因子，一般取值为 $(0,2)$ 之间。

下面是Python代码：

def jacobi(A, b, x0, tol=1e-6, maxiter=1000):
    n = len(b)
    x = x0.copy()
    for k in range(maxiter):
        x_new = np.zeros(n)
        for i in range(n):
            x_new[i] = (b[i] - np.dot(A[i,:], x) + A[i,i]*x[i]) / A[i,i]
        if np.linalg.norm(x_new - x) < tol:
            return x_new
        x = x_new
    raise ValueError("Jacobi method did not converge")
    
def sor(A, b, x0, omega, tol=1e-6, maxiter=1000):
    n = len(b)
    x = x0.copy()
    for k in range(maxiter):
        x_new = np.zeros(n)
        for i in range(n):
            x_new[i] = (1-omega)*x[i] + omega*(b[i] - np.dot(A[i,:i], x_new[:i]) - np.dot(A[i,i+1:], x[i+1:])) / A[i,i]
        if np.linalg.norm(x_new - x) < tol:
            return x_new
        x = x_new
    raise ValueError("SOR method did not converge")

现在我们可以使用这些函数来求解希尔伯特矩阵了。下面是Python代码：

# 生成希尔伯特矩阵
n = 10
H = hilbert(n)

# 生成随机向量
b = np.random.rand(n)

# 初始解
x0 = np.zeros(n)

# Jacobi迭代
x_jacobi = jacobi(H, b, x0)

# SOR迭代
omega = 1.5
x_sor = sor(H, b, x0, omega)

# 输出结果
print("Jacobi method:")
print(x_jacobi)
print("SOR method:")
print(x_sor)

这里我们将 $n$ 取为 10，然后使用随机向量作为右端向量 $b$，并将初始解 $x_0$ 设为全零向量。然后使用Jacobi迭代和SOR迭代分别求解，最后输出结果。你可以根据需要修改这些参数来进行求解。