Notebook

Simple iteration to solve a linear system¶

We want to solve the equation $Ax = b$ . Recall that every root-finding problem has equivalent, associated fixed-point iterations. If $x$ is a solution of $Ax = b$ then $x$ satisfies

$x = (I - A)x + b.$

and we want to find a fixed-point of

$g(x) = (I - A)x + b.$

In [5]:

import numpy as np
A = np.array([[1.1, .1, .3, .4], 
              [.2, .9, -.1, .33], 
              [-.3, .4, 1., .8], 
              [.1, .2, .3, 1]])
b = np.array([[1], 
              [1], 
              [1], 
              [1]])
T = np.eye(4) - A
x = np.array([[0],
              [0],
              [0],
              [0]])
for i in range(40):
    x = np.dot(T,x) + b
np.dot(A,x) - b

Out[5]:

array([[ -5.30764743e-09],
       [ -2.58408872e-09],
       [ -3.60924202e-09],
       [ -3.46897711e-09]])