TRƯỜNG ĐẠI HỌC BÁCH KHOA HÀ NỘI
KHOA TOÁN TIN ỨNG DỤNG
----- -----

TIỂU LUẬN
Đề tài: Tìm nghiệm xấp xỉ hệ phương trình tuyến tính đại số
Giáo viên hướng dẫn: Nguyễn Hữu Điển
Sinh viên thực hiện

: Lê Minh Cường

Lớp: Toán – Tin 1-k51.

Hà Nội, tháng 11 năm 2009


LỜI NÓI ĐẦU
Nội dung chủ yếu tập trung vào phần giải thuật và cách cài đặt nó trên chương trình ứng
dụng Maple, để giải quyết bài toán tìm nghiệm xấp xỉ hệ phương trình tuyến tính đại số
Ax = b bằng các phương pháp như Jacobi hay Gauss_Seidel …
Phiên bản sử dụng là Maple 12.Ở đây xin được trình bày 5 thuật toán chính và cài đặt
3 thuật toán đầu tiên bằng Maple 12.
Tuy nhiên do khả năng đọc hiểu tiếng Anh còn hạn chế nên có thể có những chỗ giải
thích chưa thoả đáng,rất mong được sự góp ý của thầy cô và bạn bè.
Xin chân thành cám ơn sự giúp đỡ chỉ bảo tận tình của phó giáo sư Nguyễn Hữu
Điển.
Mọi đóng góp xin gửi về địa chỉ email : hoặc
.
Một lần nữa xin chân thành cám ơn.
Sinh viên trình bày: Lê Minh Cường
Lớp : Toán Tin 1 – K51


1.3. Định lý Causi-Bunhia-Svac :
t
t
Cho x = ( x1 ,..., xn ) và y = ( y1 ,..., yn ) trong không gian R n .Ta có bất đẳng
thức :
n

n

xt y = ∑ xi yi ≤

∑ xi2

i =1

i =1

n

∑y

2
i

i =1

= x. y

1.4. Định nghĩa 4 : Cho x = ( x1 ,..., xn ) và y = ( y1 ,..., yn ) trong không gian

∞

chuẩn ⋅ khi ∀ε > 0 ,tồn tại số nguyên N (ε ) thỏa mãn :
x (k ) − x < ε ; ∀k ≥ N (ε )

1.6. Định lý : Dãy vecto {x (k )}k =1 hội tụ đến x trong R n với chú ý là chuẩn . ∞ ,khi
∞

và chỉ khi lim k → ∞ xi(k ) = xi ,với i = 1,..., n
1.7. Định lý :
Với x ∈ R n thì x ∞ ≤ x 2 ≤ n x

∞

2.Chuẩn Ma Trận :
2.1. Định nghĩa : Chuẩn của ma trận trên tập các ma trận vuông n chiều là một hàm số
,định nghĩa trên tập,thỏa mãn tất cả ma trận n × n là A,B và số thực
thực,ký hiệu

α
i) A ≥ 0
3


ii) A = 0 ⇔ A = O
iii) αA = α A
iv) A + B ≤ A + B
v) A.B ≤ A . B
2.1. Định lý : Chuẩn của ma trận được xác định như sau
Az

> evalm(B);

4


4. Định nghĩa 3 : Bán kính quy chiếu ρ ( A) của ma trận A được định nghĩa là

ρ ( A) = max λ với λ là trị riêng của A

Đối với λ là giá trị phức thì λ = α + βi = α 2 + β 2
5. Định lý : Nếu A là ma trận vuông n chiều thì :

( )

i) A 2 = ρ At A
ii) ρ ( A) ≤ A

6.Ví dụ :

⎡ 1 1 0⎤
Cho A = ⎢⎢ 1 2 1⎥⎥ Tìm A 2
⎢⎣− 1 1 2⎥⎦
>with(linalg);
>A:=matrix(3,3, [1,1,0,1,2,1,-1,1,2]);

