# Weighted Jacobi Not Working on 1D Poisson (Issue with Optimal $\omega$)

I've been trying to learn some numerical linear algebra, and I decided to try to implement the weighted Jacobi method to the 1D Poisson problem $$-u''(x)=f(x),\qquad u(0)=a,\ u(1)=b,$$ where we approximate the problem using finite differences:

we discretize the domain into $$x_j=jh,$$ where $$0\leq j\leq N$$ and $$h=1/N.$$ On the interior, we approximate

$$\frac{2v[j]-v[j-1]-v[j+1]}{h^2} = f[j],\qquad 1\leq j\leq N-1,$$ with $$v[0]=a,\ v[N] = b.$$ Here, $$v[j]\approx u(x_j),$$ and $$f[j]=f(x_j).$$ This can be set up as a matrix problem, to which we apply the aforementioned iterative method. In particular, it takes the form $$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.$$

Here's the issue: the optimal Jacobi weight is $$\omega=2/3,$$ and the method should only (generically) converge when $$\omega<\frac{2}{\lambda_{\max}(D^{-1}A)},$$ where $$A$$ is the Laplacian matrix and $$D=\text{diag}(A)$$ (the right is a little above 1). My code seems to break when $$\omega$$ exceeds $$2$$ instead. I have not been able to find any missing factors, and I have verified that it is not a matter of slow divergence. To try to make it break, I took a really bad initial guess. I also set $$a=b=0$$ for simplicity (stored in $$g$$ in the code).

Code (I commented out my attempt at using list comprehension for my Jacobi function since the error was orders of magnitude higher for some reason - if someone could tell me where I went wrong on that, I'd really appreciate that, too!):

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]
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


Output

4.778131193676897e-06
1.0006026348474466


EDIT: I guess it could also be my utilizaton of the direct solver (i.e. there is an error with the "true" solution), but I checked that my generation of my matrix and vectors are correct. Also, that would make it weird that the good $$\omega$$s are converging to the same thing...

• For me it worked with just increasing the iteration count by a factor of ten. But I was surprised that it worked, usually you choose $\omega \in (0, 1)$. May 3 at 12:49