In [1]:
# The LU factorization doesn't work.
#=
Here is an example where LU totally fails!
Closely related problems arise naturally
in approximating the solution of a particular
class of boundary value problems.
--- TODO -- work these out more!
https://dl.acm.org/citation.cfm?id=154806
=#
using LinearAlgebra
# define a matrix with a bad growth factor
function bad_growth(n)
    A = Matrix(1.0I,n,n)
    A = A .- tril(ones(n,n),-1)
    A[:,end] .= 1
    return A
end
bad_growth(5)
Out[1]:
5×5 Array{Float64,2}:
  1.0   0.0   0.0   0.0  1.0
 -1.0   1.0   0.0   0.0  1.0
 -1.0  -1.0   1.0   0.0  1.0
 -1.0  -1.0  -1.0   1.0  1.0
 -1.0  -1.0  -1.0  -1.0  1.0
In [2]:
##
using Printf
function example(n)
    xtrue = (2*(rand(n,1) .< 0.5).-1)./10
    A = bad_growth(n)
    b = A*xtrue
    xcomp = A\b
    @printf("\n  n = %i\n\n", n);
    @printf("    [ %25s %25s ]\n", "x true", "x computed");
    for i=1:n
        @printf("    [ %25.17f %25.17f ]\n", xtrue[i], xcomp[i]);
    end
    @printf("norm(A*x - b) = %.18e", norm(A*xcomp - b))
end
example(10)
  n = 10

    [                    x true                x computed ]
    [      -0.10000000000000001      -0.10000000000000001 ]
    [      -0.10000000000000001      -0.10000000000000003 ]
    [      -0.10000000000000001      -0.09999999999999998 ]
    [       0.10000000000000001       0.10000000000000009 ]
    [      -0.10000000000000001      -0.09999999999999987 ]
    [      -0.10000000000000001      -0.09999999999999964 ]
    [       0.10000000000000001       0.10000000000000053 ]
    [      -0.10000000000000001      -0.09999999999999787 ]
    [       0.10000000000000001       0.10000000000000142 ]
    [      -0.10000000000000001      -0.10000000000000001 ]
norm(A*x - b) = 5.092070308555669710e-15
In [3]:
##
example(20)
  n = 20

    [                    x true                x computed ]
    [      -0.10000000000000001      -0.10000000000000001 ]
    [       0.10000000000000001       0.10000000000000001 ]
    [       0.10000000000000001       0.09999999999999998 ]
    [      -0.10000000000000001      -0.09999999999999998 ]
    [       0.10000000000000001       0.10000000000000009 ]
    [      -0.10000000000000001      -0.09999999999999964 ]
    [       0.10000000000000001       0.10000000000000053 ]
    [       0.10000000000000001       0.10000000000000142 ]
    [       0.10000000000000001       0.10000000000000142 ]
    [      -0.10000000000000001      -0.09999999999999432 ]
    [      -0.10000000000000001      -0.09999999999999432 ]
    [      -0.10000000000000001      -0.09999999999996589 ]
    [       0.10000000000000001       0.10000000000002274 ]
    [      -0.10000000000000001      -0.09999999999990905 ]
    [       0.10000000000000001       0.10000000000013642 ]
    [      -0.10000000000000001      -0.09999999999945430 ]
    [      -0.10000000000000001      -0.09999999999945430 ]
    [      -0.10000000000000001      -0.09999999999854481 ]
    [      -0.10000000000000001      -0.09999999999854481 ]
    [      -0.10000000000000001      -0.10000000000000001 ]
