Bài giảng Ứng dụng các phương pháp tính số - Chương 3: Hệ phương trình đại số tuyến tính - Ngô Văn Thanh

TS. Ngô Văn Thanh, Viện Vật lý. Cao học vật lý – chuyên ngành Vật lý lý thuyết. Chương.3 Hệ phương trình đại số tuyến tính. 3.1. Phương pháp tính trực tiếp 3.1.1. Phương pháp khử Gaussian. 3.1.2. Quay và tính tỷ lệ - Pivoting and scaling. 3.1.3. Phương pháp phân ly LU. 3.1.4. Nghịch đảo ma trận. 3.1.5. Các hệ chéo bậc 3 3.2. Các phương pháp lặp. 3.2.1. Phương pháp lặp Jacobi. 3.2.2. Phương pháp lặp Gauss-Seidel. @2009, Ngô Văn Thanh - Viện Vật Lý 3.1. Phương pháp tính trực tiếp  Hệ M phương trình tuyến tính với N ẩn.  Biểu diễn dạng ma trận: @2009, Ngô Văn Thanh - Viện Vật Lý 3.1.1. Phương pháp khử Gaussian.  Xét hệ N phương trình tuyến tính với N ẩn số, dạng mở rộng ma trận:  Giải hệ phương trình bằng phương pháp khử theo hàng. Xét ví dụ đơn giản  Nhân hàng 1 cho 3 rồi cộng với hàng thứ 2, thay kết quả cho hàng 2.  Nhân hàng 1 cho -2 rồi cộng với hàng thứ 3, thay kết quả cho hàng 3. @2009, Ngô Văn Thanh - Viện Vật Lý  Hàng 1 vừa rồi gọi là phương trình pivot hay hàng pivot, phần tử a11 gọi là phần tử pivot.  Áp dụng phương pháp tương tự bắt đầu từ phương trình pivot thứ hai và phần tử pivot thứ 2.   Hệ phương trình tương đương bây giờ có dạng  Giải lần lượt từ phương trình thứ 3 đến phương trình 1. @2009, Ngô Văn Thanh - Viện Vật Lý  Trường hợp tổng quát:  Chọn hàng thứ nhất làm hàng pivot, ta có: @2009, Ngô Văn Thanh - Viện Vật Lý  Thực hiện tiếp cho hàng pivot thứ hai:  Cuối cùng ta thu được ma trận chéo, trong đó các phần tử nằm phía dưới đường chéo đều bằng 0. @2009, Ngô Văn Thanh - Viện Vật Lý 3.1.2. Pivoting and scaling.  Nếu các phần tử pivot (các phần tử trên đường chéo) rất bé thì sẽ có sai số do các phép tính làm tròn các con số sau dấu chấm thập phân.  Nếu các phần tử pivot bằng không thì không sử dụng được phương pháp khử Gaussian.  Để giải quyết các vấn đề trên.  Nếu phần tử pivot bằng 0, thực hiện tráo đổi hàng pivot cho hàng kế tiếp.  Sử dụng các phương pháp scaling Phương pháp scaling theo cột (partial pivoting). Phương pháp scaling theo hàng (scaled partial pivoting). Kết hợp cả hai phương pháp trên (full pivoting).  Cuối cùng là tìm các nghiệm bằng phương pháp thay thế ngược (backward substitution) @2009, Ngô Văn Thanh - Viện Vật Lý Phương pháp scaling theo cột (partial pivoting):  Xét phần tử pivot akk, ta có  Tìm giá trị mki lớn nhất, tương đương với việc tìm phần tử có giá trị tuyệt đối lớn nhất trong cột i: |aki|  Tráo đổi hàng pivot cho hàng tương ứng với  Ví dụ:  Phần tử |-3| lớn nhất trong các phần tử của cột 1.  Tráo đổi hàng thứ 2 cho hàng 1.  Tiếp tục thực hiện phương pháp này cho đến khi thu được ma trận rút gọn về dạng chéo. @2009, Ngô Văn Thanh - Viện Vật Lý Phương pháp scaling theo hàng (scaled partial pivoting):  Xác định hệ số tỷ lệ là phần tử lớn nhất trong từng hàng  Tính tỷ số các phần tử trong cột k :  Tìm hàng thứ p mà tỷ số này có giá trị lớn nhất.  Tráo đổi hàng pivot cho hàng thứ p đó.  Xét ví dụ:  Ta có:  Tính:  Vì a31/s3 lớn nhất nên tráo đổi hàng 1 cho hàng 3. @2009, Ngô Văn Thanh - Viện Vật Lý 3.1.3. Phương pháp phân ly LU.  Giả thiết rằng ma trận A có thể viết dưới dạng tích của hai ma trận tam giác:  Suy ra  Ta có hai phương trình ma trận: : giải phương trình này để tìm y : giải phương trình này để tìm x @2009, Ngô Văn Thanh - Viện Vật Lý  Giải phương trình:  Giải phương trình: Phân ly LU để tìm hai ma trận L và U.  Các phần tử của các ma trận thoả mãn biểu thức:  Với i < j :  Với i = j :  Với i > j :  Đây là hệ N 2 phương trình với (N 2 + N ) ẩn số. @2009, Ngô Văn Thanh - Viện Vật Lý  Thuật toán Crout:  Đầu tiên đặt N phần tử trên đường chéo của ma trận L bằng 1.  Chạy vòng lặp j = 1, 2,, N. Bước 1: chạy tiếp vòng lặp trong i = 1,2, , j. Áp dụng điều kiện cho i < j ; i = j và . Tính Bước 2: chạy vòng lặp i = j + 1, j + 2, , N . Áp dụng điều kiện cho i > j. Tính  Kết hợp 2 ma trận L và U: @2009, Ngô Văn Thanh - Viện Vật Lý @2009, Ngô Văn Thanh - Viện Vật Lý  Chương trình đơn giản: SUBROUTINE ludcmp(a,indx,d) IMPLICIT NONE REAL, DIMENSION(:,:), INTENT(INOUT) :: a INTEGER, DIMENSION(:), INTENT(OUT) :: indx REAL, INTENT(OUT) :: d REAL, DIMENSION(size(a,1)) :: vv REAL, PARAMETER :: TINY=1.0e-20 !A small number. INTEGER :: j,n,imax n=assert_eq(size(a,1),size(a,2),size(indx),’ludcmp’) d=1.0 !No row interchanges yet. vv=maxval(abs(a),dim=2) if (any(vv == 0.0)) & write(6,*) ’ There is a row of zeros.’ vv=1.0/vv @2009, Ngô Văn Thanh - Viện Vật Lý do j=1,n imax=(j-1)+imaxloc(vv(j:n)*abs(a(j:n,j))) if (j /= imax) call swap(a(imax,:),a(j,:)) d = -d vv(imax) = vv(j) end if indx(j) = imax if (a(j,j) == 0.0) a(j,j)=TINY a(j+1:n,j) = a(j+1:n,j)/a(j,j) * do i=1,j-1 a(i,j)=a(i,j)-sum(a(i,1:i-1)*a(1:i-1,j)) enddo * do i=j,n !This is i = j and i = j+1: ::N a(i,j)=a(i,j)-sum(a(i,1:j-1)*a(1:j-1,j)) enddo enddo END SUBROUTINE ludcmp @2009, Ngô Văn Thanh - Viện Vật Lý  Chương trình F90 kiểu parallel tối ưu: SUBROUTINE ludcmp(a,indx,d) IMPLICIT NONE REAL, DIMENSION(:,:), INTENT(INOUT) :: a INTEGER, DIMENSION(:), INTENT(OUT) :: indx REAL, INTENT(OUT) :: d REAL, DIMENSION(size(a,1)) :: vv REAL, PARAMETER :: TINY=1.0e-20 !A small number. INTEGER :: j,n,imax n=assert_eq(size(a,1),size(a,2),size(indx),’ludcmp’) d=1.0 !No row interchanges yet. vv=maxval(abs(a),dim=2) if (any(vv == 0.0)) & write(6,*) ’ There is a row of zeros.’ vv=1.0/vv @2009, Ngô Văn Thanh - Viện Vật Lý do j=1,n imax=(j-1) + imaxloc(vv(j:n)*abs(a(j:n,j))) !Find the pivot row (only j:n elements). if (j /= imax) then !Do we need to interchange rows? call swap(a(imax,:),a(j,:)) ! Yes, do so... d = -d !...and change the parity of d. vv(imax)= vv(j) !Also interchange the scale factor. end if indx(j) = imax if (a(j,j) == 0.0) a(j,j) = TINY a(j+1:n,j)=a(j+1:n,j)/a(j,j) !Divide by the pivot element. a(j+1:n,j+1:n)=a(j+1:n,j+1:n) - & outerprod(a(j+1:n,j),a(j,j+1:n)) !Reduce remaining submatrix. end do END SUBROUTINE ludcmp @2009, Ngô Văn Thanh - Viện Vật Lý SUBROUTINE lubksb(a,indx,b) n=assert_eq(size(a,1),size(a,2),size(indx),’lubksb’) k=0 do i=1,n j=indx(i) summ=b(j) b(j)=b(i) if (k /= 0) then summ=summ-dot_product(a(i,k:i-1),b(k:i-1)) else if (summ /= 0.0) then k = i end if !have to do the dot product above. b(i)=summ end do do i=n,1,-1 ! Now we do the backsubstitution. b(i) = (b(i)-dot_product(a(i,i+1:n),b(i+1:n)))/a(i,i) end do END SUBROUTINE lubksb @2009, Ngô Văn Thanh - Viện Vật Lý 3.1.4. Nghịch đảo ma trận. Xét hệ phương trình tuyến tính:  Nhân hai vế cho ma trận nghịch đảo A-1:  Nếu xác định được A-1 thì hệ phương trình hoàn toàn có thể tìm nghiệm dưới dạng:  Để xác định ma trận nghịch đảo A-1, ta phải thực hiện các bước sau:  Áp dụng phương pháp phân ly LU (ludcmp) cho ma trận A  Sử dụng chương trình lubksb để tìm nghiệm cho từng cột của ma trận b @2009, Ngô Văn Thanh - Viện Vật Lý REAL, DIMENSION(:,:), INTENT(INOUT) :: a INTEGER,DIMENSION(size(a,1)) :: indx INTEGER,DIMENSION(size(a,1),size(a,1)) :: y INTEGER :: i,j,n REAL :: d ... n=assert_eq(size(a,1),size(a,2),size(indx),’inverse’) y = 0.0 do i=1,n !Set up identity matrix. y(i,i)=1. enddo call ludcmp(a,indx,d) !Decompose the matrix just once. do j=1,n !Find inverse by columns. call lubksb(a,indx,y(1,j)) enddo @2009, Ngô Văn Thanh - Viện Vật Lý Tính định thức của một ma trận.  Sau khi áp dụng phân ly LU cho ma trận A, tính định thức theo công thức: REAL, FUNCTION det(a) REAL, DIMENSION(:,:), INTENT(INOUT) :: a INTEGER,DIMENSION(size(a,1)) :: indx INTEGER :: i,n REAL :: d n=size(a,1) call ludcmp(a,indx,d) !Decompose the matrix just once. do i=1,n d=d*a(i,i) enddo det = d END FUNCTION det @2009, Ngô Văn Thanh - Viện Vật Lý 3.1.5. Các hệ chéo bậc 3.  Chương trình: SUBROUTINE tridag_ser(a,b,c,r,u) IMPLICIT NONE REAL, DIMENSION(:), INTENT(IN) :: a,b,c,r REAL, DIMENSION(:), INTENT(OUT) :: u REAL, DIMENSION(size(b)) :: gam INTEGER :: n,j REAL :: bet @2009, Ngô Văn Thanh - Viện Vật Lý REAL :: bet n=size(b) bet=b(1) if (bet == 0.0) write(6,*) ‘Error at code stage 1’ u(1) = r(1)/bet do j=2,n !Decomposition and forward substitution. gam(j) = c(j-1)/bet bet = b(j)-a(j-1)*gam(j) if (bet == 0.0) write(6,*) ‘Error at code stage 2’ u(j)=(r(j)-a(j-1)*u(j-1))/bet enddo do j=n-1,1,-1 !Backsubstitution. u(j)=u(j)-gam(j+1)*u(j+1) end do END SUBROUTINE tridag_ser @2009, Ngô Văn Thanh - Viện Vật Lý 3.2. Các phương pháp lặp.  Hệ phương trình tuyến tính dạng có thể giải gầng đúng bằng phương pháp tính lặp. Nghiệm của hệ phương trình có dạng:  Trong đó H được gọi là ma trận lặp.  Ma trận A có thể viết dưới dạng  L : ma trận tam giác dưới hoàn toàn (không có thành phần chéo)  D : ma trận chéo  U : ma trận tam giác trên hoàn toàn. @2009, Ngô Văn Thanh - Viện Vật Lý 3.2.1. Phương pháp lặp Jacobi.  Định nghĩa nghiệm gần đúng:  Ma trận lặp  Biến đổi phương trình nghiệm  Suy ra  Công thức tính lặp: @2009, Ngô Văn Thanh - Viện Vật Lý K=1 ERR=TOL+1.0 DO WHILE (K TOL) DO I=1,N S = 0.0 DO J=1,N S = S-A(I,J)*X1(J) ENDDO S = (S+A(I,N+1))/A(I,I) IF (ABS(S).GT.ERR) ERR=ABS(S) X2(I) = X1(I)+S ENDDO DO I=1,N X1(I) = X2(I) ENDDO K=K+1 ENDDO @2009, Ngô Văn Thanh - Viện Vật Lý 3.2.2. Phương pháp lặp Gauss-Seidel.  Phương pháp được phát triển từ phương pháp Jacobi  Định nghĩa phương trình tính nghiệm gần đúng:  hoặc:  Trong đó ma trận lặp là  Công thức tính lặp: @2009, Ngô Văn Thanh - Viện Vật Lý