# 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 [ ]: