Download Solving Linear Systems using LU Decomposition: A Lab Report and more Study notes Calculus in PDF only on Docsity!
Lab 2: Solving Linear Systems
Adrian Down
February 11, 2006
Part 1
Discussion
General method
I found a general method for solving systems of equations using LU decom-
position on the fourth homework assignment. I use the notation,
A =
a 1 b 1
c 1
. (^) bn− 1
cn− 1 an
l 1
ln− 1 1
u 1 b 1
.. .
. (^) bn− 1
un
Multiplying the general forms of the L and U matrices and comparing each
element of the product matrix with the original A term by term, I found
u 1 = a 1
uk = ak − lk− 1 bk− 1
lk =
ck
uk
I then performed forwards and backwards elimination, using these forms
of the L and U matrices. Defining the vector UY = Z,
z 1 = f 1
zk = fk − lk− 1 zk− 1
yn =
zn
un
yk =
zk − bkyk+
uk
Comparing with the exact solution
Consider a function such that y
′′ (x) = 0 ∀x ∈ (0, 1) except at one point xk ,
where y
′′ (xk) = −1. Since y
′′ (x) = 0 for all other x, y
′ (x) must be constant
for all other x, and hence y(x) must be linear. At the particular point xk ,
there must be a discontinuous change in y
′ (x) from some initial value α to
some other value −β, where β > 0. Using the difference quotient to evaluate
the derivative of y
′ (x) near xk,
y
′′ (xk) = lim h→ 0
y
′ (xk + h) − y
′ (xk)
h
= lim h→ 0
−(α + β)
h
Hence, α + β must equal h, so that the limit does not diverge for small h.
My solution
In the code, I read in the tri-diagonal matrix of coefficients A as three vectors,
one for each of the sub-, super-, and main diagonal elements of A. I then
computed the components of the LU decomposition using the algorithms
found above, and then used these values to compute first z by forwards
elimination and finally y by backwards elimination, all using for loops.
I found that the solutions that I obtained were all piecewise linear. I
computed the slopes α and β of the two linear components, and I found that
their sum was always equal to h, as expected.
Code
n = 31; % Number of subdivisions of the interval
k = 16; % Locating of y’’(x(k)) = 1
h = 1/(n+1); % Size of each subdivision
Results
n = 31 k = 13 α = 0. 0186 β = 0. 0127 α + β = 0.0312 =
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
0
1
2
3
4
5
6
7
8
x 10
−
n = 31 k = 29 α = 0. 0029 β = 0. 0283 α + β = 0.0312 =
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
0
1
2
3
x 10
−
Part 2
Discussion
Despite the presence of the cosine term, the given differential equation
y
′′ − (cos x)y = − 1
u(1) = a(1); % Perform the LU decoposition from part 1
for j = 1:n-
l(j) = c(j)/u(j);
u(j+1) = a(j+1) - l(j)*b(j);
end
y(1) = f(1);
for j = 2:n
y(j) = f(j) - l(j-1)*y(j-1);
end
x(n) = y(n)/u(n);
for j = n-1:-1:
x(j) = (y(j) - b(j)*x(j+1)) / u(j);
end
close all;
plot([0:h:1],[0 x 0],’b.--’)
title([’Solution to d^{2}y - (cos x)y = -1, n = ’ int2str(n)]);
Results
