Giáo trình Phương pháp tính - Chương 3: Hệ phương trình đại số tuyến tính
Chöông 3
HEÄ PHÖÔNG TRÌNH ÑAÏI SOÁ
TUYEÁN TÍNH
Trong chöông naøy chuùng ta neâu leân moät soá phöông phaùp duøng ñeå
giaûi heä phöông trình ñaïi soá tuyeán tính
Ax = b,
(3.1)
raát thöôøng gaëp trong caùc baøi toaùn khoa hoïc kó thuaät. Ta chæ xeùt heä
goàm n phöông trình vôùi n aån. Do vaäy ma traän heä soá A laø ma traän
vuoâng caáp n, vaø vectô nghieäm x cuõng nhö vectô töï do b laø caùc vectô
coät n chieàu thuoäc Rn. Ta luoân giaû thieát raèng det A = 0, vaø do ñoù bao
giôø heä cuõng coù nghieäm duy nhaát x = A−1b. Tuy nhieân vieäc tìm ma
traän nghòch ñaûo A−1 ñoâi khi coøn khoù khaên gaáp nhieàu laàn so vôùi vieäc
giaûi tröïc tieáp heä phöông trình xuaát phaùt. Döôùi ñaây chuùng ta seõ xeùt
moät soá phöông phaùp thöôøng duøng ñeå giaûi heä phöông trình (3.1).
3.1 PHÖÔNG PHAÙP GAUSS
Tröôùc khi trình baøy phöông phaùp Gauss, chuùng ta xeùt moät soá
tröôøng hôïp ñôn giaûn khi ma traän heä soá A cuûa heä phöông trình (3.1)
coù daïng ñaët bieät.
36
HEÄ PHÖÔNG TRÌNH ÑAÏI SOÁ TUYEÁN TÍNH
Tröôøng hôïp ñôn giaûn nhaát laø tröôøng hôïp heä phöông trình coù ma
traän heä soá coù daïng ñöôøng cheùo:
a11
0
0
· · ·
0
0
a22 · · ·
A =
· · · · · · · · · · · ·
0
0
· · · ann
Khi aáy heä töông ñöông vôùi n phöông trình baäc nhaát aiixi = bi, ∀i = 1, n.
Vì det A = a11a22 · · ·ann = 0 neân aii = 0, ∀i. Vaø do ñoù nghieäm cuûa heä
coù theå ñöôïc vieát döôùi daïng:
bi
xi =
,
i = 1, 2, 3, · ·· , n
aii
Tröôøng hôïp thöù hai khi ma traän heä soá A coù daïng tam giaùc treân:
a11 a12 · · · a1n
0
a22 · · · a2n
· · · · · · · · · · · ·
· · · ann
A =
0
0
Vôùi giaû thieát det A = 0, ta coù aii = 0, ∀i = 1, n, vaø nghieäm cuûa heä ñöôïc
cho bôûi coâng thöùc:
bn
xn =
ann
n
(3.2)
X
1
akk
xk =
bk −
akjxj
,
k = n − 1, · · · , 1
j=k+1
Cuoái cuøng khi ma traän heä soá A coù daïng tam giaùc döôùi:
a11
0
· · ·
0
0
a21 a22 · · ·
A =
· · · · · · · · · · · ·
an1 an2 · · · ann
3.1 Phöông phaùp Gauss
37
Töông töï det A = 0 ⇒ aii = 0, ∀i = 1, n, vaø nghieäm cuûa heä coù daïng:
b1
x1 =
a11
k−1
(3.3)
X
1
akk
xk =
bk −
akjxj
,
k = 2, · · · , n
j=1
Thuaät toaùn giaûi heä phöông trình vôùi ma traän tam giaùc ñöôïc theå
hieän trong Chöông trình 3.1 vaø 3.2. Ñoái soá cuûa chöông trình goàm: N
laø soá phöông trình vaø soá aån, a laø ma traän heä soá caáp N × (N + 1), coät
thöù N + 1 laø vectô töï do. Keát quaû traû veà cuûa chöông trình laø vectô
nghieäm x .
Chöông trình 3.1. - c3upper : Ma traän heä soá tam giaùc treân.
function [x] = c3upper(N,a)
if nargin < 2, error('Haøm coù toái thieåu 2 ñoái soá');end;
x(N)=a(N,N+1)/a(N,N);
for k=N-1:-1:1
sum = 0;
for j=k+1:N
sum=sum+a(k,j)*x(j);
end;
x(k)=(a(k,N+1)-sum)/a(k,k);
end;
Chöông trình 3.2. - c3lower : Ma traän heä soá tam giaùc döôùi.
function [x] = c3lower(N,a)
if nargin < 2, error('Haøm coù toái thieåu 2 ñoái soá');end;
x(1)=a(1,N+1)/a(1,1);
for k=2:N
sum = 0;
38
HEÄ PHÖÔNG TRÌNH ÑAÏI SOÁ TUYEÁN TÍNH
for j=1:k-1
sum=sum+a(k,j)*x(j);
end;
x(k)=(a(k,N+1)-sum)/a(k,k);
end;
Baây giôø chuùng ta seõ trình baøy phöông phaùp Gauss ñeå giaûi heä
phöông trình toång quaùt daïng (3.1). Noäi dung cuûa phöông phaùp Gauss
duøng ñeå giaûi heä phöông trình ñaïi soá tuyeán tính laø söû duïng caùc pheùp
bieán ñoåi sô caáp theo haøng ñeå chuyeån veà moät heä phöông trình môùi
töông ñöông vôùi heä phöông trình cuõ maø ma traän heä soá coù daïng tam
giaùc. Caùc pheùp bieán ñoåi sô caáp thöôøng hay söû duïng laø:
• Nhaân moät haøng cho moät soá khaùc khoâng.
• Hoaùn chuyeån hai haøng cho nhau.
• Coäng moät haøng cho moät haøng khaùc ñaõ nhaân vôùi moät soá khaùc
khoâng.
Xeùt heä thoáng phöông trình sau:
a11x1
a21x1
. . .
+
+
a12x2
a22x2
. . .
+
+
. . .
. . .
. . .
. . .
+
+
a1nxn
a2nxn
=
=
b1
b2
. . .
bn
an1x1
+
an2x2
+
+
annxn
=
Do ñònh thöùc cuûa ma traän heä soá A khaùc khoâng neân moät trong caùc
soá a11, a12, . . ., a1n phaûi khaùc khoâng. Giaû söû a11 = 0. Laáy phöông trình
ak1
thöù k vôùi k = 2, n tröø cho phöông trình moät ñaõ nhaân vôùi
, ta ñöôïc
a11
moät heä môùi coù daïng nhö sau:
a11x1
+
a12x2
+
+
. . .
+
+
a1nxn
=
=
b1
(1)
(1)
a
x2
. . .
. . .
a
xn
b(21)
. . .
22
2n
. . .
(1)
n2
(1)
a
x2
+
. . .
+
a
x
=
b(n1)
nn
n
3.1 Phöông phaùp Gauss
39
(1)
22
(1)
n2
Trong caùc soá a , . . ., a phaûi coù moät soá khaùc khoâng, vì neáu ngöôïc
(1)
22
laïi thì det A = 0, traùi vôùi giaû thieát. Giaû söû a = 0. Coøn neáu chæ coù
(1)
p2
(1)
22
a
= 0 vaø a = 0 thì ta thöïc hieän pheùp hoaùn chuyeån hai phöông
trình thöù 2 vaø thöù p. Tieáp tuïc bieán ñoåi cho n − 2 phöông trình cuoái.
Vaø cöù tieáp tuïc cho ñeán phöông trình thöù n, ta ñöôïc heä phöông trình
sau
a11x1
+
a12x2
(1)
+
+
. . .
+
+
a1nxn
(1)
2n
=
=
b1
a
x2
22
. . .
. . .
. . .
a
xn
b(21)
. . .
(n−1)
a
xn
=
b(nn−1)
nn
Ñaây laø heä phöông trình coù ma traän heä soá coù daïng tam giaùc treân
vaø coù theå giaûi ñöôïc baèng coâng thöùc (3.2).
Ví duï 3.1. Xeùt heä phöông trình ñaïi soá tuyeán tính sau:
x1
2x1
x1
−
−
+
−
x2
2x2
x2
+
+
+
+
2x3
3x3
x3
−
−
x4
3x4
=
=
=
=
−8
−20
−2
4
x1
x2
4x3
+
3x4
Ma traän heä soá môû roäng coù daïng
1
2
1
1
−1
−2
1
−1
2
3
1
4
−1
−3
0
:
−8
:
:
:
−20
−2
4
A(0)
=
3
Ta thöïc hieän caùc pheùp bieán ñoåi sau: (h2 = h2 − 2h1), (h3
=
h3 − h1), (h4 = h4 − h1), khi ñoù ma traän trôû thaønh
1
0
0
0
−1
0
2
2
−1
:
:
:
:
−8
−4
6
−1 −1
A(1)
=
−1
1
4
0
2
12
(1)
22
Phaàn töû a
= 0, do ñoù ñeå tieáp tuïc, ta thöïc hieän pheùp chuyeån ñoåi
40
HEÄ PHÖÔNG TRÌNH ÑAÏI SOÁ TUYEÁN TÍNH
giöõa haøng thöù hai vaø thöù ba vaø thu ñöôïc
1
0
0
0
−1
2
0
2
−1
:
:
:
:
−8
6
−1
1
A(2)
=
−1 −1
−4
0
2
4
12
Cuoái cuøng laáy haøng thöù tö coäng cho hai laàn haøng thöù ba ta ñöôïc:
1
0
0
0
−1
2
0
2
−1
:
:
:
:
−8
6
−1
1
A(3)
=
−1 −1
−4
0
0
2
4
Vaø söû duïng coâng thöùc (3.2) ta coù theå deã daøng tìm ñöôïc x =
[−7, 3, 2, 2]T.
Thuaät toaùn giaûi heä phöông trình baèng phöông phaùp Gauss ñöôïc
theå hieän trong Chöông trình 3.3. Ñoái soá cuûa chöông trình goàm: N laø
soá phöông trình vaø soá aån, a laø ma traän heä soá caáp N × (N + 1), coät
thöù N + 1 laø vectô töï do. Keát quaû traû veà cuûa chöông trình laø vectô
nghieäm x .
Chöông trình 3.3. - c3gauss : Phöông phaùp Gauss.
function [x] = c3gauss(N,a)
if nargin < 2, error('Haøm coù toái thieåu 2 ñoái soá');end
for k=1:N
if a(k,k)==0
flag=0;
for i=k+1:N
if a(i,k)∼=0
flag=1;
for j=1:N+1
tmp=a(k,j);
a(k,j)=a(i,j);
a(i,j)=tmp;
3.1 Phöông phaùp Gauss
41
end;
break;
end;
end;
if flag==0
error('Ma traän suy bieán.');
end;
end;
for i=k+1:N
tmp=a(i,k);
for j=k:N+1
a(i,j)=a(i,j)-tmp*a(k,j)/a(k,k);
end;
end;
end;
x=c3upper(N,a);
(1)
Trong ví duï 3.1 ôû phaàn treân, ôû böôùc thöù hai, do a = 0 neân ta
22
phaûi hoaùn chuyeån hai haøng thöù hai vaø thöù ba. Ñeå traùnh tröôøng hôïp
naøy, ta coù theå caûi tieán phöông phaùp Gauss theo höôùng nhö sau. Taïi
moãi böôùc, khi choïn phaàn töû ñeå bieán ñoåi, ta seõ choïn phaàn töû coù trò
tuyeät ñoái lôùn nhaát, sao cho khoâng cuøng haøng vaø coät vôùi nhöõng phaàn
töû ñaõ choïn tröôùc. Phaàn töû nhö vaäy thöôøng ñöôïc goïi laø phaàn töû chính
hay phaàn töû troäi. Sau ñoù ta seõ bieán ñoåi ñeå cho taát caû caùc phaàn töû treân
cuøng coät cuûa phaàn töû troäi baèng khoâng. Qua n böôùc nhö vaäy ta seõ tìm
1
ñöôïc nghieäm deã daøng . Ta minh hoaï phöông phaùp naøy baèng ví duï
sau.
Ví duï 3.2. Xeùt heä phöông trình trong ví duï tröôùc coù ma traän heä soá môû roäng
1
2
1
1
−1
−2
1
−1
2
3
1
4
−1
−3
0
:
:
:
:
−8
−20
−2
4
A(0)
=
3
1Phöông phaùp naøy cuõng ñöôïc goïi laø phöông phaùp Gauss-Jordan hay Jordan
42
HEÄ PHÖÔNG TRÌNH ÑAÏI SOÁ TUYEÁN TÍNH
(0)
43
Ñaàu tieân ta seõ choïn phaàn töû chính laø phaàn töû a
= 4 vaø thöïc hieän
caùc pheùp bieán ñoåi (4h3 − h4), (4h2 − 3h4), (2h1 − h4) ta thu ñöôïc
1
5
3
1
−1
−5
5
−1
0
0
0
4
−5
−21
−3
3
:
:
:
:
−20
−92
−12
4
A(1)
=
Böôùc tieáp theo, phaàn töû chính ñöôïc choïn khoâng ñöôïc naèm treân haøng
(1)
24
thöù tö vaø coät thöùc ba. Ñoù laø phaàn töû a = −21. Tieáp tuïc thöïc hieän
caùc pheùp bieán ñoåi (21h1 − 5h2), (7h3 − h2), (7h4 + h2) ta thu ñöôïc
−4
4
−5
40
0
0
0
0
−21
0
:
:
:
:
40
−92
8
−64
5
A(2)
=
16
12 −12 28
0
Tieáp theo phaàn töû chính ñöôïc choïn khoâng ñöôïc naèm treân haøng thöù
hai, thöù tö vaø coät thöù ba, thöù tö vaø do ñoù phaàn töû chính seõ laø phaàn töû
(2)
32
a
= 40. Thöïc hieän pheùp bieán ñoåi (10h1 −h3), (8h2 +h3), (10h4 +
3h3) ta ñöôïc
−56
0
0
0
0
0
0
−168
0
:
:
:
:
392
−728
8
−616
56
A(3)
=
16 40
168
0
280
0
Cuoái cuøng phaàn töû chính khoâng cuøng naèm treân haøng vaø coät cuûa
nhöõng phaàn töû chính ñaõ ñöôïc choïn tröôùc laø phaàn töû a
Thöïc hieän caùc pheùp bieán ñoåi (h2 + h1), (7h3 + 2h1), (h4 + 3h1) ta
(3)
11
= −56.
coù ma traän cuoái cuøng
−56
0
0
280
0
0
0
0
0
−168
0
:
:
:
:
392
−336
840
0
0
0
A(3)
=
280
0
560
vaø heä phöông trình ñaàu töông ñöông vôùi heä sau
−56x1
−168x4
280x2
=
=
=
=
392
−336
840
280x3
560
3.1 Phöông phaùp Gauss
43
Töø ñaây chuùng ta cuõng suy ra ñöôïc x = [−7, 3, 2, 2]T.
Thuaät toaùn giaûi heä phöông trình baèng phöông phaùp phaàn töû troäi
ñöôïc theå hieän trong Chöông trình 3.4. Ñoái soá cuûa chöông trình goàm:
N laø soá phöông trình vaø soá aån, a laø ma traän heä soá caáp N × (N + 1),
coät thöù N + 1 laø vectô töï do. Keát quaû traû veà cuûa chöông trình laø vectô
nghieäm x .
Chöông trình 3.4. - c3jordan : Phöông phaùp phaàn töû troäi.
function [x] = c3jordan(N,a)
if nargin < 2, error('Haøm coù toái thieåu 2 ñoái soá'); end
for i=1:N, b(i)=0; end;
for k=1:N
max=0;
for i=1:N
if b(i)==0
for j=1:N
if max<abs(a(i,j))
max=abs(a(i,j)); im=i; jm=j;
end;
end;
end;
end;
if max<=0, error('Ma traän suy bieán.'); end;
b(im)=jm; max=a(im,jm);
for j=1:N+1, a(im,j)=a(im,j)/max; end;
for i=1:N
if i∼=im
max=a(i,jm);
for j=1:N+1
a(i,j)=a(i,j)-a(im,j)*max;
end;
end;
end;
44
HEÄ PHÖÔNG TRÌNH ÑAÏI SOÁ TUYEÁN TÍNH
end;
for i=1:N, x(b(i))=a(i,N+1); end;
Caùc phöông phaùp coù söû duïng caùc pheùp bieán ñoåi sô caáp cô baûn coù
öu ñieåm laø ñôn giaûn, deã laäp trình. Tuy nhieân neáu phaàn töû ñöôïc choïn
ñeå bieán ñoåi gaàn vôùi khoâng thì phöông phaùp Gauss coù theå cho keát
quaû khoâng chính xaùc. Hôn nöõa, neáu caùc pheùp toaùn coäng, tröø, nhaân,
chia ñöôïc laøm ñuùng hoaøn toaøn, thì caùc phöông phaùp treân cho chuùng
ta nghieäm ñuùng cuûa heä phöông trình. Tuy nhieân, khi thöïc hieän treân
caùc coâng cuï tính toaùn, ta vaãn gaëp phaûi sai soá laøm troøn. Cho neân caùc
phöông phaùp Gauss vaãn ñöôïc xem nhö laø caùc phöông phaùp gaàn ñuùng.
3.2 PHÖÔNG PHAÙP NHAÂN TÖÛ LU
Noäi dung cuûa phöông phaùp nhaân töû LU laø phaân tích ma traän heä
soá A thaønh tích cuûa hai ma traän L vaø U, trong ñoù L laø ma traän tam
giaùc döôùi vaø U laø ma traän tam giaùc treân. Khi ñoù vieäc giaûi heä phöông
trình (3.1) seõ ñöa veà vieäc giaûi hai heä phöông trình Ly = b vaø Ux = y
maø ma traän heä soá laø caùc ma traän tam giaùc vaø nghieäm thu ñöôïc töø caùc
coâng thöùc (3.2) vaø (3.3). Ta coù ñònh lí sau ñaây.
Ñònh lí 3.1. Neáu A laø ma traän khoâng suy bieán, thì bao giôø cuõng toàn
taïi moät ma traän P khoâng suy bieán sao cho ma traän PA phaân tích
ñöôïc thaønh tích cuûa ma traän tam giaùc döôùi L vaø ma traän tam giaùc
treân U, nghóa laø PA = LU.
Coù raát nhieàu phöông phaùp phaân tích A = LU, tuy nhieân ta thöôøng
xeùt tröôøng hôïp ma traän L coù ñöôøng cheùo chính baèng 1 vaø goïi laø
phöông phaùp Doolittle. Khi ñoù L vaø U coù daïng:
1
l21
0
1
· · ·
· · ·
0
0
u11 u12 · · · u1n
0
u22 · · · u2n
· · · · · · · · · · · ·
· · · unn
L =
vaø U =
· · · · · · · · · · · ·
ln1 ln2 · · ·
1
0
0
3.2 Phöông phaùp nhaân töû LU
45
Caùc phaàn töû cuûa hai ma traän L vaø U ñöôïc xaùc ñònh theo coâng thöùc
sau:
u1j = a1j
(1 6 j 6 n)
(2 6 i 6 n)
ai1
li1
=
u11
i−1
P
aij
uij = aij
−
likukj
(1 < i 6 j)
(1 < j < i)
(3.4)
k=1
!
j−1
X
1
ujj
lij
=
−
likukj
k=1
Ví duï 3.3. Xeùt heä phöông trình
2x1
−4x1
2x1
+
2x2
3x2
x2
−
3x3
4x3
2x3
=
=
=
9
−15
3
−
+
+
+
Ta phaân tích ma traän heä soá
2
2
−3
4
2
1
l21
l31 l32
0
1
0
0
1
u11 u12 u13
−4 −3
0
0
u22 u23
=
·
2
1
0
u33
vôùi
Do ñoù
u11 = 2; u12 = 2; u13 = −3; l21 = −2; l31 = 1;
u22 = 1; u23 = −2;l32 = −1; u33 = 3.
1
−2
1
2
0
0
1
0
0
1
9
−15
3
9
3
Ly = b ⇐⇒
y =
=⇒ y =
−1
−3
2
1
0
−3
9
3
2
1
Ux = y ⇐⇒
−2
x =
=⇒ x =
0
3
−3
−1
Ví duï 3.4. Heä phöông trình
x1
+
+
−
+
x2
x2
x2
2x2
+
+
+
+
−
3x4
x4
2x4
x4
=
=
=
=
4
1
2x1
3x1
−x1
−
−
+
x3
x3
3x3
−3
4
46
HEÄ PHÖÔNG TRÌNH ÑAÏI SOÁ TUYEÁN TÍNH
Ta coù phaân tích cuûa ma traän heä soá A = LU theo phöông phaùp
Doolittle nhö sau:
1
2
3
0
1
4
0
0
1
0
0
0
0
1
1
0
0
0
1
0
3
−1 −1 −5
A =
0
0
3
0
13
−13
−1 −3
Ta coù
1
2
3
0
1
4
0
0
1
0
0
0
0
1
y1
y2
y3
y4
4
1
Ly = b ⇐⇒
=
−3
−1 −3
4
Ta thu ñöôïc y = [4, −7, 13, −13]T. Vaø cuoái cuøng töø heä phöông trình
1
0
0
0
1
0
3
x1
x2
x3
x4
4
−7
−1 −1 −5
Ux = y ⇐⇒
=
0
0
3
0
13
13
−13
−13
ta coù nghieäm x = [−1, 2, 0, 1]T.
Thuaät toaùn phaân raõ LU ñöôïc theå hieän trong Chöông trình 3.5.
Ñoái soá cuûa chöông trình goàm: N laø soá phöông trình, a laø ma traän heä
soá caáp N × N. Keát quaû traû veà cuûa chöông trình laø ma traän tam giaùc
döôùi l vaø ma traän tam giaùc treân u .
Chöông trình 3.5. - c3LUfactor : Phöông phaùp nhaân töû LU.
function [l,u] = c3LUfactor(N,a)
if nargin < 2, error('Haøm coù toái thieåu 2 ñoái soá'); end;
l=zeros(N); u=zeros(N);
for i=1:N, l(i,i)=1; end;
for j=1:N, u(1,j)=a(1,j); end;
if u(1,1)==0, error('Khoâng theå phaân tích ñöôïc.');end;
for i=2:N, l(i,1)=a(i,1)/u(1,1); end;
for i=2:N-1
3.2 Phöông phaùp nhaân töû LU
47
for j=i:N
sum=0;
for k=1:i-1,sum=sum+l(i,k)*u(k,j);end;
u(i,j)=a(i,j)-sum;
end;
if
u(i,i)==0,error('Khoâng
theå
phaân
tích
ñöôïc.');end;
for j=i+1:N
sum=0;
for k=1:i-1,sum=sum+l(j,k)*u(k,i);end;
l(j,i)=(a(j,i)-sum)/u(i,i);
end;
end;
sum=0; for k=1:N-1,sum=sum+l(N,k)*u(k,N); end;
u(N,N)=a(N,N)-sum;
Phöông phaùp phaân raõ LU aùp duïng raát hieäu quaû trong tröôøng hôïp
ma traän heä soá coù daïng ba ñöôøng cheùo:
a11 a12
0
· · ·
0
0
0
· · ·
0
0
0
· · ·
a21 a22 a23 · · ·
0
a32 a33 · · ·
A =
· · · · · · · · · · · ·
0
0
0
0
0
0
· · · an−1,n−1 an−1,n
· · ·
an,n−1
an,n
Khi ñoù phaân raõ Doolittle cho ta
1
l21
0
0
1
l32
0
0
1
· · ·
· · ·
· · ·
0
0
0
u11 u12
0
· · ·
u22 u23 · · ·
u33 · · ·
0
0
0
0
0
0
L =
vaø U =
· · · · · · · · · · · · · · ·
· · · · · · · · · · · · · · ·
0
0
0
· · ·
1
0
0
0
· · · un,n
48
HEÄ PHÖÔNG TRÌNH ÑAÏI SOÁ TUYEÁN TÍNH
Khi aáy töø coâng thöùc (3.4) ta coù
a21
u11
u11 = a11; u12 = a12; l21
=
;
ui,i = ai,i − li,i−1ui−1,i
;
=
∀i = 2, 3, . . ., n
ai+1,i
ui,i+1 = ai,i+1; li+1,i
;
∀i = 2, 3, . . ., n − 1
ui,i
Ñeå theå hieän thuaät toaùn Doolittle trong tröôøng hôïp ma traän heä soá
coù daïng ba ñöôøng cheùo, thay vì ñöa vaøo ma traän A, ta ñöa vaøo ba
vectô: a chöùa caùc phaàn töû treân ñöôøng cheùo chính cuûa A, b laø ñöôøng
cheùo naèm treân ñöôøng cheùo chính vaø c laø ñöôøng cheùo naèm döôùi
ñöôøng cheùo chính, vectô töï do laø d . Chöông trình treân MatLab ñöôïc
theå hieän trong chöông trình 3.6. Keát quaû traû veà laø vectô nghieäm x .
Chöông trình 3.6. - c3tridiag : Ma traän ba ñöôøng cheùo.
function [x]=c3tridiag(N,a,b,c,d)
if nargin < 5, error('Ham co toi thieu 5 doi so.'); end;
u(1)=a(1);v(1)=b(1);l(1)=c(1)/u(1);y(1)=d(1);
for i=2:N-1
u(i)=a(i)-l(i-1)*v(i-1);
v(i)=b(i);
l(i)=c(i)/u(i);
y(i)=d(i)-l(i-1)*y(i-1);
end;
u(N)=a(N)-l(N-1)*v(N-1);
y(N)=d(N)-l(N-1)*y(N-1);
x(N)=y(N)/u(N);
for i=N-1:-1:1
x(i)=(y(i)-v(i)*x(i+1))/u(i);
end;
3.3 Phöông phaùp Choleski
49
3.3 PHÖÔNG PHAÙP CHOLESKI
Ñaây laø tröôøng hôïp ñaët bieät cuûa phöông phaùp nhaân töû LU, vaø ñöôïc
duøng cho tröôøng hôïp ma traän heä soá A ñoái xöùng vaø xaùc ñònh döông.
Ma traän vuoâng A ñöôïc goïi laø ñoái xöùng neáu AT = A. Coù theå noùi raèng
ma traän A laø ñoái xöùng neáu caùc phaàn töû cuûa noù ñoái xöùng vôùi nhau qua
ñöôøng cheùo chính, nghóa laø aij = aji, ∀i, j = 1, n. Coøn ma traän A laø
xaùc ñònh döông neáu ∀x ∈ Rn, x = 0 : xT Ax > 0. Ñeå kieåm tra tính xaùc
ñònh döông cuûa moät ma traän, ta thöôøng duøng ñònh lí sau ñaây.
Ñònh lí 3.2. Moät ma traän laø xaùc ñònh döông khi vaø chæ khi taát caû caùc
ñònh thöùc con chính cuûa noù ñeàu döông.
Trong ñoù ñònh thöùc con chính caáp k, 1 6 k 6 n cuûa ma traän laø ñònh
thöùc con thu ñöôïc töø k haøng vaø k coät ñaàu tieân cuûa ma traän ñoù.
1
1
−1
1
2
0
ꢀ
ꢀ
ꢀ
ꢀ
ꢀ
−1
0
4
1
1
Ví duï 3.5. Xeùt ma traän A =
coù ∆1 = |1| = 1 > 0,
ꢀ
ꢀ
ꢀ
ꢀ
ꢀ
ꢀ
ꢀ
ꢀ
ꢀ
ꢀ
ꢀ
ꢀ
ꢀ
ꢀ
ꢀ
1
2
0
−1
0
4
1
1
1
2
∆2 =
= 1 > 0 vaø ∆3 =
= 2 > 0. Do ñoù
−1
A laø xaùc ñònh döông.
Ta coù ñònh lí Choleski:
Ñònh lí 3.3. Ma traän A laø ñoái xöùng vaø xaùc ñònh döông khi vaø chæ khi
toàn taïi moät ma traän B tam giaùc döôùi, khaû ñaûo sao cho A = BBT .
Khi ñoù ma traän tam giaùc döôùi B coù theå tìm ñöôïc theo coâng thöùc
sau:
ai1
√
b11
=
a11; bi1
=
(2 6 i 6 n)
(1 < i 6 n)
b11
s
i−1
P
bii
=
aii −
b2ik
(3.5)
k=1
!
j−1
X
1
bjj
bij
=
aij
−
bikbjk
(1 < j < i)
k=1
50
HEÄ PHÖÔNG TRÌNH ÑAÏI SOÁ TUYEÁN TÍNH
Ví duï 3.6. Xeùt heä phöông trình
1
1
−1
1
2
0
−1
0
4
x1
x2
x3
1
2
3
Ax =
=
= b
1
1
−1
1
2
0
−1
0
4
Vì ma traän A =
laø ñoái xöùng vaø xaùc ñònh döông,
neân töø coâng thöùc (3.5) ta coù theå xaùc ñònh caùc heä soá bij, i > j cuûa
ma traän tam giaùc döôùi B nhö sau
√
b11 = 1, b21 = 1, b31 = −1, b22 = 1, b32 = 1, b33
=
2
1
1
−1
0
1
1
0
0
2
vaø do ñoù A = BBT vôùi B =
. Heä phöông trình
√
ꢁ
By = b
xuaát phaùt seõ töông ñöông vôùi hai heä
. Ta ñöôïc
BT x = y
1
1
1
−1
0
1
1
0
0
2
1
2
3
1
By = b ⇔
y =
x =
⇒ y =
√
3
√
2
1
1
0
0
1
1
0
−1
3
1
T
1
−1/2
B x = y ⇔
⇒ x =
.
√
3
√
3/2
2
2
Thuaät toaùn phaân raõ Choleski ñöôïc theå hieän trong Chöông trình
3.7. Ñoái soá cuûa chöông trình goàm: N laø caáp cuûa ma traän, a laø ma
traän heä soá ñoái xöùng vaø xaùc ñònh döông. Keát quaû traû veà cuûa chöông
trình laø ma traän tam giaùc döôùi b .
Chöông trình 3.7. - c3choleski : Phöông phaùp Choleski.
function [b] = c3choleski(N,a)
if nargin < 2, error('Haøm coù toái thieåu 2 ñoái soá'); end;
for i=1:N
3.3 Phöông phaùp Choleski
51
for j=1:N
if a(i,j) ∼ = a(j,i)
error('Ma traän khoâng ñoái xöùng.');
end;
end;
end;
b=zeros(N);
if a(1,1)<=0, error('Ma traän khoâng xaùc ñònh döông.');
end;
b(1,1)=sqrt(a(1,1));
for i=2:N, b(i,1)=a(i,1)/b(1,1); end;
for k=2:N
ak=0; for j=1:k-1, ak=ak+b(k,j)*b(k,j); end;
ak=a(k,k)-ak;
if ak<=0, error('Ma traän khoâng xaùc ñònh döông.');
end;
b(k,k)=sqrt(ak);
for i=k+1:N
ak=0;
for j=1:k-1
ak=ak+b(i,j)*b(k,j);
end;
b(i,k)=(a(i,k)-ak)/b(k,k);
end;
end;
Trong thöïc teá, öùng duïng cuûa phöông phaùp Choleski ñeå giaûi heä
phöông trình ñaïi soá tuyeán tính, ta chæ caàn tính ñoái xöùng vaø khoâng
caàn tính xaùc ñònh döông cuûa ma traän heä soá A. Khi ñoù caùc phaàn töû
cuûa ma traän tam giaùc B coù theå laø nhöõng soá phöùc (söû duïng ñôn vò aûo
√
i = −1). Tuy nhieân, keát quaû tính toaùn cuoái cuøng seõ cho chuùng ta
nghieäm thöïc cuûa heä phöông trình ban ñaàu.
52
HEÄ PHÖÔNG TRÌNH ÑAÏI SOÁ TUYEÁN TÍNH
Ví duï 3.7. Xeùt heä phöông trình
1
2
2
1
0
−1
0
1
x1
x2
x3
1
1
1
Ax =
=
= b
−1
Töø coâng thöùc (3.5) ta coù theå xaùc ñònh caùc heä soá bij, i > j cuûa ma
traän tam giaùc döôùi B nhö sau
√
2
2
√
√
b11 = 1, b21 = 2, b31 = −1, b22 = i 3, b32
=
, b33
=
i 3
3
Ta coù
1
1
2
0
0
0
√
1
1
1
1
√
i 3
2
−
i 3
By = b ⇔
y =
⇒ y =
2
2
√
√
−1
√
i 3
3
3
1
0
2
−1
1
√
2
1
√
√
i 3
2
i 3
−
BT x = y ⇔
x =
i 3
2
√
√
0
0
3
3
vaø ta coù nghieäm cuûa heä phöông trình laø x1 = 0, x2 = 1, x3 = 1.
3.4 CHUAÅN VECTÔ VAØ CHUAÅN MA TRAÄN
Xeùt khoâng gian tuyeán tính thöïc Rn. Chuaån cuûa vectô x ∈ Rn laø
moät soá thöïc, kyù hieäu laø kxk, thoaû caùc ñieàu kieän sau ñaây:
(i) ∀x ∈ Rn, kxk > 0, kxk = 0 ⇔ x = 0
(ii) ∀x ∈ Rn, ∀λ ∈ R, kλxk = |λ| .kxk
(iii) ∀x, y ∈ Rn, kx + yk 6 kxk + kyk. Ñieàu kieän naøy thöôøng ñöôïc goïi laø
baát ñaúng thöùc tam giaùc.
Trong Rn coù theå coù raát nhieàu chuaån, tuy nhieân chuùng ta chæ xeùt
chuû yeáu hai chuaån thöôøng duøng sau ñaây: ∀x = [x1, x2, . . ., xn]T ∈ Rn
n
X
kxk1 = |x1| + |x2| + · · · + |xn| =
|xk|
(3.6)
k=1
3.4 Chuaån vectô vaø chuaån ma traän
kxk∞ = max(|x1|, |x2| , . . ., |xn|) = max |xk|
53
(3.7)
k=1,n
Vieäc kieåm tra caùc coâng thöùc (3.6) vaø (3.7) thoaû caùc ñieàu kieän (i),
(ii), (iii) laø ñôn giaûn vaø daønh cho baïn ñoïc.
Baây giôø chuùng ta ñònh nghóa chuaån ma traän töông öùng vôùi chuaån
vectô.
Ñònh nghóa 3.1. Chuaån ma traän töông öùng vôùi chuaån vectô ñöôïc xaùc
ñònh theo coâng thöùc:
kAxk
kAk = kmxka=x1 kAxk = kmxka=x0 kxk
(3.8)
ꢂ
ꢃ
1
3
x1
x2
2
4
Ví duï 3.8. Xaùc ñònh chuaån cuûa ma traän A =
töông öùng
ꢂ
ꢃ
vôùi chuaån moät cuûa vectô. Vôùi moïi x =
kxk1 = |x1| + |x2| = 1, ta coù
∈ R2, sao cho
kAxk1 = |x1 + 2x2| + |3x1 + 4x2| 6 4 |x1| + 6 |x2| = 4 + 2 |x2| 6 6
Do ñoù kAk1 = 6.
Töø coâng thöùc (3.8) ta deã daøng suy ra ñöôïc raèng: kAxk 6 kAkkxk
Ñònh lí 3.4. Chuaån ma traän theo coâng thöùc (3.8) töông öùng vôùi chuaån
vectô ñöôïc xaùc ñònh nhö sau:
n
X
kAk1 = max
|aij|
|aij|
(3.9)
16j6n
i=1
n
X
kAk∞ = max
(3.10)
16i6n
j=1
∞
Ñònh nghóa 3.2. Xeùt daõy caùc vectô {x(m)
}
vôùi x(m) ∈ Rn. Ta noùi
m=0
daõy caùc vectô naøy hoäi tuï veà vectô x khi m → +∞ neáu vaø chæ neáu
kx(m) − xk → 0 khi m → +∞ (hoäi tuï theo chuaån).
54
HEÄ PHÖÔNG TRÌNH ÑAÏI SOÁ TUYEÁN TÍNH
Khi ñoù ta kyù hieäu lim x(m) = x vaø chuaån coù theå laáy moät chuaån
baát kyø trong caùc coângmt→hö∞ùc (3.6) hoaëc (3.7). Ta cuõng coù theå noùi daõy
vectô {x(m)} hoäi tuï veà x theo chuaån ñaõ cho. Ta coù ñònh lí sau ñaây:
Ñònh lí 3.5. Ñeå daõy caùc vectô {x(m)} hoäi tuï tôùi vectô x khi m → +∞
theo chuaån thì ñieàu kieän caàn vaø ñuû laø daõy {x(km)} hoäi tuï veà xk, ∀k = 1, n
(hoäi tuï theo toaï ñoä).
Baây giôø xeùt heä phöông trình Ax = b(det A = 0) coù nghieäm x = A−1b.
Cho b moät soá gia ∆b, khi ñoù nghieäm x töông öùng seõ coù soá gia laø ∆x,
vaø A∆x = ∆b ⇔ ∆x = A−1∆b. Ta coù
k∆xk = kA−1∆bk 6 kA−1kk∆bk
vaø
kbk = kAxk 6 kAkkxk
Töø ñaây chuùng ta deã daøng suy ra ñöôïc
k∆xk
kxk
k∆bk
kbk
6 kAk.kA−1
k
Soá
k(A) = Cond(A) = kAk.kA−1
k
(3.11)
ñöôïc goïi laø soá ñieàu kieän cuûa ma traän A. Ta coù theå chöùng toû ñöôïc raèng
1 6 k(A) 6 +∞
Soá ñieàu kieän cuûa ma traän ñaëc tröng cho tính oån ñònh cuûa heä thoáng
phöông trình ñaïi soá tuyeán tính. Giaù trò cuûa k(A) caøng gaàn vôùi 1 thì
heä caøng oån ñònh. Soá ñieàu kieän caøng lôùn thì heä caøng maát oån ñònh.
ꢂ
ꢃ
1
1
2
2.01
Ví duï 3.9. Xeùt heä phöông trình Ax = b vôùi A =
vaø
ꢂ
ꢃ
ꢂ
ꢃ
3
3.01
1
1
b =
. Deã thaáy heä coù nghieäm laø x =
. Baây giôø xeùt heä
ꢂ
ꢃ
ꢂ
ꢃ
3
3.1
−17
˜
˜
Ax˜ = b vôùi b =
. Nghieäm cuûa heä baây giôø laø x˜ =
.
10
˜
Ta nhaän thaáy k (A) = 1207.01 ꢀ 1. Do ñoù b ≈ b, nhöng x vaø x˜
∞
khaùc nhau raát xa.
Tải về để xem bản đầy đủ
30 trang
22640
Free
Bạn đang xem 20 trang mẫu của tài liệu "Giáo trình Phương pháp tính - Chương 3: Hệ phương trình đại số tuyến tính", để 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:
- giao_trinh_phuong_phap_tinh_chuong_3_he_phuong_trinh_dai_so.pdf
