矩陣程序大作業(yè)-說明文檔

LU分解、QR分解Householderreduction、Givensreduction輸入輸出參數(shù)的含義Factorization_and_Reduction.m中輸入輸出參數(shù)的含義：輸入?yún)?shù)：A：待分解的矩陣mode：工作模式的選擇LU分解reductionGivensreduction輸出參數(shù)：flag_right：分解或約減成功標(biāo)志，為0時表示失敗，為1時表示成功B1,B2,B3，在不同模式下的含義如下mode=0：B1=L：下三角陣；B2=U：上三角陣；B3=P：permutationmatrix；的單位正交基為列組成的矩陣；B2=R：上三角陣；0；滿足。由反射得到的單位正交陣(unitarymatrixororthogonalmatrix)；B2=T：上梯形陣；B3=P_T：P的轉(zhuǎn)置；PA=TA=P_T*T。(unitarymatrixororthogonalmatrix)；B2=T：上梯形陣；B3=P_T：P的轉(zhuǎn)置；A=P_T*T。試驗(yàn)步驟和原理LU分解a a a a11 12 1na22Aa2122

... a2na a ... a m1 m2 mnALU分解。和標(biāo)簽label〔用于記錄行交換后行的位置，將其初始化為label1 2 ... n。4〕j1n循環(huán)執(zhí)行以下步驟〔jj個主元：在第j列j,j以下位置中查找第一個非零元素作為第j個主元，即從ajin中查找第一個非零元素作為主元。ijjiALijlabelij個元素。AjUjL中j,j1。A0L的aajj

以下的aij

變?yōu)?時L中i,j)位置的元素l ij

ija 。jj1，其他元素0。A是否為奇異陣：查看矩陣U0元素，假設(shè)UAALU分解。PALU是否成立，假設(shè)成立，則分解正確。QR分解a a a 11 12 1n22〔1〕輸入待分解的矩陣：Aa21 a22

... a2na am1 m2

... amn1n循環(huán)執(zhí)行以下步驟：k=1uA〔A*1 *1

A的第一列；

k1

|A Q A

k

TA

Q

表示矩陣A*k i1

*k *i

i1

*k *i *kkQ*i

Qi列。u矩陣Q的第k列：Q uu*kRkRkj

Q |A*k

(jk,k1,,n)QR分解的結(jié)果是否正確：Q1，即驗(yàn)證QTQI是否成立。AQR是否成立。假設(shè)以上兩個條件均滿足，則QR分解正確。Householderreductiona a a 11 12 1n22輸入待處理的矩陣：Aa21 a22

... a2na am1 m2PI。

... a mn1n循環(huán)執(zhí)行以下步驟：以下圖紅框所示： 11 12 1n21 22 2n a’a’ ... Akm1 m2 mn令uAk*1

Ak e*1 1?I2uuTuTu

I 0 構(gòu)造m*m的矩陣R0 ? PPRPP*Ak+1k+1AkPT=P*A。Householderreduction的結(jié)果是否正確：P是否是單位正交矩陣(unitarymatrixororthogonalmatrix)，即驗(yàn)證PTPPPTI是否成立。APTT是否成立。Householderreduction分解正確。Givensreductiona a a 11 12 1n22輸入待處理的矩陣：Aa21 a22

... a2na am1 m2PI。

... amn1n循環(huán)執(zhí)行以下步驟：利用旋轉(zhuǎn)矩陣，使得原矩陣k列中k行元素以下的元素為0jk+1m循環(huán)執(zhí)行以下步驟：AkP*AkkAk 2Akkk2jkAk 2AkkkkkAk 2Akkk2jkAk 2Akkk2jk?Ic?kks?中kj列的元素，用-sjkc中jj列的元素。PP?PAkAkPAPT=P*A。Givensreduction的結(jié)果是否正確：驗(yàn)證P是否是單位正交陣(unitarymatrixororthogonalmatrix)，即驗(yàn)證PTPPPTI是否成立。APTT是否成立。Givensreduction分解正確。試驗(yàn)結(jié)果截圖當(dāng)輸入矩陣A為方陣0 20 14 例：A3 27 4 4 11 2LU分解顯示結(jié)果中變量含義：L：下三角陣；U：上三角陣；P：permutationmatrix。QR分解顯示結(jié)果中變量含義：Q：R(A)的單位正交基為列組成的矩陣；R：上三角陣。Householderreduction顯示結(jié)果中變量含義：P_H：由反射得到的單位正交陣(unitarymatrixororthogonalmatrix)；T_H：上梯形陣；P_T_H：P_H〔后綴_H用于與Givens的結(jié)果區(qū)分〕Givensreduction顯示結(jié)果中變量含義：P_G：由旋轉(zhuǎn)得到的單位正交陣(unitarymatrixororthogonalmatrix)；T_G：上梯形陣；P_T_G：P_G〔后綴_G用于與Householder的結(jié)果區(qū)分〕當(dāng)輸入矩陣A不是方陣4 3 42A2

