In [1]:
## Lecture 6, we see the Neumann series method.
In [2]:
## Recall our linear system to find the expected length of that
# simple "random walk" on the integers.
using LinearAlgebra
A = Matrix(SymTridiagonal(ones(11), -0.5*ones(10)))
Out[2]:
11×11 Array{Float64,2}:
  1.0  -0.5   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0
 -0.5   1.0  -0.5   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0
  0.0  -0.5   1.0  -0.5   0.0   0.0   0.0   0.0   0.0   0.0   0.0
  0.0   0.0  -0.5   1.0  -0.5   0.0   0.0   0.0   0.0   0.0   0.0
  0.0   0.0   0.0  -0.5   1.0  -0.5   0.0   0.0   0.0   0.0   0.0
  0.0   0.0   0.0   0.0  -0.5   1.0  -0.5   0.0   0.0   0.0   0.0
  0.0   0.0   0.0   0.0   0.0  -0.5   1.0  -0.5   0.0   0.0   0.0
  0.0   0.0   0.0   0.0   0.0   0.0  -0.5   1.0  -0.5   0.0   0.0
  0.0   0.0   0.0   0.0   0.0   0.0   0.0  -0.5   1.0  -0.5   0.0
  0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0  -0.5   1.0  -0.5
  0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0  -0.5   1.0
In [3]:
##
A[1,2] = 0
A[end,end-1] = 0
b = [0;ones(9);0]
Out[3]:
11-element Array{Float64,1}:
 0.0
 1.0
 1.0
 1.0
 1.0
 1.0
 1.0
 1.0
 1.0
 1.0
 0.0
In [4]:
##
xjulia = A\b
Out[4]:
11-element Array{Float64,1}:
  0.0              
  9.0              
 16.0              
 21.0              
 24.000000000000004
 25.000000000000004
 24.000000000000004
 21.0              
 16.0              
  9.0              
  0.0
In [5]:
##
""" Solve a linear system of equations \$A x = b\$
using a method based
on the Neumann series for the inverse. This will converge
if the spectral radius of \$(I - A)\$ is less than 1.

* `A`:   The square, non-singular matrix to solve
* `b`:   The right hand size
"""
function solve_linsys(A::Matrix,b::Vector; niter=0)
  n = size(A,1)
  @assert(n == size(A,2), "The matrix is not square")
  if niter <= 0
    niter = 10*size(A,1)
  end
  x = copy(b) # make a copy of the right hand side
  r = b - A*x # compute the residual
  iter = 0
  while norm(r) >= n*eps(1.0) && iter <= niter
    x .+= r # this works in place!
    iter += 1
    r = b - A*x
  end
  @show iter
  return x
end
x = solve_linsys(A,b)

iter = 111
Out[5]:
11-element Array{Float64,1}:
  0.0              
  8.97109409048539 
 15.94508684878328 
 20.924323346415022
 23.911148654901982
 24.906458511859263
 23.911148654901982
 20.924323346415022
 15.94508684878328 
  8.97109409048539 
  0.0
In [6]:
##
norm(b-A*x)
Out[6]:
0.0102340503651399
In [7]:
##
x = solve_linsys(A,b; niter=10000)

iter = 692
Out[7]:
11-element Array{Float64,1}:
  0.0              
  8.999999999999993
 15.999999999999988
 20.999999999999982
 23.99999999999998 
 24.99999999999998 
 23.99999999999998 
 20.999999999999982
 15.999999999999988
  8.999999999999993
  0.0
In [8]:
##
norm(b-A*x)
Out[8]:
2.3498992183808826e-15