數(shù)值分析實習(xí)報告5

數(shù)值分析實習(xí)報告學(xué)院：數(shù)學(xué)與統(tǒng)計學(xué)院專業(yè)：數(shù)學(xué)類專業(yè)班姓名：楊圣學(xué)號：320210928991編寫語言：matlab問題一：.用LU分解及列主元Gauss消去法求解線性方程組:10-35210-352-72.099999—115.90000151⑴輸出lu分解法的分解矩陣上和&解向量局；⑵輸出列主元Gauss消去法的行交換次序,解向量均;⑶比較兩種方法所得的結(jié)果.function[L,U]=LU(A)[n,n]=size(A);e=eye(n);L=e;fork=l:n-lLk=eye(n);fori=k+l:nLk(i,k)=-A(i,k)/A(k,k);endL=Lk*L;A=Lk*A;U=A;L=inv(L);endA=[10,-7,0,1；-3,2.099999,6,2；5,-1,5,-1；2,1,0,2]；b=[8,5.900001,5,1][L,U]=LU(A)x=inv(L*L))*b解得：L1 3.70E-16 -9.87E-29 6.45E+39TOC\o"1-5"\h\z-0.3 1 0 00.5 -2500000 1 -6.89E-170.2 -2400000 0.9599997 1U10 -7 0 10 -1.00E-06 6 2.30 0 15000005 5749998.4990 0 0 5.079998908X 1.20E-11 -1 1 1function[X]=rowElimination(AJb)%系數(shù)矩陣Adoublen*n,常數(shù)項bdoublen*l%解向量Xdoublen*lAb=[Ab];%doublen*(n+l)[n,n]=size(A);P=eye(n);%向量p代表置換矩陣，左乘即代表交換行后原來每行的位置L=eye(n);fori=l:n-l[aJk]=max(abs(Ab(i:nJi)));%對位置(i,i)及以下的列提取列主元a和列主元位置kifa==0error('Error:第%d步的列主元為0!\n"i);%最大值為0證明不滿秩，做不了endk=k+i-l;%1＜是相對于i的位置，所以這里需要更正為矩陣的真實位置(k,i)ifk~=i%如果原矩陣(i,i)位置不是列主元，就需要交換第i行和第k行，同時在A和p上操作tmp=Ab(iJ:);Ab(i,:)=Ab(k,:);Ab(k,:)=tmp;tmp=p(i,:);p(i,:)=p(k,:);p(k,:)=tmp;fprintf(1在第％d步迭代時交換行%€1和％心門"i,i,k);endL(i+l:n^i)=Ab(i+l:n,i)/Ab(i,i);%計算L的位置(i,i)及以下的列，這里對L來說是加號，對Lk是減號Ab(i+l:n,i:n+l)=Ab(i+l:n,i:n+l)-L(i+l:n,i)*Ab(i,i:n+l);end%至此，得到上三角增廣矩陣Abdoublen*(n+l),下面回代解這個方程組A=Ab(:,l:n);b=Ab(:,n+1);X=zeros(nJl);X(n,l)=b(n,l)/A(n,n);fori=n-1:-1:1X(i,l)=(b(i,l)-A(i,i+l:n)*X(i+l:n,l))/A(i,i);endend解得：X2.66453525910038e-16-111(3)比較：因為準(zhǔn)確值是[0,-1,1,1],列主元消去法得到的結(jié)果是le-16量級的,比單純的Gauss消去法的le-11量級要精確很多。所以，列主元消去法更為精確。問題二：2.計算下列矩陣的Choiesky分解：TOC\o"1-5"\h\z4 2 4 2-2 10 8 14 8 9 5 ,_21519問題分析：比較簡單，按照書上偽代碼直接編寫function[L]=cholesky(A)[n,n]=size(A);L=zeros(n);fori=l:nforj=l:i-lendL(i,i)=(A(i,i)-L(i/：i-l)*L(i,l：i-l”》(0.5)；endendA=[4,2,4,2；2"0,8,l；4,8,9,5；2,l,5,19]；[L]=cholesky(A)解得TOC\o"1-5"\h\z2 0321 00 00 01 03 3問題三:3.用追趕法求解下列三對角方程組:2-10—12—10-120 0 -10 0 0問題分析：比較簡單，按照書上偽代碼直接編寫function[X]=chasemethod(A^b)%X解向量，系數(shù)矩陣A,常數(shù)項b[n,n]=size(A);r=zeros(nJl);f=zeros(n-lJl);l=zeros(nJl);%LU分解fori=l:n-lf(i)=A(i,i+l);endr(l)=A(l,l);fori=2:nl(i)=A(i,i-l)/r(i-l);r(i)=A(i,i)-l(i)*f(i-l);end%追過程y=zeros(nJl);y(l)=b(l,l);fori=2:r)y(i)=b(i)-l(i)*y(i-l);end%趕過程X=zeros(nJl);X(n)=y(n)/r(n);fori=n-l:-1:1X(i)=(y(i)-f(i)*X(i+l))/r(i);endendA=[2-1000；