Math 753/853 HW3: LU decomposition, no pivoting

Problem 1. Write a function named ludecomp that computes the LU decomposition of a square matrix $A$ and returns a lower-triangular matrix $L$ and an upper-triangular matrix $U$ such that $A=LU$.

In [1]:
function ludecomp(A)
    (m,n) = size(A)
    if m  n
        println("error in ludecomp(A): A is not square")
        return (zeros(0,0),zeros(0,0))
    end
    X = copy(A)
    L = zeros(m,m)
    
    for j = 1 : m           # (loop over cols of X, zeroing subdiagonal elements of each)
        L[j,j] = 1
        for i = j+1 : m     # (loop over subdiag rows of jth col of X)
            L[i,j] = X[i,j]/X[j,j]
            for k = j : m   
                X[i,k] -= L[i,j]*X[j,k];  # (update ith row of X past diagonal)
            end
        end
    end
    (L,X)                   # X is the upper-triabgular U at this point
end
Out[1]:
ludecomp (generic function with 1 method)

Problem 2. Test your LU decomposition function on a random 5 x 5 matrix. Use Julia's randn function to generate the random matrix with normally distributed elements.

In [2]:
A = rand(5,5)
Out[2]:
5×5 Array{Float64,2}:
 0.638923  0.275339  0.117635  0.0532384  0.327422
 0.878361  0.134522  0.291026  0.709555   0.488893
 0.641023  0.662     0.74987   0.877505   0.621356
 0.725091  0.398055  0.52983   0.745761   0.179988
 0.112174  0.912097  0.95977   0.166219   0.686541
In [3]:
(L,U) = ludecomp(A);
In [4]:
L
Out[4]:
5×5 Array{Float64,2}:
 1.0        0.0       0.0        0.0     0.0
 1.37475    1.0       0.0        0.0     0.0
 1.00329   -1.58096   1.0        0.0     0.0
 1.13486   -0.350743  0.528155   1.0     0.0
 0.175567  -3.53996   1.67033   11.1497  1.0
In [5]:
U
Out[5]:
5×5 Array{Float64,2}:
  0.638923      0.275339  0.117635   0.0532384   0.327422 
  0.0          -0.244002  0.129307   0.636365    0.0387682
  0.0           0.0       0.836277   1.83016     0.354149 
 -1.11022e-16   0.0       0.0       -0.0580629  -0.365039 
 -1.38778e-17   0.0       0.0        0.0         4.24483  
In [6]:
L*U
Out[6]:
5×5 Array{Float64,2}:
 0.638923  0.275339  0.117635  0.0532384  0.327422
 0.878361  0.134522  0.291026  0.709555   0.488893
 0.641023  0.662     0.74987   0.877505   0.621356
 0.725091  0.398055  0.52983   0.745761   0.179988
 0.112174  0.912097  0.95977   0.166219   0.686541
In [7]:
A-L*U
Out[7]:
5×5 Array{Float64,2}:
 0.0          0.0  0.0          0.0           0.0        
 0.0          0.0  0.0          0.0           0.0        
 0.0          0.0  0.0          0.0           0.0        
 0.0          0.0  0.0          0.0          -5.55112e-17
 1.23512e-15  0.0  1.11022e-16  2.77556e-17   0.0        

Excellent! L is lower-triangular, U is upper-triangular (except for some nearly-zero entries), and A=L*U.

Problem 3. Write a forwardsub function that does forward substitution forward to solve $Ly=b$ for lower-triangular $L$.

In [8]:
function forwardsub(L,b)
    (m,n) = size(L)
    # should check squareness, size of b, and 1's on L diag.
    y = zeros(m)
    for i=1:m
        y[i] = b[i]
        for j=1:i-1
            y[i] -= L[i,j]*y[j]
        end
    end
    y
end
Out[8]:
forwardsub (generic function with 1 method)

Test your forwardsub function on a known $Ly=b$ problem. That is, construct a lower-triangular $L$ and a random $y$ vector and multiply them together to get a vector $b$. Then use forwardsub to solve $Ly=b$ for $y$ given $L$ and $b$, and verify that you get the $y$ you started with.

In [9]:
y = randn(5)
Out[9]:
5-element Array{Float64,1}:
  1.56557 
  0.790595
  0.510563
  0.302395
 -0.442511
In [10]:
b = L*y
Out[10]:
5-element Array{Float64,1}:
 1.56557 
 2.94287 
 0.831387
 2.07147 
 1.2581  
In [11]:
yₛ = forwardsub(L,b)
Out[11]:
5-element Array{Float64,1}:
  1.56557 
  0.790595
  0.510563
  0.302395
 -0.442511
In [12]:
norm(y-yₛ)
Out[12]:
4.710277376051325e-15

Excellent! The forwardsub function solves the lower-triangular system $Ly=b$ to roughly machine precision.

Problem 4. Write a backwardsub function that does backward substitution forward to solve $Ux=y$ for upper-triangular $U$. Test it in a similar manner as your test for forwardsub.

In [13]:
function backwardsub(U,y)
    (m,n) = size(U)
    # should check squareness, size of b, and 1's on L diag.
    x = zeros(m)
    for i = m:-1:1
        x[i] = y[i]
        for j=i+1:m
            x[i] -= U[i,j]*x[j]
        end
        x[i] /= U[i,i]
    end
    x
