310 likes | 540 Views
Lecture 11-II. Iterative approaches for solving linear systems Successive over-relaxation Post-nonlinear system Partial pivoting. Derivation of Successive over-relaxation. Successive over-relaxation method. D=diag(diag(A)); U=triu(A)-D; L=A-triu(A); R=D+w*L; T= inv(R)* ((1-w)*D-w*U);
E N D
Lecture 11-II • Iterative approaches for solving linear systems • Successive over-relaxation • Post-nonlinear system • Partial pivoting Numerical method
Derivation ofSuccessive over-relaxation Numerical method
Successive over-relaxation method D=diag(diag(A)); U=triu(A)-D; L=A-triu(A); R=D+w*L; T=inv(R)*((1-w)*D-w*U); c=w*inv(R)*b; Numerical method
Successive over-relaxation method x it_xor.m A,b A=[10 -1 2 0;-1 11 -1 3; 2 -1 10 -1;0 3 -1 8]; b=[6 25 -11 15]'; inv(A)*b Numerical method
>> it_xor(A,b) It takes 17 iterations to converge ans = 1.00000000300255 2.00000000357295 -1.00000000163021 0.99999999724547 Numerical method
Variantsuccessive over-relaxation Numerical method
Variant Successive over-relaxation D=diag(diag(A)); U=triu(A)-D; L=A-triu(A); R=w*D+D+L; T=inv(R)*(w*D-U); c=inv(R)*b; Numerical method
Variant Successive over-relaxation method x it_vxor.m A,b A=[10 -1 2 0;-1 11 -1 3; 2 -1 10 -1;0 3 -1 8]; b=[6 25 -11 15]'; inv(A)*b Numerical method
>> it_vxor(A,b,w) It takes 17 iterations to converge ans = 1.00000000300255 2.00000000357295 -1.00000000163021 0.99999999724547 Numerical method
w=0 >> it_vxor(A,b,0) It takes 10 iterations to converge ans = 0.99999999998681 1.99999999985957 -0.99999999997639 1.00000000005561 Numerical method
Large w >> it_vxor(A,b,3) It takes 90 iterations to converge ans = 0.99999998371358 1.99999997234472 -0.99999998858309 1.00000003116601 Numerical method
Post-Nonlinear System • Ax=b, A=rand(100,100),b=rand(100,1) • f(Ax)=b, where f is an arbitrary one-dimensional nonlinear function • Given A,b find f and x Numerical method
f(Ax)=b mean(abs(b-f(Ax))) • Multiple solutions mean absolute error: 0.017080 mean square error: 0.000755 f x Numerical method
f(Ax)=b mean absolute error: 0.024637 mean square error: 0.002738 • Multiple solutions f x Numerical method
Matrix inversion • Solve Gx=f • Example Gf.zip >> load Gf.mat; >> inv(G); Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 4.206041e-017. Numerical method
rank ss=cputime; for i=1:1000 rank(G); end ss2=cputime; fprintf('cputime: rank %f \n',ss2-ss); cputime: rank 1.703125 Numerical method
rank ss=cputime; for i=1:50000 rank(G); end ss2=cputime; fprintf('cputime: rank %f \n',ss2-ss); cputime: rank 82.546875 Numerical method
pinv ss=cputime; for i=1:1000 pinv(G); end ss2=cputime; fprintf('cputime: pinv %f \n',ss2-ss); cputime: pinv 7.062500 Numerical method
Inverse of submatrix ss=cputime; for i=1:50000 inv(G(1:61,1:61)); end ss2=cputime; fprintf('cputime: pinv %f \n',ss2-ss); cputime: inv of submatrix 28.937500 Time reduced from 50*7.06 to 28.93 sec Numerical method
Partial pivoting • Maximal column pivoting The pivot or pivot element is the element of a matrix, which is selected first by an algorithm (e.g. Gaussian elimination, Quicksort, Simplex algorithm), to do certain calculations with the matrix Numerical method
A=[10.^-16 59.17;5.29 -6.13]; b=[59.17 46.78]'; B=[A b]; fGauss.m B=fGauss(B) B backward.m ans = 10.0019 1.0000 x=backward(B) 0.0 1.0000 Numerical method
B = 10.^-16 59.1400 59.1700 5.2900 -6.1300 46.7800 Find Exchange row 1 and row 2 temp=B(1,:); B(1,:)=B(2,:); B(2,:)=temp; Numerical method
B = 5.2900 -6.1300 46.7800 10.^-16 59.1400 59.1700 U=fGauss(B); x=backward(U) x = 10.0019 1.0000 Numerical method
A=[30.00 591400;5.29 -6.13]; b=[591700 46.78]'; B=[A b]; Numerical method
Scaled partial pivoting Exchange row i and row p Numerical method
A = 2.1100 -4.2100 0.9210 4.0100 10.2000 -1.1200 1.0900 0.9870 0.8320 >> s=max(abs(A'))' s = 4.2100 10.2000 1.0900 Determine the first pivot by [v p]=max(A(:,1)./s); p = 3 Numerical method
Scaled partial pivoting Function [B,pr]=pss(A,b) n=length(b);B=[A b]; pr=[]; s=max(abs(A'))' Pr=1:n; for i= 1:n • Find the ith pivot row and set it to p • Exchange pr(i) and pr(p) • Exchange row i and row p of matrix B • For each row j > i • Set d to the ratio B(j,i) / B(i,i) • Set B(j,:) to B(j,:)-d*B(i,:) Return B and pr Numerical method
Example A=[2.11 -4.21 .921;4.01 10.2 -1.12;1.09 .987 .832]; b=[2.01 -3.09 4.21]'; spp.m B = 1.0900 0.9870 0.8320 4.2100 0 -6.1206 -0.6896 -6.1396 0 0.0000 -4.9209 -25.1675 pr = 3 1 2 [B,pr]=spp(A,b) Numerical method
Example A=[1.19 2.11 -100 1;14.2 -0.122 12.2 -1;0 100 -99.9 1; 15.3 0.110 -13.1 -1]; b=[1.12 3.44 2.15 4.16]'; x = 0.1768 0.0127 -0.0207 -1.1826 [B,pr]=spp(A,b); x=backward(B); Numerical method
Exercise • Implement the scaled partial pivoting method to improve the forward Gauss elimination method. • Verify your matlab codes by two examples Numerical method
Avoid row swapping • Row swapping results in time consuming • No row swapping • Use a vector v to emulate row swapping • Swap elements in v instead of swapping rows in B • How to revise spp.m to avoid row swapping ? Numerical method