Hw5soln
.pdf
School
University of Colorado, Boulder**We aren't endorsed by this school
Course
APPM 5600
Subject
Mathematics
Date
Dec 18, 2024
Pages
7
Uploaded by PrivateChinchilla4899
Copyright APPM2021APPM 5600 — HOMEWORK # 51. (30 points) LetA∈Rn×nbe a tridiagonal matrix where the diagonal entries are given byajforj= 1, . . . , n, the lower diagonal entries arebjforj= 2, . . . , nand the upper diagonalentries arecjforj= 1, . . . , n-1.(a) (10 points) Forn= 3, derive the LU factorization of the matrixA.(b) (10 points) What is the extension of the LU factorization for generaln?(c) (10 points) Theorem 8.2 in Atkinson states that Gaussian Elimination applied to atridiagonal matrix satisfying certain diagonal-dominance conditions does not requirepivoting. What is the operation count (give an exact formula) when applying GaussianElimination to a tridiagonal system without pivoting?You must explain how youderive the operation count.Soln:(a)L=100b1α1100b2α21U=α1c100α2c200α3whereα1=a1,α2=a2-b1α1c1andα3=a3-b2α2c2.(b)Lis going to only have nonzero entries on the diagonal and the lower diagonal. Thediagonal entries are 1 and the lower diagonal entries areγi=biαifori= 1, . . . , n-1whereα1=a1andαi=ai-γi-1cifori= 2, . . . , n. The matrixUonly have nonzeroentries on the diagonal and the upper diagonal. The diagonal entries areαjforj= 1, . . . , n. The upper diagonal entries arecjforj= 1, . . . , n-1.(c) One division is required to create eachγi. There aren-1 of them.It cost nothing to make the upper diagonal entries ofU.There are two operations in creatingαifori= 2, . . . , n: multipliation and subtraction.There aren-1 of each. Thus the total cost is 3(n-1).1
Copyright APPM20212. (15 points) Consider the linear system6x+ 2y+ 2z=-22x+ 2/3y+ 1/3z= 1x+ 2y-z= 0(a) (3 points) Verify that (x, y, z) = (2.6,-3.8,-5) is the exact solution.(b) (5 points) Using 4 digit floating point arithmetic with rounding, solve the system viaGaussian elimination without pivoting.(c) (5 points) Repart part (a) with partial pivoting.(d) (2 points) Which method is more accurate? i.e. stable.(Remember to do the rounding to 4 significant digits as the machine would.)Soln:Coming soon.2
Copyright APPM20213. (25 points) Consider the systemAx=bwhereA=4-10-100-14-10-100-14-10-1-10-14-100-10-14-100-10-14andb=212212.Use the ones vector asx0. (i.e.x0= [1 1 1 1 1 1]T)(a) (5 points) Use Gauss-Jacobi iteration to approximate the solution to this problem withϵ= 1e-7.(b) (5 points) Use Gauss-Siedel iteration to approximate the solution to this problem withϵ= 1e-7.(c) (5 points) Use SOR iteration withω= 1.6735 to approximate the solution to thisproblem withϵ= 1e-7.(d) (2 points) Which method converges faster? Do you expect this to always be true?(e) (5 points) Setc=ρ(B) (spectral radius). Use the following to error estimate to deriveerror bounds for the last computed approximations with all methods.∥xk+1-x∥ ≤c1-c∥xk+1-xk∥(f) (3 points) What happens if you change the parameterωfor SOR?Soln:Soln:a-c The code is here.functionHW5_prob1A = diag(4*ones(6,1))+diag(-1*ones(5,1),-1)+diag(-1*ones(5,1),1)+...diag(-1*ones(3,1),3)+diag(-1*ones(3,1),-3);L = tril(A,-1);U = triu(A,1);D = diag(4*ones(6,1));b = [2 1 2 2 1 2]’;x0 = ones(6,1);Nmax = 1000;omega =1.2;tol = 1e-7;xext = inv(A)*b;3
Copyright APPM2021[alpha,alpham1,c,iter] = Jacobi(L,U,D,b,x0,tol,Nmax);fprintf(’Jacobi:%i iterations \n’, iter);fprintf( [’Jacobi: Approximate solution: ’, ...repmat(’%d ’, 1, numel(alpha)), ’\n’], alpha);fprintf( ’Jacobi: Error bound:%1.8d\n’, c/(1-c)*norm(alpha-alpham1));[alpha,alpham1,c,iter] = Siedel(L,U,D,b,x0,tol,Nmax);fprintf(’Siedel:%i iterations \n’, iter)fprintf( [’Siedel: Approximate solution: ’, ...repmat(’%d ’, 1, numel(alpha)), ’\n’], alpha);fprintf( ’Siedel: Error bound:%1.8d\n’, c/(1-c)*norm(alpha-alpham1));[alpha,alpham1,c,iter] = SOR(L,U,D,b,omega,x0,tol,Nmax);fprintf(’SOR:%i iterations \n’, iter)fprintf( [’SOR: Approximate solution: ’, ...repmat(’%d ’, 1, numel(alpha)), ’\n’], alpha);fprintf( ’SOR: Error bound:%1.8d\n’, c/(1-c)*norm(alpha-alpham1));returnfunction [alpha,alpham1,c,iter] = Jacobi(L,U,D,b,x0,tol,Nmax)M = D;N = -(L+U);Minv = inv(M);B = Minv*N;c= norm(B,2);for j =1:Nmaxxk = Minv*N*x0+Minv*b;if norm(xk-x0)/norm(xk)<tolalpha = xk;alpham1 = x0;iter = j;returnendx0 = xk;end4
Copyright APPM2021alpha = NaN;iter = j;returnfunction [alpha,alpham1,c,iter] = Siedel(L,U,D,b,x0,tol,Nmax)M = D+L;N = -U;Minv = inv(M);B = Minv*N;c = max(abs(eig(B)));for j =1:Nmaxxk = Minv*N*x0+Minv*b;if norm(xk-x0)/norm(xk)<tolalpha = xk;alpham1 = x0;iter = j;returnendx0 = xk;endalpha = NaN;iter = j;returnfunction [alpha,alpham1,c,iter] = SOR(L,U,D,b,omega,x0,tol,Nmax)M = D+omega*L;N = -omega*U+(1-omega)*D;Minv = inv(M);B = Minv*N;c = max(abs(eig(B)));for j =1:Nmaxxk = Minv*N*x0+omega*Minv*b;5
Copyright APPM2021if norm(xk-x0)/norm(xk)<tolalpha = xk;alpham1 = x0;iter = j;returnendx0 = xk;endalpha = NaN;iter = j;return(d) Jacobi takes 37 iterations, Gauss-Siedel takes 21 iterations and SOR takes 45 iterations.So Gauss-Siedel is the fastest. I’m guessing that the parameter for SOR is not optimal.(e) For Jacobi,c= 0.6830. For Gauss-Siedel,c= 0.6123. For SOR,c=∥B∥= 1.0724. Ifyou use spectral radiusc= 0.7257. The Jacobi error bound is 5.23856206e-07. TheGauss-Siedel error bound is 1.76418450e-07. And the SOR error bound is7.73284731e-07. So we are getting roughly the accuracy that we are asking for fromthe algorithm.(f) The number of iterations to convergences changes. It appears thatω∼1.2 is nearlyoptimal for this problem.c= 0.3124 which is much smaller than before. Also the errorbound is better. It is 1.10178857e-07.6
Copyright APPM20214. (30 points) The linear system of equations[1-a-a1]x=bwhereais a real number can under certain conditions can be solved by the iterative method[10-ωa1]xk+1=[1-ωωa01-ω]xk+ωb.(a) (15 points)For which values ofais the method convergent forω= 1?(b) (15 points) Fora= 0.5, find the value ofω∈ {0.8,0.9,1.0,1.1,1.2,1.3}which minimizesthe spectral radius of the matrix[10-ωa1]-1[1-ωωa01-ω].Soln:(a)B=[10-ωa1]-1[1-ωωa01-ω]Whenω= 1, this gets simplified.B=[10-a1]-1[0a00]=[0a0a2]We what the eigenvalues ofBto be less than 1 which means that we need|a|<1.(b) A table with the values is listed below.ωρ(B)0.80.47600.90.37591.00.25001.10.10001.20.20001.30.3000Thus for this problem the optimal parameter choice isω= 1.1. The iterations willconverge the fastest.7