In [1]:
## We'll look at a Jacobi type iteration on this system
A = [1 2; -2 1.0]
b = [2.0; 1.0]
A\b
Out[1]:
2-element Array{Float64,1}:
 -0.0
  1.0
In [2]:
##
M = zeros(Int,2)
M[1] = 2 # use row 2 to update x_1
M[2] = 1 # use row 1 to update x_2
function update(x)
  xnew = similar(x)
  for i=1:length(x)
    row = M[i]
    rhs = A[row,:]'*x - A[row,i]*x[i] # there are ways to avoid subtraction :)
    rhs = b[row] - rhs
    xnew[i] = (rhs/A[row,i])
  end
  return xnew
end
update([1,0.5])
Out[2]:
2-element Array{Float64,1}:
 -0.25
  0.5 
In [3]:
##
x = [1.0,0.5]
for iter=1:10
  global x = update(x)
  @show x
end
x = [-0.25, 0.5]
x = [-0.25, 1.125]
x = [0.0625, 1.125]
x = [0.0625, 0.96875]
x = [-0.015625, 0.96875]
x = [-0.015625, 1.0078125]
x = [0.00390625, 1.0078125]
x = [0.00390625, 0.998046875]
x = [-0.0009765625, 0.998046875]
x = [-0.0009765625, 1.00048828125]
In [4]:
## The choice of 1->2 and 2->1 was arbitrary, what about 1-1, 2-2?
M[1] = 1
M[2] = 2
x = [1.0,0.5]
for iter=1:20
  global x = update(x)
  @show x
end
x = [1.0, 3.0]
x = [-4.0, 3.0]
x = [-4.0, -7.0]
x = [16.0, -7.0]
x = [16.0, 33.0]
x = [-64.0, 33.0]
x = [-64.0, -127.0]
x = [256.0, -127.0]
x = [256.0, 513.0]
x = [-1024.0, 513.0]
x = [-1024.0, -2047.0]
x = [4096.0, -2047.0]
x = [4096.0, 8193.0]
x = [-16384.0, 8193.0]
x = [-16384.0, -32767.0]
x = [65536.0, -32767.0]
x = [65536.0, 131073.0]
x = [-262144.0, 131073.0]
x = [-262144.0, -524287.0]
x = [1.048576e6, -524287.0]
In [5]:
##
# I'm in a hurry... and tried this for my homework.
M[1] = 2
M[2] = 1
function update2(x)
  for i=1:length(x)
    row = M[i]
    rhs = A[row,:]'*x - A[row,i]*x[i]
    rhs = b[row] - rhs
    x[i] = rhs/A[row,i]
  end
  return x
end
x = [1.0,0.5]
for i=1:10
  global x = update2(x)
  @show x
end
x = [-0.25, 1.125]
x = [0.0625, 0.96875]
x = [-0.015625, 1.0078125]
x = [0.00390625, 0.998046875]
x = [-0.0009765625, 1.00048828125]
x = [0.000244140625, 0.9998779296875]
x = [-6.103515625e-5, 1.000030517578125]
x = [1.52587890625e-5, 0.9999923706054688]
x = [-3.814697265625e-6, 1.0000019073486328]
x = [9.5367431640625e-7, 0.9999995231628418]
In [6]:
## This is the Gauss Seidel iteration Note that
## we didn't allocate an extra vector of memory.