我一直在尝试学习一些数值线性代数，我决定尝试对一维泊松问题实现加权雅可比方法我们使用有限差分来近似问题：

$$-u''(x)=f(x),\qquad u(0)=a,\ u(1)=b,$$

我们将域离散化为其中和在内部，我们近似$x_j=jh,$$0\leq j\leq N$$h=1/N.$

$$\frac{2v[j]-v[j-1]-v[j+1]}{h^2} = f[j],\qquad 1\leq j\leq N-1,$$ v这里，并且这可以设置为矩阵问题，我们应用上述迭代方法。特别是，它的形式为 $v[0]=a,\ v[N] = b.$$v[j]\approx u(x_j),$$f[j]=f(x_j).$ $$v_{new}[j] = (1-\omega)v_{old}[j]+\frac{\omega}{2}\left(v_{old}[j-1]+v_{old}[j+1]\right),\qquad 1\leq j\leq N-1.$$

这就是问题所在：最佳 Jacobi 权重是并且该方法应该（通常）在其中A是拉普拉斯矩阵，D=\text{diag}(A)（右边略高于 1）。当\omega超过2时，我的代码似乎中断了。我一直无法找到任何缺失的因素，并且我已经验证这不是缓慢发散的问题。为了让它打破，我做了一个非常糟糕的初步猜测。为简单起见，我还设置了 a=b=0（存储在代码中的g中）。$\omega=2/3,$

$$\omega<\frac{2}{\lambda_{\max}(D^{-1}A)},$$ $A$$D=\text{diag}(A)$$\omega$$2$$a=b=0$$g$

代码（我注释掉了我对 Jacobi 函数使用列表理解的尝试，因为由于某种原因，错误要高几个数量级 - 如果有人能告诉我我哪里出错了，我也非常感谢！）：

import numpy as np
from scipy.sparse import diags

def Jacobi(x, f, omega, N):
    
    dx2 = 1/(N**2)
    x_new = x[:]
    for i in range(1,N):
        x_new[i] = (1-omega)*x[i] \
                        + 0.5*omega*(x[i-1] + x[i+1] + dx2*f[i])
    # x_new[1:N] = [omega_min1*x[i] \
    #                + omega_div2*(x[i-1] + x[i+1] + dx2*f[i])\
    #                    for i in range(1,N)] 
    return x_new                            

# test Jacobi
g = [0,0] # boundary
N = 2**6  # dimension 

x = 100*np.random.rand(N+1) # initial guess, make it bad
x[0]    = g[0] # fill in BBC's
x[N]    = g[1]

f = np.random.rand(N+1) # forcing; for indexing convenience, elongate f
f[0]    = 0 
f[N]    = 0

omega = 1.5 #relaxation parameter 

# Set up matrix problem for direct solve
band        = [-1*np.ones(N-2), 2*np.ones(N-1), -1*np.ones(N-2)]
offset      = [-1,0,1]
A           = diags(band, offset).toarray()
ff          = f[1:N]
x_true      = np.zeros(N+1)
x_true[0]   = g[0]
x_true[N]   = g[1]
gg          = np.zeros(N-1)
gg[0]       = g[0]
gg[N-2]     = g[1]
RHS         = np.add(ff, N*N*gg)
x_true[1:N] = np.linalg.solve(N*N*A, RHS) # do direct solve 


for i in range(0,2500): # Jacobi loop
    y = Jacobi(x, f, omega, N)
    x  = y[:]

print(np.linalg.norm(np.subtract(x,x_true))) # error 


# find relaxation threshold 
D = np.zeros([N-1,N-1])
for i in range(0,N-1):
    D[i,i]  = 2
D_inv = np.linalg.inv(D)
prod = D_inv.dot(A)
[eigs, evec] = np.linalg.eig(prod)
print(2/max(eigs)) # should not work above this omega

输出

4.778131193676897e-06
1.0006026348474466

编辑：我想这也可能是我对直接求解器的利用（即“真实”解决方案存在错误），但我检查了我生成的矩阵和向量是否正确。此外，这会让人很奇怪，好的$\omega$ s 正在收敛到同一件事......