norm(A*x - b) = 4.537761004057923937e-12
In [4]:
##
example(60)
  n = 60

    [                    x true                x computed ]
    [       0.10000000000000001       0.10000000000000001 ]
    [       0.10000000000000001       0.10000000000000003 ]
    [      -0.10000000000000001      -0.09999999999999998 ]
    [      -0.10000000000000001      -0.09999999999999998 ]
    [      -0.10000000000000001      -0.10000000000000009 ]
    [       0.10000000000000001       0.10000000000000009 ]
    [      -0.10000000000000001      -0.09999999999999964 ]
    [      -0.10000000000000001      -0.09999999999999964 ]
    [      -0.10000000000000001      -0.10000000000000142 ]
    [       0.10000000000000001       0.10000000000000142 ]
    [       0.10000000000000001       0.09999999999999432 ]
    [       0.10000000000000001       0.09999999999999432 ]
    [       0.10000000000000001       0.10000000000002274 ]
    [       0.10000000000000001       0.10000000000002274 ]
    [       0.10000000000000001       0.09999999999990905 ]
    [       0.10000000000000001       0.09999999999990905 ]
    [      -0.10000000000000001      -0.10000000000036380 ]
    [       0.10000000000000001       0.09999999999854481 ]
    [      -0.10000000000000001      -0.10000000000218279 ]
    [      -0.10000000000000001      -0.10000000000582077 ]
    [      -0.10000000000000001      -0.10000000000582077 ]
    [       0.10000000000000001       0.09999999997671694 ]
    [      -0.10000000000000001      -0.10000000003492460 ]
    [       0.10000000000000001       0.09999999986030161 ]
    [      -0.10000000000000001      -0.10000000009313226 ]
    [      -0.10000000000000001      -0.10000000055879354 ]
    [       0.10000000000000001       0.09999999962747097 ]
    [       0.10000000000000001       0.09999999776482582 ]
    [       0.10000000000000001       0.09999999776482582 ]
    [       0.10000000000000001       0.09999999403953552 ]
    [       0.10000000000000001       0.09999999403953552 ]
    [      -0.10000000000000001      -0.10000002384185791 ]
    [      -0.10000000000000001      -0.10000002384185791 ]
    [      -0.10000000000000001      -0.10000014305114746 ]
    [      -0.10000000000000001      -0.10000014305114746 ]
    [       0.10000000000000001       0.09999942779541016 ]
    [       0.10000000000000001       0.09999942779541016 ]
    [      -0.10000000000000001      -0.10000228881835938 ]
    [      -0.10000000000000001      -0.10000228881835938 ]
    [      -0.10000000000000001      -0.10000610351562500 ]
    [       0.10000000000000001       0.09999084472656250 ]
    [       0.10000000000000001       0.09997558593750000 ]
    [       0.10000000000000001       0.09997558593750000 ]
    [      -0.10000000000000001      -0.10009765625000000 ]
    [      -0.10000000000000001      -0.10009765625000000 ]
    [       0.10000000000000001       0.09960937500000000 ]
    [       0.10000000000000001       0.09960937500000000 ]
    [       0.10000000000000001       0.09765625000000000 ]
    [      -0.10000000000000001      -0.10156250000000000 ]
    [       0.10000000000000001       0.09375000000000000 ]
    [      -0.10000000000000001      -0.10937500000000000 ]
    [       0.10000000000000001       0.06250000000000000 ]
    [       0.10000000000000001       0.06250000000000000 ]
    [      -0.10000000000000001      -0.25000000000000000 ]
    [       0.10000000000000001       0.00000000000000000 ]
    [       0.10000000000000001       0.00000000000000000 ]
    [       0.10000000000000001       0.00000000000000000 ]
    [       0.10000000000000001       0.00000000000000000 ]
    [      -0.10000000000000001       0.00000000000000000 ]
    [       0.10000000000000001       0.10000000000000001 ]
norm(A*x - b) = 1.121021362485523198e+00
In [5]:
##
n = 60
xtrue = (2*(rand(n,1) .< 0.5).-1)./10
A = bad_growth(n)
b = A*xtrue
Q,R = qr(A)
xcomp = R\(Q'*b)
@show norm(A*xcomp - b)
@show xcomp
norm(A * xcomp - b) = 6.048264093909924e-15
xcomp = [0.09999999999999716; 0.100000000000001; -0.09999999999999952; -0.09999999999999978; 0.10000000000000027; 0.09999999999999942; 0.1000000000000001; -0.0999999999999999; -0.09999999999999994; 0.10000000000000006; 0.10000000000000003; 0.10000000000000028; 0.09999999999999984; 0.10000000000000021; 0.09999999999999985; 0.10000000000000007; 0.1000000000000001; -0.09999999999999999; -0.09999999999999966; -0.1000000000000004; 0.09999999999999935; 0.10000000000000046; 0.09999999999999985; -0.09999999999999992; -0.09999999999999974; -0.0999999999999999; -0.0999999999999998; -0.10000000000000014; 0.09999999999999983; 0.09999999999999974; 0.10000000000000021; 0.09999999999999987; -0.09999999999999981; -0.09999999999999994; -0.09999999999999973; -0.10000000000000017; -0.10000000000000009; -0.09999999999999999; -0.09999999999999994; -0.09999999999999985; 0.09999999999999976; -0.09999999999999996; 0.09999999999999985; 0.10000000000000012; 0.10000000000000012; 0.09999999999999991; 0.09999999999999998; -0.09999999999999999; 0.09999999999999998; -0.09999999999999998; -0.10000000000000002; 0.09999999999999996; 0.09999999999999987; 0.10000000000000009; -0.1; 0.09999999999999995; -0.09999999999999988; -0.09999999999999999; 0.09999999999999996; 0.09999999999999938]
Out[5]:
60×1 Array{Float64,2}:
  0.09999999999999716
  0.100000000000001  
 -0.09999999999999952
 -0.09999999999999978
  0.10000000000000027
  0.09999999999999942
  0.1000000000000001 
 -0.0999999999999999 
 -0.09999999999999994
  0.10000000000000006
  0.10000000000000003
  0.10000000000000028
  0.09999999999999984
  ⋮                  
  0.09999999999999998
 -0.09999999999999998
 -0.10000000000000002
  0.09999999999999996
  0.09999999999999987
  0.10000000000000009
 -0.1                
  0.09999999999999995
 -0.09999999999999988
 -0.09999999999999999
  0.09999999999999996
  0.09999999999999938
In [6]:
## Added in class: We can fix this by addring random noise
A2 = A + 1e-12*randn(60,60)
b = A2*xtrue
A2\b
Out[6]:
60×1 Array{Float64,2}:
  0.10000000000000002
  0.10000000000000016
 -0.09999999999999998
 -0.1000000000000001 
  0.09999999999999996
  0.09999999999999999
  0.1                
 -0.10000000000000005
 -0.09999999999999984
  0.09999999999999981
  0.1000000000000001 
  0.10000000000000013
  0.10000000000000012
  ⋮                  
  0.10000000000000006
 -0.09999999999999995
 -0.09999999999999992
  0.09999999999999999
  0.0999999999999999 
  0.1000000000000001 
 -0.10000000000000017
  0.09999999999999995
 -0.09999999999999978
 -0.10000000000000012
  0.09999999999999998
  0.10000000000000002