Dùng hàm chuyển vị transpose để tính At
>B:=transpose(A);

Dùng hàm nhân multiply để tính At A
>C:=multiply(B,A);

Newton,
NumberOfSignificantDigits, PolynomialInterpolation, Quadrature,
RateOfConvergence,
RelativeError, RemainderTerm, Roots, RungeKutta, Secant, SpectralRadius,
Steffensen, Taylor,
TaylorPolynomial, UpperBoundOfRemainderTerm, VectorLimit]
Gói lệnh Student hỗ trợ cho việc dạy và học toán.
• Từ Maple 8, gói lệnh Student được phát triển từ gói lệnh student trước đó nhằm hỗ trợ
cho việc dạy và học toán ở đại học và phổ thông. Khai thác khả năng của gói lệnh này sẽ
đem đến cho giáo viên rất nhiều công cụ hỗ trợ mới trong phương pháp dạy học. Có thể
nói rằng gói lệnh này đã đề cập đến tất cả các nội dung toán học của đại học và phổ
6


thông, cung cấp nhiều lệnh và thủ tục cho các phép toán và algorithm xuất hiện trong
chương trình giảng dạy, cung cấp nhiều công cụ tương tác dưới dạng Maplet và hỗ trợ
việc làm
từng bước các phép toán cơ bản của vi tích phân.
• Gói lệnh Student có 3 gói lệnh con là
•

Calculus1

•

LinearAlgebra

•

Precalculus


> Norm(A);
> Norm(A,infinity);
Tính giá trị lớn nhất của chuẩn ma trận :
> for i from 1 to 4 do
R[i]:=Row(A,i);
r[i]:=add(abs(R[i][j]),j=1..4);
od;
supnorm:=max(r[1],r[2],r[3],r[4]);

8


> Norm(A,2);

> evalf(%);
> Norm(A,Euclidean);

> evalf(%);

III.Giải hệ phương trình tuyến tính bằng kỹ thuật lặp:
1.Phép lặp Jacobi :
Để giải quyết bài toán Ax = b với xấp xỉ gần đúng ban đầu x ( 0)
- Đầu vào : số phương trình và số nghiệm là n ; ma trận A với các thành phần
aij , 1 ≤ i, j ≤ n ;vecto b với thành phần b j ,1 ≤ j ≤ n ;các thành phần
XOi ,1 ≤ i ≤ n của vecto XO = x (0 ) ;sai số cho phép TOL;số lần lặp lớn nhất N
- Đầu ra : các nghiệm xấp xỉ x1 , x2 ,..., xn hoặc thông báo “Số lượng tối đa của
phép lặp đã bị vượt quá”.
Thuật toán :
Bước 1 : Đặt k = 1


XO = x (0 ) , TOL, N
k=1

k≤N

F
Output(‘Max…’)
STOP

T
Tính các xi

x − XO < TOL

T

Output(‘ x1 ,..., xn ’)
STOP

F
k = k +1
Tính các XOi = xi

> restart;
> with(LinearAlgebra):
> libname:="c:/nalib",libname;
Chương trình con cho giải thuật Jacobi
>jacobi:=proc(A::Matrix,b::Vector,x0::Vector,tol::positive,
n::nonnegint,v::name)

> S := S-A[I,J]*X0[J];
> od;
> for J from I+1 to N do
> S := S-A[I,J]*X0[J];
od;
> S := evalf((S+B[I])/A[I,I]);
> if abs(S-X0[I]) > ERR then
> ERR := abs(S-X0[I]);
> fi;
> X1[I] := S;
> od;
printf(`%3d
[`,K);
for I from 1 to N do
printf(`% 12.8f`,X1[I]);
if IN then printf(`,`) fi;
od;
printf(`]\n`);
> if ERR OK := TRUE;
> fi;
> K := K+1;
> for I from 1 to N do
> X0[I] := X1[I];
> od;
> od;
> if OK = FALSE then
> printf(`Maximum Number of Iterations Exceeded.\n`);
> else
> printf(`\nThe solution vector is\n`,args[6]);

