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


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.