Báo cáo bài tập lớn môn Phương pháp tính - Đề tài 5: Giải hệ Ax=b bằng phương pháp Jacobi - Trịnh Thanh Hiệu
ĐẠI HỌC QUỐC GIA TP. HỒ CHÍ MINH
TRƯỜNG ĐẠI HỌC BÁCH KHOA
…………..o0o…………..
BÁO CÁO BÀI TẬP LỚN
MÔN PHƯƠNG PHÁP TÍNH
Đề tài 5: Giải hệ Ax=b bằng phương pháp Jacobi
Giảng viên hướng dẫn: Hoàng Hải Hà
Lớp L02, nhóm 3
Danh sách thành viên
Trịnh Thanh Hiệu
Huỳnh Gia Huy
1611114
1711482
1710107
1711694
1710124
Ngô Trung Huy
Nguyễn Phước Khang
Nguyễn Quốc Khang
BÀI TẬP LỚN
PHƯƠNG PHÁP TÍNH
ĐỀ TÀI 5
MỤC LỤC
I. LỜI NÓI ĐẦU ...................................................................................2
II. ĐỀ TÀI .............................................................................................3
PHẦN 1. CƠ SỞ LÝ THUYẾT ......................................................................3
PHẦN 2. THỰC HIỆN ...................................................................................4
Công cụ sử dụng..............................................................................4
Một số hàm được dùng....................................................................4
Giải thuật.........................................................................................4
PHẦN 3. TÍNH NĂNG VÀ VÍ DỤ ..............................................................11
Các tính năng của chương trình ....................................................11
Một số tính năng khác...................................................................11
Ví dụ..............................................................................................12
III. TÀI LIỆU THAM KHẢO.............................................................15
1
BÀI TẬP LỚN
PHƯƠNG PHÁP TÍNH
ĐỀ TÀI 5
I. LỜI NÓI ĐẦU
Thân chào Thầy cô và các bạn sinh viên!
Đây là báo cáo Bài tập lớn môn học Phương pháp tính do nhóm 3 thực hiện.
Nội dung là giải hệ Ax b bằng phương pháp Jacobi dưới sự hướng dẫn của cô Hoàng
Hải Hà bằng phần mềm MATLAB.
MATLAB là ngôn ngữ lập trình do MathWorks phát triển, cho phép người dùng xây
dựng ma trận, vẽ đồ thị hàm số hay biểu đồ dữ liệu, thực hiện các phép toán, tạo các
giao diện người dùng, liên kết các chương trình máy tính được viết trên nhiều ngôn ngữ
khác nhau, bao gồm cả C, C++, Java, và FORTRAN, phân tích dữ liệu, phát triển các
thuật toán, tạo ra các mô hình và ứng dụng.
MATLAB được tích hợp nhiều lệnh và các hàm toán học, giúp người dùng thực hiện
tính toán các con số, vẽ đồ thị và thực hiện các phương pháp số.
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 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 đến cô Hoàng Hải Hà đã tận
tình hướng dẫn và dạy bảo chúng em trong học kì 2 năm học 2018 - 2019 này.
2
BÀI TẬP LỚN
PHƯƠNG PHÁP TÍNH
ĐỀ TÀI 5
II. ĐỀ TÀI
ĐỀ TÀI 5: Giải hệ Ax b bằng phương pháp JACOBI
1. Kiểm tra sự hội tụ của nghiệm
0
2. Chọn vectơ x tùy ý.
3. Tính vectơ nghiệm x
n
.
4. Đánh giá sai số tiên nghiệm và hậu nghiệm theo cả hai chuẩn.
5. Đánh giá tính ổn định của hệ.
n
6. Tìm chỉ số
n
nhỏ nhất để nghiệm x có sai số nhỏ hơn
cho trước.
PHẦN 1. CƠ SỞ LÝ THUYẾT
- Để 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 ... -a
1n
12
a21 0 ... 0
0
0 ... -a2n
D L U
... ... ... ...
... ... ... ...
an1 an2 ... 0
0
0 ... 0
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 và
Do aii 0,i 1,2,...,n nên det D 0 như vậy tồn tại D1
Khi đó ta có:
Ax b
(D L U)x b
Dx (L U)x b
x D1 *(L U)x D1 *b
1
Tj D *(L U)
Đặt
c D1 *b
j
m
m1
Khi đó thành lập công thức có dạng x Tj x cj
- Kiểm tra tính hội tụ:
_
Nếu Tj 1 thì nghiệm của hệ hội tụ về
x
3
BÀI TẬP LỚN
PHƯƠNG PHÁP TÍNH
ĐỀ TÀI 5
- Công thức đánh giá sai số:
❖ Đánh giá sai số tiên nghiệm
m
_
T
m
1
0
x x
x x
1 T
• Đánh giá sai số hậu nghiệm
_
T
m
m
m1
x x
x x
1 T
PHẦN 2. THỰC HIỆN
• Công cụ sử dụng: Matlab R2016a
• Một số hàm được dùng:
Tên hàm Chức năng
Ví dụ
norm
inv
zeros
Tính chuẩn vectơ và chuẩn ma trận
Tính nghịch đảo của vectơ và ma trận inv(A)
Tạo ma trận 0 A = zeros(5,5)
norm(A,1), norm(A,'inf')
for i = 1:N
…
Lệnh for Vòng lặp
end
If a == 0
….
end
Lệnh if Lệnh điều kiện
clear;clc Xóa dữ liêu, xóa màn hình
• Giải thuật
Đề bài yêu cầu ta giải hệ Ax=b bằng phương pháp lặp JACOBI gồm 6 ý nhỏ. Trong 6 ý
nhỏ trên; ý 1, 2, 5 ta có thể thực hiện ngay khi nhập ma trận A,b; các ý 3, 4 ta cần nhập
số lần lặp (n) và ý thứ 6 ta cần nhập
và chuẩn cần dùng.
Vì vậy chương trình sẽ đi theo 2 hướng cho người dùng lựa chọn: một là tìm Xn, đánh
giá sai số; hai là tìm chỉ số n nhỏ nhất để nghiệm Xn có sai số nhỏ hơn
cho trước.
code = input('Ban muon chuong trinh thuc hien dieu gi? \n
1: Tim Xn, danh gia sai so \n
2: Tim chi so n nho nhat de nghiem Xn co sai so nho hon eps cho truoc\n
0: Thoat\n
Nhap: ');
✓ Ý thứ 1: Kiểm tra sự hội tụ của nghiệm.
- Để kiểm tra sự hội tụ của nghiệm, ta phải tìm được Tj và xét nếu Tj 1 thì nghiệm
_
của hệ hội tụ về
x.
- Tìm Tj: Xét ma trận có phải ma trận đường chéo nghiêm ngặt hay không (det(A)
khác 0, các phần tử trên đường chéo khác 0).
4
BÀI TẬP LỚN
PHƯƠNG PHÁP TÍNH
ĐỀ TÀI 5
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;
-
Tìm các ma trận D, L, U như cơ sở lý thuyết:
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;
-
Tính Tj, cj:
Tj = inv(D)*(L+U);
cj = inv(D)*b;
-
Tính Tj và xét điều kiện hội tụ của nghiệm:
if norm(Tj,'inf') < 1
disp('Nghiem cua he hoi tu ');
else
disp('Nghiem cua he khong hoi tu ');
end;
0
✓ Ý thứ 2: Chọn vectơ x tùy ý.
0
-
x có thể là vectơ được nhập vào hoặc random, nhập 0 để chọn vectơ 0 hoặc 1 để
tạo vectơ ngẫu nhiên.
if X0 == 0
X0 = zeros(N,1);
end;
if X0 == 1
X0 = rand(N,1);
end;
n
✓ Ý thứ 3,4: Tính vectơ nghiệm x - đánh giá sai số tiên nghiệm và hậu nghiệm
theo cả hai chuẩn.
-
Nhập số lần lặp tối đa, nếu số lần lặp nhỏ hơn 1 thì yêu cầu người dùng nhập lại.
if code == 1
maxlap = input('Nhap so lan lap: ');
while maxlap < 1
maxlap = input('So lan lap phai lon hon 0, moi ban nhap lai: ');
end;
end;
-
-
Gán giá trị ban đầu: số lần đã lặp n = 0 và Xn = X0, codec = 0 (sử dụng sai số hậu
nghiệm chuẩn 1 để so sánh với eps)
Tính sai số theo:
5
BÀI TẬP LỚN
PHƯƠNG PHÁP TÍNH
ĐỀ TÀI 5
❖ Tiên nghiệm chuẩn 1
m
_
T
m
1
0
x x
x x
1
1
1 T
1
1
❖ Tiên nghiệm chuẩn vô cùng
m
_
T
m
1
0
x x
x x
1 T
❖ Hậu nghiệm chuẩn 1
_
T
m
m
m1
x x
x x
1
1
1 T
1
1
❖ Hậu nghiệm chuẩn vô cùng
_
T
m
m
m1
x x
x x
1 T
%Gan cac gia tri ban dau.
n = 0;
Xn=X0;
for j = 1:maxlap
Xn2 = Xn;
Xn = Tj*Xn2 + cj;
n = n + 1;
%sai so tien nghiem chuan 1
TienNgChuan1 = abs((norm(Tj,1)^n)*norm(X1-X0,1)/(1-norm(Tj,1)));
%sai so tien nghiem chuan vo cung
TienNgChuanVoCung = abs((norm(Tj,'inf')^n)*norm(X1-X0,'inf')/(1-
norm(Tj,'inf')));
%sai so hau nghiem chuan 1
HauNgChuan1 = abs(norm(Tj,1)*norm(Xn-Xn2,1)/(1-norm(Tj,1)));
%sai so hau nghiem chuan vo cung
HauNgChuanVoCung = abs(norm(Tj,'inf')*norm(Xn-Xn2,'inf')/(1-
norm(Tj,'inf')));
if codec == 0
saiso = HauNgChuan1;
end;
if saiso < eps
break;
end;
end;
✓ Ý thứ 5: Đánh giá tính ổn định của hệ.
-
Số điều kiện k A 1 * A1 càng gần 1 thì nghiệm của hệ càng ổn định, ngược lại
1
k càng lớn thì nghiệm của hệ càng không ổn định
%tim so dieu kien
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;
n
✓ Ý thứ 6: Tìm chỉ số
- Để tìm , ta cần nhập
n
nhỏ nhất để nghiệm x có sai số nhỏ hơn
và điều kiện cần dùng
cho trước.
n
6
BÀI TẬP LỚN
PHƯƠNG PHÁP TÍNH
ĐỀ TÀI 5
codec = 0;
if code == 2
eps = input('Moi ban nhap eps: ');
codec = input('Ban muon su dung dieu kien gi??\n
2: Xn - Xn-1, chuan vo cuc\nNhap: ');
1: Xn - Xn-1, chuan
1\n
end;
- Sau đó so sánh với eps để thoát khỏi vòng lặp
if codec == 1
saiso = norm(Xn-Xn2,1);
end;
if codec == 2
saiso = norm(Xn-Xn2,'inf');
end;
if saiso < eps
break;
end;
7
BÀI TẬP LỚN
PHƯƠNG PHÁP TÍNH
ĐỀ TÀI 5
• Mã nguồn
% -------------------------------------------------------------------------
% De tai 5: Giai he Ax = b bang phuong phap lap Jacobi
% ---------------------------******------------------------
% 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
% -------------------------------------------------------------------------
function Jacobi(N,A,b,X0)
clc;
disp('------------------------------------------------');
disp('Giai he Ax = b bang phuong phap lap Jacobi');
disp('----------------------******--------------------');
if nargin == 0
N = input('Nhap N: ');
A = input('Nhap ma tran A: ');
b = input('Nhap ma tran b: ');
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 Jacobi');
disp('----------------------******--------------------');
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;
8
BÀI TẬP LỚN
PHƯƠNG PHÁP TÍNH
ĐỀ TÀI 5
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;
Tj = inv(D)*(U+L);
cj = inv(D)*b;
% Xet tinh hoi tu
if norm(Tj,'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,
2: Tim chi so n nho nhat de nghiem Xn co sai so nho
0: Thoat\nNhap: ');
danh gia sai so \n
hon eps cho truoc\n
if code == 1
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 = Tj*X0+cj;
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: ');
Xn=X0;
for j = 1:maxlap
Xn2 = Xn;
Xn = Tj*Xn2 + cj;
n = n + 1;
%sai so tien nghiem chuan 1
TienNgChuan1 = abs((norm(Tj,1)^n)*norm(X1-X0,1)/(1-norm(Tj,1)));
%sai so tien nghiem chuan vo cung
TienNgChuanVoCung = abs((norm(Tj,'inf')^n)*norm(X1-X0,'inf')/(1-
norm(Tj,'inf')));
9
BÀI TẬP LỚN
PHƯƠNG PHÁP TÍNH
ĐỀ TÀI 5
%sai so hau nghiem chuan 1
HauNgChuan1 = abs(norm(Tj,1)*norm(Xn-Xn2,1)/(1-norm(Tj,1)));
%sai so hau nghiem chuan vo cung
HauNgChuanVoCung = abs(norm(Tj,'inf')*norm(Xn-Xn2,'inf')/(1-
norm(Tj,'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;
10
BÀI TẬP LỚN
• Test case
STT N
PHƯƠNG PHÁP TÍNH
ĐỀ TÀI 5
Số
Sai
số
Kết
quả
A
b
X0
lần
lặp
Yêu cầu
[10,-1,2,0;
-1,11,-
1,3;2,-
1,10,-1;
0,3,-1,8]
0.9890
2.0114
-
1.0103
1.0214
[6;25;-
11;15]
5
Tính x
1
1
4
2
[0;0;0;0]
[0.7;0.4]
5
Tính chỉ số n nhỏ
nhất để
[9,-7;-3,7]
[2;5]
0.06
n = 7
n
n1
x x
0.06
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 Jacobi.
- Giao diện trình bày dễ sử dụng.
- Độ chính xác cao.
✓ Tiêu cực:
- Việc nhập dữ 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
• Chọn vectơ x tùy ý.
n
• Tính vectơ nghiệm x
.
• Đánh giá sai số tiên nghiệm và 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 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.
• 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.
Ví dụ:
a. Ví dụ 1: Trong Giáo trình Phương Pháp Tính – Lê Thái Thanh trang 56 có bài:
11
BÀI TẬP LỚN
PHƯƠNG PHÁP TÍNH
ĐỀ TÀI 5
Giải hệ sau bằng phương pháp lặp Jacobi.
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
3 1
0
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 Jacobi.m):
>>Jacobi(4,[10,-1,2,0; -1,11,-1,3;2,-1,10,-1; 0,3,-1,8],[6;25;-11;15],0)
Số lần lặp: 5
Sau đây là màn hình khi chạy chương trình:
------------------------------------------------
Giai he Ax = b bang phuong phap lap Jacobi
----------------------******--------------------
N =
4
A =
10
-1
2
-1
11
-1
3
2
-1
10
-1
0
3
-1
8
0
b =
6
25
-11
15
12
BÀI TẬP LỚN
PHƯƠNG PHÁP TÍNH
ĐỀ TÀI 5
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 =
0.9890
2.0114
-1.0103
1.0214
TienNgChuan1 =
0.8648
TienNgChuanVoCung =
0.1420
HauNgChuan1 =
0.2349
HauNgChuanVoCung =
0.0577
Ban muon tiep tuc?
So bat ky: Tiep tuc
0: Thoat
Nhap: 0
****************CHUONG TRINH KET THUC*********************
>>
Kết quả:
Xn =
0.9890
2.0114
-1.0103
1.0214
13
BÀI TẬP LỚN
PHƯƠNG PHÁP TÍNH
ĐỀ TÀI 5
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 ví dụ này, ta xác định được:
16 2
A
6 13
2
b
5
0.4
X 0
0.8
Sai số: 0.007
Để giải hệ này, ta nhập vào Matlab ở ô Comman Window (Set Path tại thư mục
chứa file Jacobi.m):
>>Jacobi(2,[16,2;-6,13],[2;5],[0.4;0.8])
Khi hỏi sai số, ta nhập 0.007
Đây là màn hình khi ta chạy chương trình
------------------------------------------------
Giai he Ax = b bang phuong phap lap Jacobi
----------------------******--------------------
N =
2
A =
16
-6
2
13
b =
2
5
X0 =
0.4000
0.8000
Nghiem cua he hoi tu
So dieu kien: 1.900000
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
14
BÀI TẬP LỚN
PHƯƠNG PHÁP TÍNH
ĐỀ TÀI 5
Nhap: 2
Moi ban nhap eps: 0.007
Ban muon su dung dieu kien gi??
1: Xn - Xn-1, chuan 1
2: Xn - Xn-1, chuan vo cuc
Nhap: 2
n =
5
Ban muon tiep tuc?
So bat ky: Tiep tuc
0: Thoat
Nhap: 0
****************CHUONG TRINH KET THUC*********************
>>
Kết quả :
n =
5
III. 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
15
16 trang
20400
Free
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 5: Giải hệ Ax=b bằng phương pháp Jacobi - Trịnh Thanh Hiệu", để 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:
- bao_cao_bai_tap_lon_mon_phuong_phap_tinh_de_tai_5_giai_he_ax.doc
PPT_ChuDe5.pdf
