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, 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.2. Các phương pháp lặ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

@2009, Ngô Văn Thanh -Viện Vật Lý

3.2.1. Phương pháp lặp Jacobi. 3.2.2. Phương pháp lặp Gauss-Seidel.

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.

@2009, Ngô Văn Thanh -Viện Vật Lý

Biểu diễn dạng ma trận:

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

@2009, Ngô Văn Thanh -Viện Vật Lý

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.



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

@2009, Ngô Văn Thanh -Viện Vật Lý

Giải lần lượt từ phương trình thứ 3 đến phương trình 1.



Trường hợp tổng quát:

@2009, Ngô Văn Thanh -Viện Vật Lý

Chọn hàng thứ nhất làm hàng pivot, ta có:



Thực hiện tiếp cho hàng pivot thứ hai:

@2009, Ngô Văn Thanh -Viện Vật Lý

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.

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).

@2009, Ngô Văn Thanh -Viện Vật Lý

Cuối cùng là tìm các nghiệm bằng phương pháp thay thế ngược (backward substitution)

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ụ:

@2009, Ngô Văn Thanh -Viện Vật Lý

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.



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:

@2009, Ngô Văn Thanh -Viện Vật Lý

Vì a31/s3 lớn nhất nên tráo đổi hàng 1 cho hàng 3.

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

@2009, Ngô Văn Thanh -Viện Vật Lý

: giải phương trình này để tìm x



Giải phương trình:

Phân ly LU để tìm hai ma trận L và U.



Giải phương trình:



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 :

@2009, Ngô Văn Thanh -Viện Vật Lý

Đây là hệ N 2 phương trình với (N 2 + N ) ẩn số.





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.

. Tính 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à

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ý



!No row interchanges yet.

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 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ý

Chương trình đơn giản:

imax=(j-1)+imaxloc(vv(j:n)*abs(a(j:n,j))) if (j /= imax)

do j=1,n

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

@2009, Ngô Văn Thanh -Viện Vật Lý

END SUBROUTINE ludcmp



!No row interchanges yet.

write(6,*) ’ There is a row of zeros.’

vv=1.0/vv

@2009, Ngô Văn Thanh -Viện Vật Lý

Chương trình F90 kiểu parallel tối ưu:

imax=(j-1) + imaxloc(vv(j:n)*abs(a(j:n,j)))

do j=1,n

!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... !...and change the parity of d.

d = -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)

a(j+1:n,j+1:n)=a(j+1:n,j+1:n) - !Divide by the pivot element. &

outerprod(a(j+1:n,j),a(j,j+1:n))

!Reduce remaining submatrix.

end do

@2009, Ngô Văn Thanh -Viện Vật Lý

END SUBROUTINE ludcmp

…

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))

k = i

else if (summ /= 0.0) then

!have to do the dot product above.

end if 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

@2009, Ngô Văn Thanh -Viện Vật Lý

END SUBROUTINE lubksb

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:

@2009, Ngô Văn Thanh -Viện Vật Lý

Á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

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.

enddo call ludcmp(a,indx,d)

y(i,i)=1.

do j=1,n !Decompose the matrix just once. !Find inverse by columns.

call lubksb(a,indx,y(1,j))

@2009, Ngô Văn Thanh -Viện Vật Lý

enddo



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)

@2009, Ngô Văn Thanh -Viện Vật Lý

enddo det = d END FUNCTION det

3.1.5. Các hệ chéo bậc 3.



Chương trình:

@2009, Ngô Văn Thanh -Viện Vật Lý

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

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)

@2009, Ngô Văn Thanh -Viện Vật Lý

end do END SUBROUTINE tridag_ser