Bài giảng Applied numerical methods (Ứng dụng phương pháp tính số): Chương 3 - TS. Ngô Văn Thanh

Chia sẻ: Little Little | Ngày: | Loại File: PDF | Số trang:28

Thêm vào BST

Báo xấu

110
lượt xem 14
download

Download Vui lòng tải xuống để xem tài liệu đầy đủ

Chương 3 - Hệ phương trình đại số tuyến tính. Chương này cung cấp cho người học những kiến thức về phương pháp tính trực tiếp và các phương pháp lặp. Nội dung cụ thể gồm: Phương pháp khử Gaussian, quay và tính tỷ lệ - Pivoting and scaling, phương pháp phân ly LU, nghịch đảo ma trận, các hệ chéo bậc 3, phương pháp lặp Jacobi, phương pháp lặp Gauss-Seidel.

Chủ đề:

Bình luận(0) Đăng nhập để gửi bình luận!

Lưu

Nội dung Text: Bài giảng Applied numerical methods (Ứng dụng phương pháp tính số): Chương 3 - TS. 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ý