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



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

  • $\begingroup$ 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)$. $\endgroup$
    – Dan Doe
    May 3 at 12:49


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.