end
Out[13]:
backwardsub (generic function with 1 method)
In [14]:
x = randn(5)
Out[14]:
5-element Array{Float64,1}:
 -0.230904 
  0.0744537
 -0.870543 
  1.80921  
  0.567712 
In [15]:
y = U*x
Out[15]:
5-element Array{Float64,1}:
  0.0527648
  1.04259  
  2.78418  
 -0.312285 
  2.40984  
In [16]:
xₛ = backwardsub(U,y)
Out[16]:
5-element Array{Float64,1}:
 -0.230904 
  0.0744537
 -0.870543 
  1.80921  
  0.567712 
In [17]:
norm(x-xₛ)
Out[17]:
2.269856336252304e-15

Excellent! The backsub function solves the upper-triangular system to roughly machine precision.

Problem 5. Write an lusolve function that solves $Ax=b$ for $x$ using ludecomp, forwardsub, and backwardsub. Test it on a problem of your choosing.

In [18]:
function lusolve(A,b)
    (L,U) = ludecomp(A)
    y = forwardsub(L,b)
    x = backwardsub(U,y)
end
Out[18]:
lusolve (generic function with 1 method)
In [19]:
A = randn(32,32)
Out[19]:
32×32 Array{Float64,2}:
 -1.14203      1.29426     -0.827464   …  -0.135607    -0.175939 
  1.72502     -1.25533      1.28829        0.345307    -0.428494 
 -0.713377    -0.794327     0.20598        0.2692       1.17688  
 -1.26138      0.265357    -2.03715       -0.954992     0.258151 
  0.504717    -0.687678     0.0894129      1.0473      -0.338715 
 -1.09387      0.108578     0.792859   …   0.411414     0.787691 
 -0.058428     0.887398    -1.31249        1.14785     -0.882837 
 -0.446303    -0.14511      0.195519      -1.67517      1.35324  
  1.53379     -1.79658     -0.401472      -2.02185     -0.399987 
  0.157721    -0.362915    -1.24123        1.42261      0.706034 
  0.764117    -0.740615     0.556777   …  -0.364943    -0.450574 
 -0.982284    -0.46654     -0.973068      -0.25046     -0.791488 
  2.13121      0.30669      0.227247      -1.92988     -0.323416 
  ⋮                                    ⋱   ⋮                     
  0.364539    -1.15344     -0.866293   …   0.858022    -1.06618  
  0.0549567   -0.385022    -0.155702       0.567154     0.601677 
 -0.81295      0.890332     0.0748725      0.37526      1.13853  
 -0.251428    -0.900898     2.36577       -0.700241     0.396922 
 -0.00672333   0.439518     0.112386       2.35113      0.908582 
 -0.631065    -0.281876     0.249062   …   1.5095      -0.630645 
  0.906668    -0.00807386   0.525609      -1.78912      0.0233268
  0.633036     0.0722389   -1.28179       -0.796893    -1.53184  
  0.74268      0.634719    -1.3821        -1.41824      1.53221  
  0.837863     0.618835    -0.611904      -0.37337     -1.26904  
  0.015245    -0.491955     0.914065   …  -0.00297802   0.474122 
 -0.232295    -0.917687     0.101832      -1.5402      -1.37999  
In [20]:
x = randn(32)
Out[20]:
32-element Array{Float64,1}:
  0.143011 
 -0.937668 
  0.336905 
 -1.35552  
  0.948976 
  1.54148  
  1.01796  
 -0.505017 
  0.0198406
  1.35139  
 -0.0896732
  0.785742 
  0.723138 
  ⋮        
  0.0340658
  0.885109 
  0.591093 
  1.22493  
  0.179583 
  0.377693 
  0.395388 
 -1.21301  
 -0.759428 
 -0.0218676
 -0.588359 
  0.401985 
In [21]:
b = A*x
Out[21]:
32-element Array{Float64,1}:
   7.42468 
  -0.773313
  -8.18681 
   0.514739
  -0.877526
  12.238   
  -1.78862 
  -1.25269 
   8.74012 
 -10.3829  
   0.277544
   2.46259 
   7.70911 
   ⋮       
  -4.39365 
  -2.04003 
  10.3682  
  17.9381  
  -0.425204
   0.354712
   6.73621 
  -1.61831 
  -0.763436
   2.3343  
   0.620416
  -2.28625 
In [22]:
xₛ = lusolve(A,b)
Out[22]:
32-element Array{Float64,1}:
  0.143011 
 -0.937668 
  0.336905 
 -1.35552  
  0.948976 
  1.54148  
  1.01796  
 -0.505017 
  0.0198406
  1.35139  
 -0.0896732
  0.785742 
  0.723138 
  ⋮        
  0.0340658
  0.885109 
  0.591093 
  1.22493  
  0.179583 
  0.377693 
  0.395388 
 -1.21301  
 -0.759428 
 -0.0218676
 -0.588359 
  0.401985 
In [23]:
norm(xₛ - x)
Out[23]:
2.9342094044273965e-12

Excellent! My lusolve function found $x$ to 12 digits accuracy.

In [ ]: