Báo cáo bài tập lớn môn Phương pháp tính - Đề tài 6: Giải hệ bằng Ax=b phương pháp Gauss-Seidel

ĐẠI HỌC QUỐC GIA TP. HỒ CHÍ MINH  
TRƯỜNG ĐẠI HỌC BÁCH KHOA  
…………..o0o…………..  
BÁO CÁO BTL  
PHƯƠNG PHÁP TÍNH  
Giáo viên hướng dẫn: Hoàng Hải Hà  
A x b  
Đề tài 6: Giải hệ  
bằng  
phương pháp Gauss-Seidel  
Lớp L06, Nhóm 15  
Danh sách thành viên  
Lê Hoàng Dương  
Đặng Lê Thanh Hiếu  
Thái Hải Lâm  
Huỳnh Minh Thuận  
Nguyễn Duy Bảo  
Thị Thúy Quỳnh  
1710900  
1711274  
1711905  
1710315  
1710592  
1712922  
Bài tập lớn  
PHƯƠNG PHÁP TÍNH  
Nhóm 15 – Đề tài 6  
Lời nói đầu  
Thân chào Thầy cô và các bạn sinh viên!  
Đây quyển báo cáo Bài tập lớn do Nhóm 15 thực hiện.  
Nội dung là giải hệ Ax b bằng phương pháp Gauss-Seidel dưới sự hướng dẫn  
của cô ThS. Hoàng Hải Hà.  
BÀI BÁO CÁO GỒM CÁC PHẦN  
Nhóm chúng em đã cố gắng trình bày nổi bật các ý chính, cụ thể các hàm và cung  
cấp TestCase để bạn đọc thể dễ dàng hiểu rõ và đánh giá.  
Thay mặt cả lớp, Chúng em gửi lời cảm ơn chân thành nhất cô ThS. Hoàng Hải  
đã tận tình hướng dẫn dạy bảo chúng em trong học kì 1 năm học 2018 này.  
1
Bài tập lớn  
PHƯƠNG PHÁP TÍNH  
Nhóm 15 – Đề tài 6  
ĐỀ TÀI  
Ax b  
ĐỀ TÀI 6: Giải hệ  
bằng phương pháp Gauss-Seidel  
Kiểm tra sự hội tụ của nghiệm  
0
Chọn vectơ x  tùy ý.  
n
Tính vectơ nghiệm x  
.
Đánh giá sai số tiên nghiệm hậu nghiệm theo cả hai chuẩn.  
Đánh giá tính ổn định của hệ.  
n
Tìm chỉ số  
n
nhỏ nhất để nghiệm xcó sai số nhỏ hơn  
cho trước.  
PHẦN 1. CƠ SỞ THUYẾT  
-
Trong giải tích số, phương pháp Gauss-Seidel hay còn gọi phương  
pháp lặp Gauss-Seidel, phương pháp Liebmann hay phương pháp tự sửa sai là  
một phương pháp lặp được sử dụng để giải một hệ phương trình tuyến  
tính tương tự như phương pháp Jacobi. Nó được đặt tên theo hai nhà toán  
phương pháp này có thể áp dụng cho bất kỳ ma trận nào không chứa phần  
tử 0 (không) trên các đường chéo, nhưng tính hội tụ chỉ xảy ra nếu ma trận  
-
Để giải hệ Ax b ta phân tích  
a
a12 ... a1n  
a
0 ... 0  
   
   
11  
11  
a21 a22 ... a2n  
0
a22 ... 0  
   
A   
   
   
   
... ... ... ...  
... ... ... ...  
an1 an2 ... ann  
0
0 ... ann  
0
0 ... 0  
0 -a ... -a1n  
12  
   
   
a21 0 ... 0  
0
0 ... -a2n  
   
   
   
   
... ... ... ...  
... ... ... ...  
an1 -an2 ... 0  
0
0 ... 0  
D L U  
Với điều kiên giả sử  
aii 0,i 1,2,...,n  
A
là ma trận đường chéo trội nghiêm ngặt tức det A 0  
Do aii 0,i 1,2,...,n nên det D 0 như vậy tồn tại D1 cũng tồn tại  
(D L)1  
Khi đó ta có:  
2
   
Bài tập lớn  
PHƯƠNG PHÁP TÍNH  
Nhóm 15 – Đề tài 6  
Ax b  
(D L U)x b  
(D L)x Ux b  
x (D L)1 *Ux (D L)1b  
Đặt  
Tg (D L)1 *U  
cg (D L)1b  
Khi đó thành lập công thức dạng  
m
m1  
xTg x  cg  
-
Kiểm tra tính hội tụ:  
_
Nếu Tg 1 thì nghiệm của hệ hội tụ về  
Công thức đánh giá sai số:  
x
-
Đánh giá sai số tiên nghiệm  
m
_
T
m
1
0
xx   
x  x   
1T  
Đánh giá sai số hậu nghiệm  
_
T
m
m
m1  
xx   
xx  
1T  
3
Bài tập lớn  
PHƯƠNG PHÁP TÍNH  
Nhóm 15 – Đề tài 6  
PHẦN 2. HIỆN THỰC  
Công cụ sử dụng: Matlab 2016a  
Một số hàm được dùng:  
Tên hàm  
Chức năng  
dụ  
norm  
inv  
zeros  
Tính chuẩn vectơ chuẩn ma trận  
Tính nghịch đảo của vectơ và ma trận  
Tạo ma trận 0  
norm(A,1), norm(A,'inf')  
int(A)  
A = zeros(5,5)  
for i = 1:N  
Lệnh for  
Vòng lặp  
end  
If a == 0  
….  
Lệnh if  
Lệnh điều kiện  
end  
clear;clc  
Xóa dữ liêu, xóa màn hình  
Source Code  
% -------------------------------------------------------------------------  
% De tai 6: Giai he Ax = b bang phuong phap lap GaussSeidel  
% ---------------------------******------------------------  
% INPUT:  
%
%
%
%
%
N la cap cua ma tran he so  
Cac ma tran A,b la ma tran he so cua he Ax = b  
X0 là vectơ lap ban dau (nhap 0 de chon vecto 0, nhap 1 de chon random)  
eps là sai so (gia tri mac dinh là 1.0E-6)  
maxlap là so lan lap toi da cho phep (gia tri mac dinh la 100)  
% OUTPUT:  
%
%
%
%
%
%
Xn la vecto nghiem  
TienNgChuan1 la sai so tien nghiem chuan 1  
TienNgChuanVoCung la sai so tien nghiem chuan vo cung  
HauNgChuan1 la sai so hau nghiem chuan 1  
HauNgChuanVoCung la sai so hau nghiem chuan vo cung  
n la so lan lap thoa man yeu cau  
% TEST:  
% Test 1  
%
GaussSeidel(4,[10,-1,2,0; -1,11,-1,3;2,-1,10,-1; 0,3,-1,8],[6;25;-  
11;15],0)  
%
%
%
%
%
%
%
%
%
%
N = 4  
A = [10,-1,2,0; -1,11,-1,3;2,-1,10,-1; 0,3,-1,8]  
b = [6;25;-11;15]  
X0 = 0 (auto X0 = [0;0;0;0])  
so lan lap: 5  
Ket qua: Xn =  
1.0001  
2.0000  
-1.0000  
1.0000  
% Test 2  
%
%
%
%
%
GaussSeidel(2,[9,-7;-3,7],[2;5],[0.7;0.4])  
N = 2  
A = [9,-7;-3,7]  
b = [2;5]  
X0 = [0.7;0.4]  
4
 