6
[ 1.00319865, 1.99224126, -0.99452174, 0.99443374]
7
[ 0.99812847, 2.00230688, -1.00197223, 1.00359431]
8
[ 1.00062513, 1.99867030, -0.99903558, 0.99888839]
9
[ 0.99967415, 2.00044767, -1.00036916, 1.00061919]
10
[ 1.00011860, 1.99976795, -0.99982814, 0.99978598]
11
[ 0.99994242, 2.00008478, -1.00006833, 1.00010850]
12
[ 1.00002214, 1.99995896, -0.99996916, 0.99995967]
13
[ 0.99998973, 2.00001582, -1.00001256, 1.00001924]
14
[ 1.00000410, 1.99999268, -0.99999444, 0.99999250]
15
[ 0.99999816, 2.00000292, -1.00000230, 1.00000344]
16
[ 1.00000075, 1.99999868, -0.99999899, 0.99999862]
17
[ 0.99999967, 2.00000054, -1.00000042, 1.00000062]
18
[ 1.00000014, 1.99999976, -0.99999982, 0.99999975]
The solution vector is

> v:='v';


0.99888839]
1.00061919]
0.99978598]

2.Thuật toán Lặp Gauss-Seidel :
Giải bài toán Ax = b với xấp xỉ gần đúng là x (0 )
- Đầu vào : Số lượng phương trình và nghiệm n ;các phần tử aij ,1 ≤ i, j ≤ n của
ma trận A;các thành phần bi ,1 ≤ i ≤ n của vecto b ;các thành phần
XOi ;1 ≤ i ≤ n của vecto XO = x (0 ) ;sai số cho phép TOL;Số lần lặp lớn nhất N.

Đầu ra : các nghiệm xấp xỉ x1 , x2 ,..., xn hoặc thông báo “Số lượng tối đa của
phép lặp đã bị vượt quá”.
Thuật toán :
Bước 1 : Đặt k = 1
Bước 2 : Trong khi (k ≤ N ) thì thực hiện từ Bước 3 đến Bước 6
Bước 3 : Cho i chạy từ 1 đến n , ta tính các giá trị
-

i −1

xi =

− ∑ aij x (j0 ) −
j =1

n

∑a

j = i +1

T
Tính các xi

x − XO < TOL

T

Output(‘ x1 ,..., xn ’)
STOP

F
k = k +1
Tính các XOi = xi

> restart;
> with(LinearAlgebra):
> libname:="c:/nalib",libname;

Giải thuật Gauss-Seidel :
>gaussseidel:=roc(A::Matrix,b::Vector,x0::Vector,tol::
positive,n::nonnegint,v::name)
local AA, B, OK, N, I, J, X0, X1, TOL, NN, K, ERR, S,
X;
with(LinearAlgebra);
N:=RowDimension(A);
AA:=A;
B:=b;
X0:=x0;
TOL:=tol;
NN:=n;

> S := evalf((S+B[I])/A[I,I]);
> if abs(S-X0[I]) > ERR then
> ERR := abs(S-X0[I]);
> fi;
> X1[I] := S;
> od;
printf(`%3d
[`,K);
for I from 1 to N do
printf(`% 12.8f`,X1[I]);
if IN then printf(`,`) fi;
od;
printf(`]\n`);
STEP 4
> if ERR OK := TRUE;
> fi;
STEP 5
> K := K+1;
STEP 6
> for I from 1 to N do
> X0[I] := X1[I];
> od;
> od;
> if OK = FALSE then
> printf(`Maximum Number of Iterations Exceeded.\n`);
STEP 7
> else
> printf(`\nThe solution vector is\n`,args[6]);
> v:=evalm(X1);

[ 1.00000836, 2.00000117, -1.00000274, 0.99999922]
7
[ 1.00000067, 2.00000002, -1.00000021, 0.99999996]
8
[ 1.00000004, 2.00000000, -1.00000001, 1.00000000]
The solution vector is

Từ đây ta thấy giải thuật Gauss-Seidel có tốc độ hội tụ nhanh hơn giải thuật
Jacobi.

3.Thuật toán SOR :
Giải bài toán Ax = b với tham số ω và xấp xỉ gần đúng ban đầu x (0 )
- Đầu vào : Số phương trình và số ẩn n ; các phần tử aij ,1 ≤ i, j ≤ n của ma trận
A;các thành phần bi ,1 ≤ i ≤ n của vecto b ;các thành phần XOi ;1 ≤ i ≤ n của vecto

XO = x (0 ) ; tham số ω ; sai số cho phép TOL ; số lần lặp lớn nhất N.
- Đầu ra : các nghiệm xấp xỉ x1 , x2 ,..., xn hoặc thông báo “Số lượng tối đa của
phép lặp đã bị vượt quá”.
Thuật toán :
16


Bước 1 : Đặt k = 1
Bước 2 : Trong khi (k ≤ N ) thì thực hiện lặp từ Bước 3 đến Bước 6
Bước 3 : Tính x = ( xi ) với i = 1,…, n

⎛

i −1


Thuật toán SOR :
> SOR :=
proc(A::Matrix,b::Vector,w::numeric,x0::Vector,tol::po
sitive,n::nonnegint,v::name)
local AA, B, OK, N, I, J, X0, X1, TOL, NN, K, ERR, S,
X, W;
with(LinearAlgebra);
N:=RowDimension(A);
AA:=A;
B:=b;
X0:=x0;
TOL:=tol;
NN:=n;
X1:=Vector(N);
W:=w;
> K := 0;
> OK := FALSE;
printf(` n
x\n`);
printf(` -\n`);
printf(`%3d
[`,K);
for I from 1 to N do
printf(`% 12.8f`,X0[I]);
if IN then printf(`,`) fi;
od;
printf(`]\n`);
K:=1;
> while OK = FALSE and K ERR := 0;

> od;
> od;
> if OK = FALSE then
> printf(`Maximum Number of Iterations Exceeded.\n`);
> else
> printf(`\nThe solution vector is\n`);
> v:=evalm(X1);
fi;
end;
Nhập dữ liệu và tính toán :
> A:=Matrix(3,3,[[4,3,0],[3,4,-1],[0,-1,4]]);

> b:=Vector(3,[24,30,-24]);

> v:='v';
> > v:=SOR(A,b,1,Vector(3,[1,1,1]),10^(-6),100,v);
n
x
0
[ 1.00000000, 1.00000000, 1.00000000]
1
[ 5.25000000, 3.81250000, -5.04687500]
2
[ 3.14062500, 3.88281250, -5.02929688]
18


3
4
5

[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[

3.08789062,
3.05493164,
3.03433228,
3.02145767,
3.01341104,
3.00838190,
3.00523869,
3.00327418,
3.00204636,
3.00127898,
3.00079936,

3.99993647,
3.99996029,
3.99997518,
3.99998449,
3.99999031,
3.99999394,
3.99999622,
3.99999764,
3.99999852,
3.99999908,

-5.01831055]
-5.01144409]
-5.00715256]
-5.00447035]
-5.00279397]
-5.00174623]
-5.00109139]
-5.00068212]
-5.00042632]
-5.00026646]
-5.00016654]
-5.00010408]
-5.00006505]
-5.00004066]
-5.00002541]
-5.00001588]
-5.00000993]
-5.00000620]
-5.00000388]

7
[ 3.00004980, 4.00025858, -5.00034865]
8
[ 2.99974513, 4.00006534, -4.99989242]
9
[ 3.00000246, 4.00001498, -5.00002222]
10
[ 2.99998534, 4.00000306, -4.99999349]
11
[ 3.00000080, 4.00000052, -5.00000146]
12
[ 2.99999931, 4.00000006, -4.99999961]
13
[ 3.00000012, 4.00000000, -5.00000010]
The solution vector is
19


Thuật toán SOR ngắn gọn hơn nhờ vào tham số Omega .. tốc độ làm việc và hội
tụ của thuật toán SOR có thể nói là nhanh hơn khi ta tăng Omega.
n,

A(aij ), b(b j ),ϖ

,

(0 )

XO = x , TOL, N
k=1

t
Đầu ra : xấp xỉ xx = ( xxi ,…, xxn ) hoặc thông báo về số lần lặp bị vượt quá và

một xấp xỉ COND gần đến K ∞ ( A) .
Thuật toán :
Bước 0 : Giải hệ Ax = b tìm nghiệm x1 ,…, xn bằng phép loại trừ Gauss , bằng cách
giữ lại các hệ số nhân mij , j = i + 1, i + 2,…, n, i = 1,2,…, n − 1 và
Bước 1 : Đặt k = 1
Bước 2 : Trong khi (k ≤ N ) thì thực hiện từ Bước 3 đến Bước 9
Bước 3 : Tính vecto r = (ri )
n

ri = bi − ∑ aij x j
j =1

Bước 4 : Giải hệ tuyến tính Ay = r sử dụng phép loại trừ Gauss như ở Bước 0.
Bước 5 : Tính các xxi với i = 1,…, n

xxi = xi + yi
20


Bước 6 : Nếu k = 1 thì đặt COND =

Bước 7 : Nếu x − xx

∞

y



k=1

F

k≤N

T
v ≤ TOL

T

Output(‘

x1 ,.., xn ; r1 ,.., rn ’)
STOP

F
u = Av = Av (k )
t = tk =

α

; x = x (k ) = x (k −1) + tv (k )

n

∑v u
j =1


α =β

k = k +1

22


5.Thuật toán Gradient điều kiện kết hợp :
Giải bài toán Ax=b với ma trận điều kiện C −1 và xấp xỉ gần đúng x (0 )
o Đầu vào : số lượng phương trình là số chưa biết n ; ma trận A với các
thành phần aij , 1 ≤ i, j ≤ n ; vecto b với thành phần b j ,1 ≤ j ≤ n ; ma
trận hệ số điều kiện C −1 với các thành phần là γ ij ,1 ≤ i, j ≤ n ; vecto xấp
xỉ đầu x (0 ) với các xi ,1 ≤ i ≤ n và số lần lặp lớn nhất là N,sai số cho
phép TOL
o Đầu ra : các nghiệm xấp xỉ x1 , x2 ,..., xn và phần dư r1 , r2 ,..., rn hoặc thông
báo “Số lượng tối đa của phép lặp đã bị vượt quá”.
Thuật toán :
Bước 1 : Tính
r (0 ) = b − Ax (0 )
w = w(0 ) = C −1r
v = v (1) = C −t w
n

α = ∑ w2j
j =1

Bước 2 : Đặt k = 1
Bước 3 : Kiểm tra điều kiện
Trong khi k ≤ N thì làm từ Bước 4 đến Bước 7
Bước4 : If v ≤ TOL then

printf (Nghiệm là : x1 , x2 ,..., xn )
printf (với phần dư là : r1 , r2 ,..., rn )
STOP
Else (Sang bước 7)
Bước 7 : Tính

s = sk =

β
α

v = v (k + ) = C −1w + sv (k )
23


α = β ; (cập nhật α )
k = k +1
Quay lại Bước 3
Bước 8 : Kiểm tra điều kiện
If (k〉 n ) then
Printf(Số lượng tối đa của phép lặp đã bị vượt quá ! )
STOP.

24


n,

A(aij ), b(b j ) , C −1



; x = x (k ) = x (k −1) + tv (k )

n

∑v u
j =1

j

j

r = r (k ) = r ( k −1) − tu

w = w(k ) = C −1r (k )
n

β = ∑ w2j
j =1

β < TOL

T

x1 ,.., xn ; r1 ,.., rn ’)
STOP

F

s = sk =