Thứ Năm, 20 tháng 2, 2014
Tài liệu CHƯƠNG 2: MA TRẬN ppt
62
[]
[]
[]
[]
[][]
[]
[
]
[]
[]
[]
[][]
⎛⎞
′
′′ ′ ′
= − =− =−
⎜⎟
⎝⎠
T
TT
AU
UU
AH A E A U A VU
Trongđó:
[]
[]
[]
′
=
AU
V
Q
(8)
Dovậy:
[]
[]
[] []
[][]
[]
[][]
()
⎛⎞
′′
=− −
⎜⎟
⎝⎠
T
T
UU
HA H E A VU
Q
[]
[][]
[
]
[
]
[]
[][]
()
′′
=− − −
T
TT
UU
AVU AVU
Q
[]
[][]
[][]
[
]
(
)
[][][]
(
)
[]
′
′
=− − +
TTT
T
UUA UUVU
AVU
[]
[]
[
]
[
]
[
]
[
]
[
]
′
=−−+
TT T
AVU UV2gUU
Trongđó:
[][]
=
T
UV
g
2Q
(9)
Đặt:[W]=[V]‐g[U](10)
Tathấyngayphépbiếnđổicódạng:
[]
[]
[]
[]
[]
[
]
[
]
[
]
′′
=− −
TT
HA H A WU UW(11)
Thuậttoáncóthểtómlạinhưsau:
‐Cho[A’]làmatrậnvuôngcấp(n‐i)cóđượctừphầndướibênphải
củamatrận[A]
‐Đặt
++
=
⎡⎤
⎡⎤
⎣⎦
⎣⎦
L
T
i 1,i i 2,i n,i
Xa a a
‐Tính
[]
X.Chok=
[]
X nếux1>0vàk=‐
[
]
X nếux1<0
‐Cho
−
=+
⎡⎤
⎡⎤
⎣⎦
⎣⎦
L
T
12 ni
Ukxx x
‐Tính
[]
=
2
U
Q
2
‐Tính
[]
[]
[]
′
=
AU
V
Q
‐Tính
[][]
=
T
UV
g
2Q
63
‐Tính[W]=[V]‐g[U]
‐Tính
[]
[]
[][]
[
]
[
]
′
=− −
TT
AAWU UW
‐Đặt
++
==−
i,i 1 i 1,i
aa k
Taxâydựnghàm
housetrans()đểthựchiệnthuậttoántrên:
functionA=housetrans(A)
%BiendoiHouseholdermatranAthanhmatran
%bađườngchéodang[c\d\c].
%Decocvaddungd=diag(A),c=diag(A,1).
n=size(A,1);
fork=1:n‐2
u=A(k+1:n,k);
uMag=sqrt(dot(u,u));
ifu(1)<0;
uMag=‐uMag;
end
u(1)=u(1)+uMag;
A(k+1:n,k)=u;%LuuuvaophanduoicuaA.
H=dot(u,u)/2;
v=A(k+1:n,k+1:n)*u/H;
g=dot(u,v)/(2*H);
v=v‐g*u;
A(k+1:n,k+1:n)=A(k+1:n,k+1:n)‐v*uʹ‐
u*vʹ;
A(k,k+1)=‐uMag;
end
k=zeros(n);
fori=1:n
k(i,i)=A(i,i);
end
fori=1:n‐1
k(i,i+1)=A(i,i+1);
k(i+1,i)=A(i,i+1);
end
A=k;
64
Để tính ma trận bađường chéo theo phép biếnđổi Householder ta dùng
chươngtrình
cthousetrans.m:
clearall,clc
a=[1234;2935;3337;4576];
b=householder(a)
d=diag(b)
c=diag(b,1)
§3.BIẾNĐỔITHÀNHMATRẬNHESSENBERG
Nếumatrận[A]làmatrậnđốixứng,phươngphápHouseholdercóthể
đượcsửdụngđểbiếnđổinóthànhmatrậnđồngdạngđốixứngbađường
chéo.Nếumatrận[A]khôngđốixứ
ng,phươngphápHouseholderbiếnđổi
matrận[A]thànhmatrậnđồngdạngHessenberg.
MatrậnHessenberglàmatrậncódạng:
[]
11 12 13 1,n
21 22 23 2n
32 33 2n
nn
aaa a
aaa a
0a a a
H
000 a
⎡⎤
⎢⎥
⎢⎥
⎢⎥
=
⎢⎥
⎢⎥
⎢⎥
⎣⎦
L
L
L
MMML M
L
TathựchiệnphépbiếnđổiHouseholdertrênmatrận[A]vàcóđược:
[Q][H][Q’]=[A]
trongđó[Q]làmatrậntrựcgiao(tagọiđâylàphântíchHessenbergmatrận
[A]).
Thuậttoáncóthểtóm
lạinhưsau:
‐Cho[Q]làmatrậnđơnvịcấpn
‐Đặt
+
=
⎡⎤
⎡⎤
⎣⎦
⎣⎦
L
T
i2,i n,i
X0a a
‐Tính
[]
X.Choα=
[]
X nếuai+2,i>0vàα=‐
[
]
X nếuai+2,i<0
‐Cho
−
=α+
⎡⎤
⎡⎤
⎣⎦
⎣⎦
L
T
2ni
U0 x x
‐Tính
[]
β=
2
U
2
‐Tính
[] []
[]
[]
′
=−
β
UU
PE
65
‐Tính
[]
[][]
=QQP
‐Tính
[
]
[][ ][]
=APAP
Taxâydựnghàm
hessenberg()đểthựchiệnphépphântíchtrên:
function[H,Q]=hessenberg(a)
[n,n]=size(a);
q=eye(n);
fork=1:n‐2
alfa=0;
forj=k+1:n
alfa=alfa+a(j,k)^2;
end
alfa=sign(a(k+1,k))*sqrt(alfa);
u=zeros(1,n);
u(k+1:n)=a(k+1:n,k);
u(k+1)=u(k+1)+
alfa;
beta=.5*u*uʹ;
p=eye(n);
fori=1:n
p(i,1:n)=p(i,1:n)‐(u(i)*u(1:n))/beta;
end
q=q*p;
a=p*a*p;
end
H=a;
Q=q;
Đểphântíchmatrậntadùngchươngtrìnhcthessenberg.m:
clearall,clc
a=[1234;5674;648 9;3579];
[H,Q]=hessenberg(a)
§4.PHÂNTÍCHMATRẬNTHEOPHƯƠNGPHÁPDOOLITTLE
66
Mộtmatrậnkhôngsuybiến[A]gọilàphântíchđượcthànhtíchhaima
trận[L]và[R]nếu:
[A]=[L][R]
Việcphântíchnày,nếutồntại,làkhôngduynhất.
Nếuma
trận[L]cócácphầntửnằmtrênđườngchéochínhbằng1,tacó
phépphântíchDoolittle.
Nếumatrận[R]cócácphầntửnằmtrênđườngchéochínhb ằng1,ta
cóphépphântíchCrout.
Nếu
[R]=[L]
T
(hay[L]=[R]
T
)tacóphépphântíchCholeski.
Vớimatrậnbậc3,[L]và[R]códạng:
[] []
11 12 13
21 22 23
31 32 33
100 rrr
Ll 10 R0rr
ll1 00r
⎡⎤ ⎡ ⎤
⎢⎥ ⎢ ⎥
==
⎢⎥ ⎢ ⎥
⎢⎥ ⎢ ⎥
⎣⎦ ⎣ ⎦
Đểtìml
ijvàrijtathựchiệnphépnhân.Saukhinhântacó:
[]
11 12 13
11 21 12 21 22 13 21 23
11 31 12 31 22 32 13 31 23 32 33
rr r
Arlrlr rlr
rl rl rl rl rl r
⎡⎤
⎢⎥
=++
⎢⎥
⎢⎥
+++
⎣⎦
BâygiờtathựchiệnphépkhửGaussđốivớiphươngtrìnhtrên.Đầutiênta
chọnhàngthứnhấtlàmtrụvàthựchiênphépbiếnđổi:
hàng2‐l
21×hàng1(khửa21)→hàng2
hàng3‐l
31×hàng1(khửa31)→hàng3
kếtquảtacó:
[]
11 12 13
12223
22 32 23 32 33
rr r
A0r r
0rl rl r
⎡⎤
⎢⎥
=
⎢⎥
⎢⎥
+
⎣⎦
Sauđótalấyhàngthứhailàmtrụvàthựchiệnbiếnđổi:
hàng3‐l32
×hàng2(khửa32)→hàng3
vàcó:
[]
11 12 13
22223
33
rrr
A0rr
00r
⎡⎤
⎢⎥
=
⎢⎥
⎢⎥
⎣⎦
Nhưvậytathấyngayrằngmatrận[R]làmatrậncóđượckhithựchiện
loạitrừGausstiếnma trận[A]vàcácphầntửcủa[L]làcácnhânt
ửdùngkhi
67
loạitrừa ij.Điềuđócónghĩalàđểtìmmatrận[L]và[R]tadùngphépkhử
Gausstiến.Taxâydựnghàm
doolittle()đểthựchiệnloạiphântíchDoolittle.
function[l,r]=doolittle(A)
%PhantichmatranAthanhA=L*U
n=size(A,1);
u=zeros(n);
fork=1:n‐1
fori=k+1:n
ifA(i,k)~=0.0
lambda=A(i,k)/A(k,k);
A(i,k+1:n)=A(i,k+1:n)‐lambda*A(k,k+1:n);
A(i,k)=
lambda;
end
end
end
l=tril(A);
fori=1:n
l(i,i)=1;
end
l=triu(A);
fori=1:n
l(i,i)=A(i,i);
end
§5.PHÂNTÍCHMATRẬNTHEOPHƯƠNGPHÁPCROUT
TươngtựnhưthuậttoánDoolittle,tacóthểphântíchmatrận[A]theo
thuậttoánCroutthànhtíchcủamatrận[L]và[R].Cácmatrậnbậc3theo
Croutcódạng:
[] []
11 12 13
21 22 23
31 32 33
l00 1rr
Lll 0 R01r
lll 001
⎡⎤ ⎡⎤
⎢⎥ ⎢⎥
==
⎢⎥ ⎢⎥
⎢⎥ ⎢⎥
⎣⎦ ⎣⎦
Đểtìml
ijvàrijtathựchiệnphépnhân.Saukhinhântacó:
68
[]
11 11 12 11 13
21 21 12 22 21 13 22 23
31 31 12 32 31 13 32 23 33
llr lr
Allrllrlr
llrllrlrl
⎡⎤
⎢⎥
=++
⎢⎥
⎢⎥
+++
⎣⎦
Nhưvậy:
a
11=1.r11+0.0+0.0=r11;
a
12=r12;a13=r13
a
21=l21r11;
a
22=l21r12+r22;a23=l31r11
a
31=l31r11;a32=l31r12;
a
33=l31r13+l32r23+r33
Mộtcáchtổngquáttacó:
vớij>i: l
ij=rji=0
vớii=1: r
1j=a1j(j=1tớin)
l
j1=aj1/r11(j=1tớin)
vớii=2tớin
∑
−
=
−=
1i
1k
kjik
ijij
rlar (j=itớin)
ii
1i
1k
kijk
ji
ji
r
rla
l
∑
−
=
−
= (j=itớin)
Taxâydựnghàm
crout()đểphântíchmatrậntheothuậttoánCrout:
function[l,r]=crout(a)
n=size(a,1);
l=zeros(n);
r=zeros(n);
fori=1:n
r(1,i)=a(1,i);
l(i,i)=1.;
l(i,1)=a(i,1)/a(1,1);
end
fork=2:n
r(k,k:n)=a(k,k:n)‐l(k,1:k)*r(1:k,k:n);
if
k~=n
69
fori=1:n
l(i,k)=(a(i,k)‐l(i,1:k‐1)*r(1:k‐1,k))/r(k,k);
end
end
end
§6.PHÂNTÍCHMATRẬNTHEOPHƯƠNGPHÁPCHOLESKI
ThuậttoánCholeskichophépphântíchmatrận[A]thànhtíchhaima
trận:
[A]=[L][L]
T
.
Thuậttoánnàyđòihỏi:
‐[A]làmatrậnthực,đốixứng
‐[A]làmatrậnxácđịnhdương
Tavuông[A]cấp3theothuậttoánCholeski:
11 12 13 11 11 21 31
21 22 23 21 22 22 32
31 32 33 31 32 33 33
aaa l 00lll
aaa ll 00ll
aaa lll 00l
⎡⎤⎡⎤⎡⎤
⎢⎥⎢⎥⎢⎥
=
⎢⎥⎢⎥⎢⎥
⎢⎥⎢⎥⎢⎥
⎣⎦⎣⎦⎣⎦
Saukhithựchiệnphépnhântacó:
2
11 12 13 11 11 21 11 31
22
21 22 23 11 21 21 22 21 31 22 32
222
31 32 33 11 31 21 31 22 32 31 32 33
a a a l ll ll
aaa llll llll
aaa lllllllll
⎡⎤
⎡⎤
⎢⎥
⎢⎥
=+ +
⎢⎥
⎢⎥
⎢⎥
⎢⎥
+++
⎣⎦
⎣⎦
Vếphảilàmatrậnđốixứng.Cânbằngcácphầntửcủahaimatrậntacó:
11 11 21 21 11 31 31 11
222
22 22 21 32 32 21 31 22 33 33 31 32
l a la/l la/l
lall(all)/llall
== =
=− =− =−−
Tổngquát,vớimatrậncấpn,tacó:
[][]
()
j
T
i1 j1 i 2 j2 ik jk
ij
k1
LL ll ll ll i j
=
=++⋅⋅⋅+= ≥
∑
Cânbằngvớiphầntửcủamatrận[A]tacó:
j
ij ik jk
k1
a l l i j,j 1, ,n j 1,2, ,n
=
==+=
∑
Domatrận[L]làmatrậntamgiáctráinênđốivớicộtthứnhấttacó:
11 11 i1 i1 11
la la/l==
Đốivớicộtkhác,rútl
ijrakhỏitổngtacó:
70
j
1
ij ik jk ij jj
k1
allll
−
=
=+
∑
Nếui=j(phầntửtrênđườngchéo)thì:
j1
2
jj jj jk
k1
l a l j 2,3, ,n
−
=
=− =
∑
vàphầntửnằmngoàiđườngchéo:
j1
ij ij ik jk
k1
jj
1
l a l l j 2, 3, , n i j 2, j 3, ,n
l
−
=
⎛⎞
=− = =++
⎜⎟
⎝⎠
∑
Dựavàothuậttoántrêntaxâydựnghàm
choleski()
functionL=choleski(A)
%PhantichmatranathanhA=LL’.
%Cuphap:L=choleski(A)
f=posdef(A);
iff==0
error(ʹMatrankhongxacdinhduong!ʹ);
return
end
n=size(A,1);
forj=1:n
temp=A(j,j)‐
dot(A(j,1:j‐1),A(j,1:j‐1));
iftemp<0.0
error(ʹMatrankhongxacdinhduongʹ)
end
A(j,j)=sqrt(temp);
fori=j+1:n
A(i,j)=(A(i,j)‐dot(A(i,1:j‐1),A(j,1:j‐1)))/A(j,j);
end
end
L=tril(A);
functionf=posdef(M)
%Kiemtralieu
matranMcoxacdinhduonghaykong
isposdef=true;
71
fori=1:length(M)
if(det(M(1:i,1:i))<=0)
isposdef=false;
break;
end
end
f=isposdef;%0neusai,1neudung
§7.PHÂNTÍCHQRBẰNGTHUẬTTOÁNHOUSEHOLDER
Chomatrận[A],phântíchQRcủanóchota:
[A]=[Q]*[R]
Trongđó[Q]làmatrậntrựcgiaovà[R]làmatrậntamgiácphải.
TadùngbiếnđổiHouseholderđểtìmcácmatrận[Q]và[R].
[][][]
[]
[
]
−−
⋅
⋅⋅ =
n1 n2 1
HH HAR(1)
Nhưvậy:
[]
[][][]
(
)
[]
[] [ ][ ]
[]
[][ ][ ]
[][][]
−
−−
−− − −
−−
= ⋅⋅⋅ = ⋅⋅⋅
= ⋅⋅⋅ =
1
11
n1 n2 1 1 n2 n1
1n2n1
AHH H RH H HR
HHHRQR
(2)
TíchcủatấtcảcácmatrậnHouseholder:
[]
[][ ][ ]
−−
= L
1n2n1
QH H H(3)
khôngnhữngđốixứngmàcòntrựcgiaonhưmỗimatrận[H
k]:
[][]
[][ ][ ]
(
)
[][ ][ ]
[ ][ ] [][][ ][ ]
[]
−− −−
−− −−
=⋅⋅⋅ ⋅⋅⋅
= ⋅⋅⋅ ⋅⋅⋅ =
T
T
1n2n11n2n1
TT T
n1 n2 1 1 n2 n1
QQ H H H H H H
HH HHHH E
Taxâydựnghàm
qrdecom()đểphântíchmatrận:
function[Q,R]=qrdecom(A)
%PhantichQR
n=size(A,1);
R=A;
Q=eye(n);
fork=1:n‐1
H=householder(R(:,k),k);
R=H*R;%Pt.(1)
Q=Q*H;%Pt.(3)
Đăng ký:
Đăng Nhận xét (Atom)
Không có nhận xét nào:
Đăng nhận xét