Bài tập lớn  
PHƯƠNG PHÁP TÍNH  
Nhóm 15 – Đề tài 6  
%
%
esp = 0.06 ( chuan 1)  
Ket qua: n = 5  
% Test 3  
%
%
%
%
%
%
%
%
%
GaussSeidel(2,[11,5;-3,11],[2;4],[0.9;0.2])  
N = 2  
A = [11,5;-3,11]  
b = [2;4]  
X0 = [0.9;0.2]  
so lan n: 3  
Ket qua: Xn =  
0.0159  
0.3680  
% Test 4  
%
%
%
%
%
%
%
GaussSeidel(2,[15,3;6,13],[6;2],[0.2;0.2])  
N = 2  
A = [15,3;6,13]  
b = [6;2]  
X0 = [0.2;0.2]  
esp = 0.007 ( chuan 1)  
Ket qua: n = 3  
% -------------------------------------------------------------------------  
function GaussSeidel(N,A,b,X0)  
clc;  
disp('------------------------------------------------');  
disp('Giai he Ax = b bang phuong phap lap GaussSeidel');  
disp('----------------------******--------------------');  
if nargin == 0  
N = input('Nhap N: '); if N == 0 return; end;  
A = input('Nhap ma tran A: '); if A == 0 return; end;  
b = input('Nhap ma tran b: '); if b == 0 return; end;  
X0 = input('Nhap X0: ');  
end;  
if nargin == 1  
A = input('Nhap ma tran A: '); if A == 0 return; end;  
b = input('Nhap ma tran b: '); if b == 0 return; end;  
X0 = input('Nhap X0: ');  
end;  
if nargin == 2  
b = input('Nhap ma tran b: '); if b == 0 return; end;  
X0 = input('Nhap X0: ');  
end;  
if nargin == 3  
X0 = input('Nhap X0: ');  
end;  
maxlap = 100;  
eps = 1.0E-6;  
% xu li X0  
if X0 == 0  
X0 = zeros(N,1);  
end;  
if X0 == 1  
X0 = rand(N,1);  
end;  
code = 3;  
while code ~= 0  
clc;  
disp('------------------------------------------------');  
disp('Giai he Ax = b bang phuong phap lap GaussSeidel');  
disp('----------------------******--------------------');  
5
Bài tập lớn  
PHƯƠNG PHÁP TÍNH  
Nhóm 15 – Đề tài 6  
N
A
b
X0  
% Xet ma tran co phai ma tran duong cheo nghiem ngat hay khong?  
if det(A) == 0, disp('Ma tran da nhap khong phai ma tran duong cheo nghiem  
ngat.'); return; end;  
for i=1:N  
if A(i,i) == 0, disp('Ma tran da nhap khong phai ma tran duong cheo  
nghiem ngat.');return; end;  
end;  
D = zeros(N,N);  
for i=1:N D(i,i)= A(i,i); end;  
L = zeros(N,N);  
for i=2:N  
for j=1:i-1  
L(i,j) = -A(i,j);  
end;  
end;  
U = zeros(N,N);  
for i=N-1:-1:1  
for j=N:-1:i+1  
U(i,j) = - A(i,j);  
end;  
end;  
Tg = inv(D-L)*U;  
cg = inv(D-L)*b;  
% Xet tinh hoi tu  
if norm(Tg,'inf') < 1  
disp('Nghiem cua he hoi tu ');  
else  
disp('Nghiem cua he khong hoi tu ');  
end;  
k1 = norm(A,1)*norm(inv(A),1);  
fprintf('So dieu kien: %f\n',k1);  
if k1<15 disp('He on dinh'); else disp('He khong on dinh'); end;  
code = input('Ban muon chuong trinh thuc hien dieu gi? \n  
1: Tim Xn,  
danh gia sai so \n  
hon eps cho truoc\n  
if code == 1  
2: Tim chi so n nho nhat de nghiem Xn co sai so nho  
0: Thoat\nNhap: ');  
maxlap = input('Nhap so lan lap: ');  
while maxlap < 1  
maxlap = input('So lan lap phai lon hon 0, moi ban nhap lai: ');  
end;  
end;  
n = 0;  
X1 = Tg*X0+cg;  
codec = 0;  
if code == 2  
eps = input('Moi ban nhap eps: ');  
codec = input('Ban muon su dung dieu kien gi??\n  
1: Xn - Xn-1, chuan  
1\n  
end;  
2: Xn - Xn-1, chuan vo cuc\nNhap: ');  
6
Bài tập lớn  
PHƯƠNG PHÁP TÍNH  
Nhóm 15 – Đề tài 6  
Xn=X0;  
for j = 1:maxlap  
Xn2 = Xn;  
Xn = Tg*Xn2 + cg;  
n = n + 1;  
%sai so tien nghiem chuan 1  
TienNgChuan1 = abs((norm(Tg,1)^n)*norm(X1-X0,1)/(1-norm(Tg,1)));  
%sai so tien nghiem chuan vo cung  
TienNgChuanVoCung = abs((norm(Tg,'inf')^n)*norm(X1-X0,'inf')/(1-  
norm(Tg,'inf')));  
%sai so hau nghiem chuan 1  
HauNgChuan1 = abs(norm(Tg,1)*norm(Xn-Xn2,1)/(1-norm(Tg,1)));  
%sai so hau nghiem chuan vo cung  
HauNgChuanVoCung = abs(norm(Tg,'inf')*norm(Xn-Xn2,'inf')/(1-  
norm(Tg,'inf')));  
if codec == 0  
saiso = HauNgChuan1;  
end;  
if codec == 1  
saiso = norm(Xn-Xn2,1);  
end;  
if codec == 2  
saiso = norm(Xn-Xn2,'inf');  
end;  
if saiso < eps  
break;  
end;  
end;  
% Output  
if code == 1  
Xn  
codes = input('Ban co muon xuat sai so khong? \n  
1: Co\n  
2:  
Khong\nNhap: ');  
if codes == 1  
TienNgChuan1  
TienNgChuanVoCung  
HauNgChuan1  
HauNgChuanVoCung  
end;  
code = input('Ban muon tiep tuc?\n  
So bat ky: Tiep tuc\n  
0:  
Thoat\nNhap: ');  
end;  
if code == 2  
n
code = input('Ban muon tiep tuc?\n  
Thoat\nNhap: ');  
So bat ky: Tiep tuc\n  
0:  
end;  
end;  
disp('****************CHUONG TRINH KET THUC*********************');  
return;  
7
Bài tập lớn  
PHƯƠNG PHÁP TÍNH  
Nhóm 15 – Đề tài 6  
Test case  
Số lần  
lặp  
Sai  
số  
STT  
N
A
b
X0  
Yêu cầu  
Kết quả  
[10,-1,2,0; -  
1,11,-1,3;2,-  
1,10,-1; 0,3,-  
1,8]  
1.0001  
2.0000  
-1.0000  
1.0000  
5
Tính x   
1
4
[6;25;-11;15]  
[0;0;0;0]  
5
Tính chỉ số n nhỏ nhất để  
n
n1  
1
3
4
2
2
2
[9,-7;-3,7]  
[11,5;-3,11]  
[15,3;6,13]  
[2;5]  
[2;4]  
[6;2]  
[0.7;0.4]  
[0.9;0.2]  
[0.2;0.2]  
0.06  
n = 5  
x  x  
0.06  
1
0.0159  
0.3680  
3
Tính x   
3
Tính chỉ số n nhỏ nhất để  
n
n1  
0.007  
n = 3  
x  x  
0.007  
1
Một số đánh giá:  
Tích cực:  
-
-
-
Code đã giải quyết hầu hết các vấn đề về phương pháp Gauss - Seidel  
Giao diện trình bày dễ sử dụng  
Độ chính xác cao  
Tiêu cực:  
-
-
Việc nhập liệu dễ sai sót  
Code chưa thật sự tối ưu  
PHẦN 3. TÍNH NĂNG VÀ VÍ DỤ  
Các tính năng của chương trình:  
Kiểm tra sự hội tụ của nghiệm  
0
x  
Chọn vectơ tùy ý.  
n
x  
Tính vectơ nghiệm  
.
Đánh giá sai số tiên nghiệm hậu nghiệm theo cả hai chuẩn.  
Đánh giá tính ổn định của hệ.  
Tìm chỉ số nhỏ nhất để nghiệm  
n
n
x  
có sai số nhỏ hơn cho trước.  
Một số tính năng khác:  
Kiểm tra ma trận nhập vào có phải ma trận đường chéo nghiêm ngặt hay không  
Nếu nhập vào số lần lập lặp < 1 thì chương trình sẽ yêu cầu nhập lại  
Chương trình thiết kế thể tự nhập hoặc nhập dưới dạng gọi hàm.  
Cho phép người dùng nhập nhanh vectơ X0 với: Nhập 0 để chọn vectơ 0 hoặc 1 để tạo  
vectơ ngẫu nhiên  
dụ  
a. Ví dụ 1:  
Trong Giáo trình Phương Pháp Tính – Lê Thái Thanh trang 59 có bài:  
Giải hệ sau bằng phương pháp lặp Gauss-Seidel  
8
 
Bài tập lớn  
PHƯƠNG PHÁP TÍNH  
Nhóm 15 – Đề tài 6  
10x x 2x  
6  
1
2
3
x 11x x 3x 25  
1
2
3
4
2x1 x2 10x3 x4 11  
3x2 x3 8x4 15  
Từ hệ ta có:  
10 1 2  
0
1 11 1 3  
A   
b   
2 11 0 1  
0
3 1 8  
6
25  
11  
15  
0
   
   
0
   
X 0   
   
0
   
0
   
Để giải hệ này, ta nhập vào Matlab ô Comman Window (Set Path tại thư mục chứa file  
GaussSeidel.m):  
>>GaussSeidel(4,[10,-1,2,0; -1,11,-1,3;2,-1,10,-1; 0,3,-1,8],[6;25;-  
11;15],0)  
Hoặc chạy chương trình(f5) và nhập từng bước:  
N = 4  
A = [10,-1,2,0; -1,11,-1,3;2,-1,10,-1; 0,3,-1,8]  
b = [6;25;-11;15]  
X0 = 0 (auto X0 = [0;0;0;0])  
Số lần lặp: 5  
Ta được kết quả:  
Xn =  
1.0001  
2.0000  
-1.0000  
1.0000  
Sau đây là màn hình khi chạy chương trình:  
------------------------------------------------  
Giai he Ax = b bang phuong phap lap GaussSeidel  
----------------------******--------------------  
N =  
4
A =  
9
Bài tập lớn  
PHƯƠNG PHÁP TÍNH  
Nhóm 15 – Đề tài 6  
10  
-1  
2
-1  
11  
-1  
3
2
-1  
10  
-1  
0
3
-1  
8
0
b =  
6
25  
-11  
15  
X0 =  
0
0
0
0
Nghiem cua he hoi tu  
So dieu kien: 3.137255  
He on dinh  
Ban muon chuong trinh thuc hien dieu gi?  
1: Tim Xn, danh gia sai so  
2: Tim chi so n nho nhat de nghiem Xn co sai so nho hon eps cho truoc  
0: Thoat  
Nhap: 1  
Nhap so lan lap: 5  
Xn =  
1.0001  
2.0000  
-1.0000  
1.0000  
Ban co muon xuat sai so khong?  
1: Co  
2: Khong  
Nhap: 1  
TienNgChuan1 =  
0.1756  
TienNgChuanVoCung =  
0.0202  
HauNgChuan1 =  
0.0012  
HauNgChuanVoCung =  
10  
Bài tập lớn  
PHƯƠNG PHÁP TÍNH  
Nhóm 15 – Đề tài 6  
4.2279e-04  
Ban muon tiep tuc?  
So bat ky: Tiep tuc  
0: Thoat  
Nhap: 0  
****************CHUONG TRINH KET THUC*********************  
>>  
Kết quả:  
Xn =  
1.0001  
2.0000  
-1.0000  
1.0000  
b. Ví dụ 2  
Trong đề thi giữa kì PPT của Trường Đại Học Bách Khoa năm 2017 có câu  
Với dụ này, ta xác định được:  
9 7  
3 7  
A   
b   
2
   
   
5
   
0.7  
0.4  
X 0   
Sai số: 0.06  
Để giải hệ này, ta nhập vào Matlab ô Comman Window (Set Path tại thư mục chứa file  
GaussSeidel.m):  
>>GaussSeidel(2,[9,-7;-3,7],[2;5],[0.7;0.4])  
Hoặc chạy chương trình (f5) và nhập từng bước:  
N = 2  
A = [9,-7;-3,7]  
b = [2;5  
X0 = [0.7;0.4]  
Khi hỏi sai số, ta nhập 0.06  
Kết quả: n = 5  
Đây là màn hình khi ta chạy chương trình  
------------------------------------------------  
Giai he Ax = b bang phuong phap lap GaussSeidel  
----------------------******--------------------  
N =  
2
11  
Bài tập lớn  
PHƯƠNG PHÁP TÍNH  
Nhóm 15 – Đề tài 6  
A =  
9
-3  
-7  
7
b =  
2
5
X0 =  
0.7000  
0.4000  
Nghiem cua he hoi tu  
So dieu kien: 5.333333  
He on dinh  
Ban muon chuong trinh thuc hien dieu gi?  
1: Tim Xn, danh gia sai so  
2: Tim chi so n nho nhat de nghiem Xn co sai so nho hon eps cho truoc  
0: Thoat  
Nhap: 2  
Moi ban nhap eps: 0.06  
Ban muon su dung dieu kien gi??  
1: Xn - Xn-1, chuan 1  
2: Xn - Xn-1, chuan vo cuc  
Nhap: 1  
n =  
5
Ban muon tiep tuc?  
So bat ky: Tiep tuc  
0: Thoat  
Nhap: 0  
****************CHUONG TRINH KET THUC*********************  
>>  
Kết quả :  
n =  
5
TÀI LIỆU THAM KHẢO  
Giáo trình Phương Pháp Tính – Lê Thái Thanh – Nhà xuất bản ĐHQG TP.HCM  
12  
 
doc 13 trang yennguyen 31/03/2022 6620
Bạn đang xem tài liệu "Báo cáo bài tập lớn môn Phương pháp tính - Đề tài 6: Giải hệ bằng Ax=b phương pháp Gauss-Seidel", để tải tài liệu gốc về máy hãy click vào nút Download ở trên

File đính kèm:

  • docbao_cao_bai_tap_lon_mon_phuong_phap_tinh_de_tai_6_giai_he_ba.doc