jacobi.py

#
# Jacobi iterative scheme
# Vincent Legat - 2018
# Ecole Polytechnique de Louvain
#

from numpy import *
from scipy.linalg import norm

def g(x):
  y = copy(x)
  y[0] =  (-3*x[1]+6)/5
  y[1] =  -(6*x[0]-4)/8
  return y

def jacobi(x,tol,nmax):
  n = 0; delta = tol+1;  
  x = array(x,dtype=float)  
  while (norm(delta) > tol and n < nmax):
    n = n+1
    xold = x.copy()
    x = g(x)
    delta = x - xold
    print(" Estimated error %14.7e at iteration %d" % (norm(delta),n))
  return x

print(jacobi([0,0],10e-6,50))