Jacobi Iteration¶

In [1]:
import numpy as np
from copy import deepcopy
import matplotlib.pyplot as plt

Choose a matrix system to solve¶

In [2]:
A = np.array([[2,1],[-3,-5]], dtype=float)
print(A)

[[ 2.  1.]
 [-3. -5.]]
In [3]:
b = np.array([1,-1], dtype=float)
In [4]:
np.linalg.solve(A,b)
Out[4]:
array([ 0.57142857, -0.14285714])
In [5]:
y0_left = (b[0] - A[0,0]*(-1))/A[0,1]
y0_right = (b[0] - A[0,0]*(1))/A[0,1]
y1_left = (b[1] - A[1,0]*(-1))/A[1,1]
y1_right = (b[1] - A[1,0]*(1))/A[1,1]
In [6]:
plt.figure()
plt.plot([-1,1],[y0_left, y0_right]);
plt.plot([-1,1],[y1_left, y1_right]);

Jacobi Step¶

In [7]:
def Jacobi(A,b,xold):
    '''
     xnew = Jacobi(A, b, xold)
     
     Takes one Jacobi step for a 2x2 system.
    '''
    xnew = np.zeros_like(b)
    xnew[0] = (b[0] - A[0,1]*xold[1]) / A[0,0]
    xnew[1] = (b[1] - A[1,0]*xold[0]) / A[1,1]
    return xnew
In [8]:
x0 = np.array([0.5,0.5])
x1 = Jacobi(A,b,x0)
x2 = Jacobi(A,b,x1)
x3 = Jacobi(A,b,x2)
print(f'{x1}\n{x2}\n{x3}')

[ 0.25 -0.1 ]
[0.55 0.05]
[ 0.475 -0.13 ]

Solve the system iteratively¶

In [9]:
x = np.array([0.5,0.5])
x_hist = [x]
for k in range(10):
    x = Jacobi(A,b,x)
    x_hist.append(x)
x_hist = np.array(x_hist)
print(x_hist)

[[ 0.5       0.5     ]
 [ 0.25     -0.1     ]
 [ 0.55      0.05    ]
 [ 0.475    -0.13    ]
 [ 0.565    -0.085   ]
 [ 0.5425   -0.139   ]
 [ 0.5695   -0.1255  ]
 [ 0.56275  -0.1417  ]
 [ 0.57085  -0.13765 ]
 [ 0.568825 -0.14251 ]
 [ 0.571255 -0.141295]]
In [10]:
plt.figure()
plt.plot([-1,1],[y0_left, y0_right]);
plt.plot([-1,1],[y1_left, y1_right]);
plt.plot(x_hist[:,0], x_hist[:,1], '-o', ms=3, lw=0.5);
plt.plot(x_hist[0,0], x_hist[0,1], 'ko');
plt.title('The black dot shows the first estimate');

NOT diagonally dominant¶

In [11]:
A = np.array([[2,1],[-2,0.9]], dtype=float)
print(A)
x = np.array([0.5,0.5])
x_hist = [x]
for k in range(10):
    x = Jacobi(A,b,x)
    x_hist.append(x)
x_hist = np.array(x_hist)
print(x_hist)

[[ 2.   1. ]
 [-2.   0.9]]
[[ 0.5         0.5       ]
 [ 0.25        0.        ]
 [ 0.5        -0.55555556]
 [ 0.77777778  0.        ]
 [ 0.5         0.61728395]
 [ 0.19135802  0.        ]
 [ 0.5        -0.68587106]
 [ 0.84293553  0.        ]
 [ 0.5         0.76207895]
 [ 0.11896052  0.        ]
 [ 0.5        -0.84675439]]
In [12]:
plt.figure()
y0_left = (b[0] - A[0,0]*(-1))/A[0,1]
y0_right = (b[0] - A[0,0]*(1))/A[0,1]
y1_left = (b[1] - A[1,0]*(-1))/A[1,1]
y1_right = (b[1] - A[1,0]*(1))/A[1,1]
plt.plot([-1,1],[y0_left, y0_right]);
plt.plot([-1,1],[y1_left, y1_right]);
plt.plot(x_hist[:,0], x_hist[:,1], '-o', ms=3, lw=0.5);
plt.plot(x_hist[0,0], x_hist[0,1], 'ko');
In [ ]: