jacobi method python physis code example

Example: jacobi method in python

import numpy as np
from numpy.linalg import *

def jacobi(A, b, x0, tol, maxiter=200):
    """
    Performs Jacobi iterations to solve the line system of
    equations, Ax=b, starting from an initial guess, ``x0``.

    Terminates when the change in x is less than ``tol``, or
    if ``maxiter`` [default=200] iterations have been exceeded.

    Returns 3 variables:
        1.  x, the estimated solution
        2.  rel_diff, the relative difference between last 2
            iterations for x
        3.  k, the number of iterations used.  If k=maxiter,
            then the required tolerance was not met.
    """
    n = A.shape[0]
    x = x0.copy()
    x_prev = x0.copy()
    k = 0
    rel_diff = tol * 2

    while (rel_diff > tol) and (k < maxiter):

        for i in range(0, n):
            subs = 0.0
            for j in range(0, n):
                if i != j:
                    subs += A[i,j] * x_prev[j]

            x[i] = (b[i] - subs ) / A[i,i]
        k += 1

        rel_diff = norm(x - x_prev) / norm(x)
        print(x, rel_diff)
        x_prev = x.copy()

    return x, rel_diff, k

# Main code starts here
# ---------------------
GL = 1.6
d = 0.8
A = np.array([
    [1.0,      0,      0,      0,   0],
    [ GL, -(d+1),    1.0,      0,   0],
    [  0,      d, -(d+1),    1.0,   0],
    [  0,      0,      d, -(d+1), 1.0],
    [  0,      0,      0,      0, 1.0]])
b = [0.5,      0,      0,      0, 0.1]
x0 = np.zeros(5);

tol = 1E-9
maxiter = 200
x, rel_diff, k = jacobi(A, b, x0, tol, maxiter)
if k == maxiter:
    print(('WARNING: the Jacobi iterations did not '
           'converge within the required tolerance.'))
print(('The solution is %s; within a tolerance of %g, '
        'using %d iterations.' % (x, rel_diff, k)))
print('Solution error = norm(Ax-b) = %g' % \
            norm(np.dot(A,x)-b)) 
print('Condition number of A = %0.5f' % cond(A))
print('Solution from built-in functions = %s' % solve(A, b))