14 14 01 7 151LU分解QR分解顯示結(jié)果中變量含義：Q：R(A)的單位正交基為列組成的矩陣；R：上三角陣。Householderreduction顯示結(jié)果中變量含義：P_H：由反射得到的單位正交陣(unitarymatrixororthogonalmatrix)；T_H：上梯形陣；P_T_H：P_H〔后綴_H用于與Givens的結(jié)果區(qū)分〕Givensreduction顯示結(jié)果中變量含義：P_G：由旋轉(zhuǎn)得到的單位正交陣(unitarymatrixororthogonalmatrix)；T_G：上梯形陣；P_T_G：P_G〔后綴_G用于與Householder的結(jié)果區(qū)分〕程序代碼RUN_FR.mclc;clearall;%趙曉梅2023280146280662023.11.19%供給幾個待分解的矩陣%A=[4,-3,4;2,-14,-3;-2,14,0;1,-7,15];1——非方陣A=[0,-20,-14;3,27,-4;4,11,-2];%待分解矩陣2——方陣fprintf(”待處理的矩陣為\n”);A%CommandWindow中顯示待處理矩陣fprintf(”**********************************************\n”);[L,U,P,flag_right_0]=Factorization_and_Reduction(A,0%LU分解ifflag_right_0==1%假設(shè)分解正確，顯示結(jié)果fprintf(”LU分解結(jié)果：\n”);L%下三角陣；U%上三角陣；%permutationmatrix；endfprintf(”**********************************************\n”);[Q,R,B3,flag_right_1]=Factorization_and_Reduction(A,1%QR分解ifflag_right_1==1%假設(shè)分解正確，顯示結(jié)果fprintf(”QR分解結(jié)果：\n”);%R(A)的單位正交基為列組成的矩陣；R%上三角陣；endfprintf(”**********************************************\n”);[P_H,T_H,P_T_H,flag_right_2]=Factorization_and_Reduction(A,2%Householderreductionifflag_right_2==1%Householderreduction正確，顯示結(jié)果(”顯示rn〔后綴H用于與s區(qū)分”P_H%由反射得到的單位正交陣；T_H%上梯形陣；P_T_H%P_H的轉(zhuǎn)置；endfprintf(”**********************************************\n”);[P_G,T_G,P_T_G,flag_right_3]=Factorization_and_Reduction(A,3%Givensreductionifflag_right_3==1%Givensreduction正確，顯示結(jié)果(”顯示sn〔后綴G用于與r區(qū)分”P_G%由反射得到的單位正交陣；T_G%上梯形陣；P_T_G%P_G的轉(zhuǎn)置；endFactorization_and_Reduction.m%------------------------------------%函數(shù)名稱：Factorization_and_Reduction%LU分解，QR分解，Householdreduction，Givensreduction%輸入?yún)?shù)：A：待分解的矩陣% mode：工作模式的選擇% mode=0：函數(shù)進(jìn)展LU分解% mode=1：函數(shù)進(jìn)展QR分解% mode=2Householderreduction% mode=3Givensreduction%輸出參數(shù)：flag_right：分解或約減成功標(biāo)志，為0時表示失敗，為1時表示成功% B1,B2,B3，在不同模式下的含義如下% mode=0：B1=L：下三角陣；% B2=U：上三角陣；% B3=P：permutationmatrix；% PA=LU。% mode=1：B1=Q：R(A)的單位正交基為列組成的矩陣；% B2=R：上三角陣；% B3=00；% A=QR。% mode=2：B1=P：由反射得到的單位正交陣(unitarymatrixororthogonalmatrix)；% B2=T：上梯形陣；% B3=P_T：P的轉(zhuǎn)置；% PA=TA=P_T*T。% mode=3：B1=P：由旋轉(zhuǎn)得到的單位正交陣(unitarymatrixororthogonalmatrix)；% B2=T：上梯形陣；% B3=P_T：P的轉(zhuǎn)置；% PA=TA=P_T*T。%根本信息：趙曉梅2023280146280662023.11.19%------------------------------------function[B1,B2,B3,flag_right]=Factorization_and_Reduction(A,mode)B1=0;B2=0;B3=0;flag_right=0;%當(dāng)不滿足分解或約減條件時，為保證有輸出，將參數(shù)初0[m,n]=size(A);%A的行數(shù)和列數(shù)%--------------------------------------------------------------------ifmode==0%LU分解fprintf(”LU分解……\n”);A0=A;%AP*A=L*UA0A的原值ifm~=nfprintf(”輸入矩陣不是方陣，不能進(jìn)展LU分解\n”);endL=zeros([n,n]);%L,U,P初始化U=zeros([n,n]);P=zeros([n,n]);label=(1:1:n);%標(biāo)簽，標(biāo)記進(jìn)展的行變換t=zeros(1,n);%用于行交換時的中間變量forj=1:n%jj個主元fori=j:nifA(i,j)~=0%j0的主元break;endend行交換

ifi~=j%0的主元不在位置〔j,j〕,ij行進(jìn)t=A(j,:);A(i,:)=t;l=label(j);%label相應(yīng)位置進(jìn)展交換label(i)=l;t=L(j,:);%L的相應(yīng)位置交換L(i,:)=t;

endU(j,:)=A(j,:);%AjUj行L(j,j)=1;%L的〔j,j〕1fori=(j+1):n%將主元以下變成0L(i,j)=A(i,j)/A(j,j);endfori=1:nP(i,label(i))=1;endfori=1:nifU(i,i)==0;%U0U為奇異陣，A也為奇異陣flag_right=0;endendifflag_right==0;fprintf(”A是奇異矩陣不能進(jìn)展LU分解\n”);else%LU分解是否正確fprintf(”驗(yàn)證結(jié)果是否正確……\n”);ifP*A0==L*Ufprintf(”P*A=L*U，LU分解正確進(jìn)展\n”);elsefprintf(”P*AL*U，LU分解錯誤\n”);endend

%--------------------------------------------------------------------ifmode==1%QR分解fprintf(”QR分解……\n”);% [m,n]=size(A);Q=zeros([m,n]);%初始化Q、R的大小R=zeros([n,n]);fork=1:nu=A(:,k);ifk>1fori=1:k-1%AkQ1~k-1列上的投影得uend

u=u-Q(:,i)”*A(:,k)*Q(:,i);v=norm(u);Q(:,k)=u/v;%uQk列forj=k:n%RkAjQk列構(gòu)成的向量上的投影的大小R(k,j)=Q(:,k)”*A(:,j);endend%驗(yàn)證分解結(jié)果是否正確fprintf(”驗(yàn)證結(jié)果是否正確……\n”);I=zeros([n,n]);fori=1:nI(i,i)=1;endeps=10^(-10);ifnorm(Q”*Q-I)<eps%Q的各列是否正交且各列的模為1QQ是否是單位陣flag_right_1=1;fprintf(”1.QQ是單位陣,Q各列相互正交且各列的模為1\n”);elseflag_right_1=0;fprintf(”1Q的轉(zhuǎn)置乘以Q不是單位陣,Q1\n”);endifnorm(A-Q*R)<epsfprintf(”2.A=Q*R\n”);elseflag_right_2=0;fprintf(”2.A=Q*R\n”);endfprintf(”所以，QR分解錯誤\n”);flag_right=0;elsefprintf(”所以，QR分解正確\n”);flag_right=1;endend%--------------------------------------------------------------------ifmode==2%Householderreductionfprintf(”Householderreduction……\n”);I0=zeros([m,m]);fori=1:mI0(i,i)=1;endP=I0;%P初始化為單位矩陣Ak=A;fork=1:n[mm,nn]=size(Ak);%Akkk列元素下方和右方的元素組成的矩陣ifmm==1break;ende(1)=1;u=Ak(:,1)-norm(Ak(:,1))*e;I=zeros([m-k+1,m-k+1]);fori=1:m-k+1I(i,i)=1;endR_hat=I-2*u*u”/(u”*u);%獲得能使Ak矩陣第一列只有第一個元素不為0的矩陣R=I0;%R_hatR矩陣的右下方，此時的RP*Akk行以下的元素全為0fori=(m-m_R+1):mforj=(m-n_R+1):mR(i,j)=R_hat(i-m+m_R,j-m+n_R);endendP=R*P;%RP相乘Ak=P*A;%PAkk列以下的元素構(gòu)成的矩陣更AkAk_hat=zeros(m-k,n-k);fori=k+1:mforj=k+1:nAk_hat(i-k,j-k)=Ak(i,j);endendAk=Ak_hat;endT=P*A;%Householderreduction結(jié)果是否正確fprintf(”驗(yàn)證結(jié)果是否正確……\n”);I=zeros([m,m]);fori=1:mI(i,i)=1;endeps=10^(-10);ifnorm(P”*P-I)<eps&&norm(P*P”-I)<eps%PP的轉(zhuǎn)置，即PPP是否是單位正交矩陣flag_right_1=1;fprintf(”1.P_T*P=IP*P_T=I,P確實(shí)是單位正交矩陣(unitarymatrixororthogonalmatrix)〔P_TP的轉(zhuǎn)置〕\n”);elseflag_right_1=0;fprintf(”1.P_T*P=IP*P_T=I，P不是是單位正交矩陣(unitarymatrixororthogonalmatrix)〔P_TP的轉(zhuǎn)置〕\n”);endifnorm(A-P”*T)<epsfprintf(”2.A=P_T*T(P_TP的轉(zhuǎn)置)\n”);elseflag_right_2=0;fprintf(”2.A=P_T*T(P_TP的轉(zhuǎn)置)\n”);endifflag_right_1==0||flag_right_2==0fprintf(”所以，Householderreduction錯誤\n”);flag_right=0;elsefprintf(”所以，Householderreduction正確\n”);flag_right=1;B1=P;B2=T;endend%-----------------------------------------------------------------------ifmode==3%Givensreductionfprintf(”Givensreduction……\n”);T=zeros