Ma 449: Numerical Applied Mathematics Model Solutions to Homework Assignment 4 Prof. Wickerhauser Exercise 6 of Section 3.4, p.138. >> X=[1 1 1 ; 1 2 4 ; 1 3 9 ] X = 1. 1. 1. 1. 2. 4. 1. 3. 9. >> Y=[6;5;2] Y = 6. 5. 2. >> ABC=inv(X)*Y = ABC = 5. 2. -1. Exercise 11 of Section 3.4, p.138. >> A=[1 2 0 0 ; 2 3 -1 0 ; 0 4 2 3 ; 0 0 2 -4] A = 1 2 0 0 2 3 -1 0 0 4 2 3 0 0 2 -4 >> b=[7 9 10 12]' b = 7 9 10 12 >> A\b ans = 1 3 2 -2 You may also solve this by forward elimination and back substitution. Algorithm 1 of Section 3.4, p.140. MATLAB Program: % Tridiagonal matrix solver % Indexing: % % Diagonals of the matrix: Right-hand side: % D(1) C(1) B(1) % A(1) D(2) C(2) B(2) % A(2) D(3) C(3) . % ... . % A(n-2) D(n-1) C(n-1) . % A(n-1) D(n) B(n) % function X = trisol(A,D,C,B) n=length(B); % Forward elimination: for i=1:n-1 A(i) = A(i)/D(i); % Overwrite the eliminated coefficient A(i) with L D(i+1) = D(i+1) - A(i)*C(i); % M(i+1,i+1) -= M(i+1,i)*M(i,i+1)/M(i,i) B(i+1) = B(i+1) - A(i)*B(i); % B(i+1) -= M(i+1,i)*B(i)/M(i,i) end % Back substitution: X(n) = B(n)/D(n); % X(n) = B(n)/M(n,n) for i=n-1:-1:1 X(i)=(B(i)-C(i)*X(i+1))/D(i); % X(i)=(B(i)-M(i,i+1)*X(i+1))/M(i,i) end INPUT: subd=[1,1,1]; diag=[4,5,6,7]; supd=[-1,-1,-1]; rhs=[1,0,0,0]; x=trisol(subd, diag, supd, rhs) OUTPUT: x = 0.2385 -0.0462 0.0075 -0.0011 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Exercise 4(a,b) of Section 3.5, p.153. >> A=[4 2 1; 2 5 -2; 1 -2 7]; [L,U,P]=lu(A) L = U = P = 1.0000 0 0 4.0000 2.0000 1.0000 1 0 0 0.5000 1.0000 0 0 4.0000 -2.5000 0 1 0 0.2500 -0.6250 1.0000 0 0 5.1875 0 0 1 >> B=[1 -2 7; 4 2 1; 2 5 -2]; [L,U,P]=lu(B) L = U = P = 1.0000 0 0 4.0000 2.0000 1.0000 0 1 0 0.5000 1.0000 0 0 4.0000 -2.5000 0 0 1 0.2500 -0.6250 1.0000 0 0 5.1875 1 0 0 Exercise 2(a) of Section 3.6, p.165. MATLAB PROGRAM [from the textbook's website, in "jacobi.m"]: function X=jacobi(A,B,P,delta, max1) % Input - A is an N x N nonsingular matrix % - B is an N x 1 matrix % - P is an N x 1 matrix; the initial guess % - delta is the tolerance for P % - max1 is the maximum number of iterations % Output - X is an N x 1 matrix: the jacobi approximation to % the solution of AX = B % NUMERICAL METHODS: Matlab Programs % (c) 2004 by John H. Mathews and Kurtis D. Fink % Complementary Software to accompany the textbook: % NUMERICAL METHODS: Using Matlab, Fourth Edition % ISBN: 0-13-065248-2 % Prentice-Hall Pub. Inc. % One Lake Street % Upper Saddle River, NJ 07458 N = length(B); for k=1:max1 for j=1:N X(j)=(B(j)-A(j,[1:j-1,j+1:N])*P([1:j-1,j+1:N]))/A(j,j); end err=abs(norm(X'-P)); relerr=err/(norm(X)+eps); P=X'; if (err> A=[8 -3; -1 4]; B=[10; 6]; P=[0; 0]; delta=1e-12; max1=3; >> X=jacobi(A,B,P,delta,max1) X = 1.9297 1.9531 >> X=jacobi(A,B,P,delta,1) X = 1.2500 1.5000 >> X=jacobi(A,B,P,delta,2) X = 1.8125 1.8125 >> X=jacobi(A,B,P,delta,3) X = 1.9297 1.9531 >> X=jacobi(A,B,P,delta,4) X = 1.9824 1.9824 This looks like it is converging to (2,2), which indeed solves the system. Exercise 2(b) of Section 3.6, p.165. MATLAB PROGRAM [from the textbook's website, in "gseid.m"]: function X=gseid(A,B,P,delta, max1) % Input - A is an N x N nonsingular matrix % - B is an N x 1 matrix % - P is an N x 1 matrix; the initial guess % - delta is the tolerance for P % - max1 is the maximum number of iterations % Output - X is an N x 1 matrix: the gauss-seidel approximation to % the solution of AX = B % NUMERICAL METHODS: Matlab Programs % (c) 2004 by John H. Mathews and Kurtis D. Fink % Complementary Software to accompany the textbook: % NUMERICAL METHODS: Using Matlab, Fourth Edition % ISBN: 0-13-065248-2 % Prentice-Hall Pub. Inc. % One Lake Street % Upper Saddle River, NJ 07458 N = length(B); for k=1:max1 for j=1:N if j==1 X(1)=(B(1)-A(1,2:N)*P(2:N))/A(1,1); elseif j==N X(N)=(B(N)-A(N,1:N-1)*(X(1:N-1))')/A(N,N); else %X contains the kth approximations and P the (k-1)st X(j)=(B(j)-A(j,1:j-1)*X(1:j-1)'-A(j,j+1:N)*P(j+1:N))/A(j,j); end end err=abs(norm(X'-P)); relerr=err/(norm(X)+eps); P=X'; if (err> A=[8 -3; -1 4]; B=[10; 6]; P=[0; 0]; delta=1e-12; max1=3; >> X=gseid(A,B,P,delta,max1) X = 1.9934 1.9984 >> X=gseid(A,B,P,delta,1) X = 1.2500 1.8125 >> X=gseid(A,B,P,delta,2) X = 1.9297 1.9824 >> X=gseid(A,B,P,delta,3) X = 1.9934 1.9984 >> X=gseid(A,B,P,delta,4) X = 1.9994 1.9998 This also seems to be converging to the solution (2,2). Exercise 6(a,b) of Section 3.7, p.181. MATLAB PROGRAMS: Same as for Exercise 2(a,b) in this section. INPUT/OUTPUT: >> A=[2 8 -1; 5 -1 1; -1 1 4]; B=[11; 10; 3]; P=[0;0;0]; delta=1e-12; max1=3; >> X=jacobi(A,B,P,delta,1) X = 5.5000 -10.0000 0.7500 >> X=jacobi(A,B,P,delta,2) X = 45.8750 18.2500 4.6250 >> X=jacobi(A,B,P,delta,3) X = -65.1875 224.0000 7.6562 >> X=gseid(A,B,P,delta,1) X = 5.5000 17.5000 -2.2500 >> X=gseid(A,B,P,delta,2) X = -65.6250 -340.3750 69.4375 >> X=gseid(A,B,P,delta,3) X = 1.0e+03 * 1.4017 7.0680 -1.4158 Neither method seems to produce a convergent sequence. Exercise 6(a) of Section 3.7, p.182. MATLAB PROGRAM: % Initial value p0 = -0.3; q0 = -1.3; x=p0; y=q0; % Jacobi (repaste to iterate): x1 = (y-x*x*x+3*x*x+3*x)/7; y1 = (y*y+2*y-x-2)/2; x=x1,y=y1 OUTPUT: >> % Initial value p0 = -0.3; q0 = -1.3; x=p0; y=q0; >> % Jacobi (repaste to iterate): x1 = (y-x*x*x+3*x*x+3*x)/7; y1 = (y*y+2*y-x-2)/2; x=x1,y=y1 x = -0.2719 y = -1.3050 >> % Jacobi (repaste to iterate): x1 = (y-x*x*x+3*x*x+3*x)/7; y1 = (y*y+2*y-x-2)/2; x=x1,y=y1 x = -0.2684 y = -1.3176 >> % Jacobi (repaste to iterate): x1 = (y-x*x*x+3*x*x+3*x)/7; y1 = (y*y+2*y-x-2)/2; x=x1,y=y1 x = -0.2696 y = -1.3154 Exercise 6(a,b) of Section 3.7, p.182. MATLAB PROGRAMS: % Initial value p0 = -0.3; q0 = -1.3; x=p0; y=q0; % Gauss-Seidel (repaste to iterate): x = (y-x*x*x+3*x*x+3*x)/7, y = (y*y+2*y-x-2)/2 OUTPUT: >> % Initial value p0 = -0.3; q0 = -1.3; x=p0; y=q0; >> % Gauss-Seidel (repaste to iterate): x = (y-x*x*x+3*x*x+3*x)/7, y = (y*y+2*y-x-2)/2 x = -0.2719 y = -1.3191 >> % Gauss-Seidel (repaste to iterate): x = (y-x*x*x+3*x*x+3*x)/7, y = (y*y+2*y-x-2)/2 x = -0.2704 y = -1.3139 >> % Gauss-Seidel (repaste to iterate): x = (y-x*x*x+3*x*x+3*x)/7, y = (y*y+2*y-x-2)/2 x = -0.2694 y = -1.3160 The intermediate results are different, but both methods seem to be convergent. Algorithm 4(b) of Section 3.7, p.184. MATLAB PROGRAM: format long P=[0; 0; 0] % Iterate by pasting these lines in the Matlab console window: x=P(1); y=P(2); z=P(3); P0 = [x; y; z]; F = [ x*(x-1)+2*y*y+y*z-10; 5*x-6*y+z; z-x*x-y*y ]; dxf1 = 2*x-1; dxf2 = 5; dxf3 = -2*x; dyf1 = 4*y+z; dyf2 = -6; dyf3 = -2*y; dzf1 = y; dzf2 = 1; dzf3 = 1; JF = [ dxf1 dxf2 dxf3 ; dyf1 dyf2 dyf3 ; dzf1 dzf2 dzf3 ]; P = P0 - inv(JF')*F, DP = norm(P-P0), OUTPUT: Starting at point P=[0;0;0], the iteration does not converge. Starting at point P=[0;1;2], the iteration converges as follows: >> format long g >> P=[0;1;2] P = 0 1 2 P = 1.88888888888889 2.11111111111111 3.22222222222222 DP = 2.50924217569694 P = 1.23809131423 1.59918697512171 3.4046652795803 DP = 0.847873382863593 P = 1.11428720305211 1.51353016837175 3.50974499496998 DP = 0.183592736984666 P = 1.11010355911148 1.51098143348823 3.51537080537199 DP = 0.00745980353681243 P = 1.11009927672436 1.51097910954278 3.51537827363488 DP = 8.91709103560239e-06 P = 1.11009927672023 1.51097910954074 3.5153782736433 DP = 9.59606706410858e-12 P = 1.11009927672023 1.51097910954074 3.5153782736433 DP = 4.44089209850063e-16 P = 1.11009927672023 1.51097910954074 3.5153782736433 DP = 4.44089209850063e-16