import numpy as np
from copy import deepcopy
import matplotlib.pyplot as plt
A = np.array([[2,1],[-3,-5]], dtype=float)
print(A)
[[ 2. 1.] [-3. -5.]]
b = np.array([1,-1], dtype=float)
np.linalg.solve(A,b)
array([ 0.57142857, -0.14285714])
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.figure()
plt.plot([-1,1],[y0_left, y0_right]);
plt.plot([-1,1],[y1_left, y1_right]);
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
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 ]
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]]
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');
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]]
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');