Solving Linear Systems using LU Decomposition: A Lab Report, Study notes of Calculus

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 =      a1 b1 c1 . . . . . . . . . . . . bn−1 cn−1 an      =      1 l1 . . . . . . . . . ln−1 1           u1 b1 . . . . . . . . . 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 u1 = a1 uk = ak − lk−1bk−1 lk = ck uk I then performed forwards and backwards elimination, using these forms 1 of the L and U matrices. Defining the vector UY = Z, z1 = f1 zk = fk − lk−1zk−1 yn = zn un yk = zk − bkyk+1 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 = −1 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 2 n = 31 k = 29 α = 0.0029 β = 0.0283 α + β = 0.0312 = 1 32 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 0.5 1 1.5 2 2.5 3 x 10 −3 Part 2 Discussion Despite the presence of the cosine term, the given differential equation y′′ − (cos x)y = −1 5 is still linear in y. The equation can be easily written in matrix form using the results of part 1 above,            −2 − h2 cos x1 1 1 −2 − h2 cos x2 1 1 −2 − h2 cos x3 1 1 . . . . . . . . . . . . 1 1 −2 − h2 cos xn−1 1 1 −2 − h2 cos xn            ·            y1 y2 y3 ... ... yn−1 yn            = −h2            1 1 1 ... ... 1 1            The coefficient matrix is again tri-diagonal, and the system can be easily solved using the same algorithm developed above. Code n = 100; % Number of subdivisions of the interval h = 1/(n+1); % Size of each subdivision xbin = [0:h:1]; % Subdivide the interval from 0 to 1 a = -2*ones(1, n); % Fill the coefficient matricies for j = 1:n a(j) = a(j) + h^2*cos(xbin(j+1)); end b = ones(1, n-1); c = ones(1, n-1); f = h^2*ones(1, n); l = zeros(1, n-1); u = zeros(1, n); x = zeros(1, n); y = zeros(1, n); 6 u(1) = a(1); % Perform the LU decoposition from part 1 for j = 1:n-1 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: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 7