數(shù)值線性代數(shù)大作業(yè)報告

包頭師范學(xué)院數(shù)學(xué)科學(xué)學(xué)院2011級信息與計算科學(xué)專業(yè)數(shù)值線性代數(shù)實(shí)驗大報告指導(dǎo)老師：趙國忠姓名：1108300001劉帥1108300004王敏1108300032郭蒙實(shí)驗名稱：16題P75上機(jī)習(xí)題實(shí)驗?zāi)康模壕幹仆ㄓ玫淖映绦?，完成?xí)題的計算任務(wù)實(shí)驗內(nèi)容與要求：P75上機(jī)習(xí)題先用熟悉的計算機(jī)語言將算法2.5.1編制成通用的子程序，然后再用所編制的子程序完成下面兩個計算任務(wù)：估計5到20階Hilbert矩陣的無窮范數(shù)條件數(shù)。設(shè)A=先隨機(jī)地選取x∈R，并計算出b=Ax;然后再用列主元Gauss消去法求解該方程組，假定計算解為.試對n從5到30估計計算解的精度，并且與真實(shí)的相對誤差作比較。實(shí)驗原理：（1）矩陣范數(shù)（martixnorm）是數(shù)學(xué)上向量范數(shù)對矩陣的一個自然推廣。利用for循環(huán)和cond（a）Hilbert求解Hilbert矩陣的無窮范數(shù)，再利用norm(a,inf)求矩陣的無窮范數(shù)條件數(shù)。（2）本題分為4步來求解。先運(yùn)用rand隨機(jī)選取x∈R，輸入A矩陣，編制一個M文件計算出b。第二步用列主元高斯消去法求解出方程的解X2。第三步建立M文件:soluerr.m估計計算解的精度。第四步,建立M文件:bijiao.m,與真實(shí)相對誤差作比較。實(shí)驗過程：(1)程序：clearforn=5:20fori=1:nforj=1:na(i,j)=1/(i+j-1);endendc=cond(a);f=norm(c,inf);fprintf('n=%3.0f\nnorm(c,inf)%e\n',n,f)end運(yùn)行結(jié)果：n=5norm(c,inf)4.766073e+005n=6norm(c,inf)1.495106e+007n=7norm(c,inf)4.753674e+008n=8norm(c,inf)1.525758e+010n=9norm(c,inf)4.931542e+011n=10norm(c,inf)1.602467e+013n=11norm(c,inf)5.224376e+014n=12norm(c,inf)1.698855e+016n=13norm(c,inf)3.459404e+017n=14norm(c,inf)4.696757e+017n=15norm(c,inf)2.569881e+017n=16norm(c,inf)7.356249e+017n=17norm(c,inf)4.362844e+017n=18norm(c,inf)1.229633e+018n=19norm(c,inf)9.759023e+017n=20norm(c,inf)1.644051e+018（2）程序：M文件：matrix1.mfunction[a,b,x1]=matrix1(n)formatlongA1=-1*ones(n,n)A2=tril(A1)fori=1:nA2(i,i)=1endA2(:,n)=1a=A2x1=rand(n,1)b=A2*x1end運(yùn)行結(jié)果：>>A1=-1-1-1-1-1-1-1-1-1-1-1-1-1-1-1-1-1-1-1-1-1-1-1-1-1A2=-10000-1-1000-1-1-100-1-1-1-10-1-1-1-1-1A2=10000-1-1000-1-1-100-1-1-1-10-1-1-1-1-1A2=10000-11000-1-1-100-1-1-1-10-1-1-1-1-1A2=10000-11000-1-1100-1-1-1-10-1-1-1-1-1A2=10000-11000-1-1100-1-1-110-1-1-1-1-1A2=10000-11000-1-1100-1-1-110-1-1-1-11A2=10001-11001-1-1101-1-1-111-1-1-1-11a=10001-11001-1-1101-1-1-111-1-1-1-11x1=0.8147236863931790.9057919370756190.1269868162935060.9133758561390190.632359246225410b=1.4470829326185890.723427496907850-0.961169560949882-0.301767337397875-2.128519049675914a=10001-11001-1-1101-1-1-111-1-1-1-11b=1.4470829326185890.723427496907850-0.961169560949882-0.301767337397875-2.128519049675914x1=0.8147236863931790.9057919370756190.1269868162935060.9133758561390190.632359246225410M文件：LZYgauss.mfunction[x2]=LZYgauss(a,b)formatlongn=length(a);x2=zeros(n,1);a=[ab];fork=1:n-1max=k;fori=k+1:nifa(i,k)>a(max,k)max=i;endendtemp=a(k,k:n+1);a(k,k:n+1)=a(max,k:n+1);a(max,k:n+1)=temp;fori=k+1:na(i,k)=-a(i,k)/a(k,k);a(i,k+1:n+1)=a(i,k+1:n+1)+a(i,k)*a(k,k+1:n+1);endendx2(n,1)=a(n,n+1)/a(n,n);fori=n-1:-1:1sum=0;forj=i+1:nsum=sum+x2(j,1)*a(i,j);endx2(i,1)=(a(i,n+1)-sum)/a(i,i);end運(yùn)行結(jié)果：>>LZYgauss(a,b)ans=0.8147236863931790.9057919370756190.1269868162935060.9133758561390200.632359246225410估計計算解的精度：M文件:soluerr.mfunction[x,error]=soluerr(a,b)formatlong%估計計算解的精度%算法：列主元Gauss消去法，其中%A---系數(shù)矩陣%b-右端項%index---index=0表示計算成??；index=1表示計算成功%輸出結(jié)果：error--本算法給出的計算解的估計%normA--逆矩陣無窮范數(shù)估計%rnorm--計算解的殘量[n,m]=size(a);nb=length(b);ifn~=merror('Therowsandcolumnsofmatrixamustbeequal!');return;endifm~=nberror('Thecolumnsofamustbeequalthedimensionofb!');return;endindex=1;%列主元矩陣三角分解[L,U,u,index_col]=Gauss_col(a);%解下三角方程組Ly=Pb[y,index_low]=Gauss_low(L,b(u));%解上三角方程組Ux=y[x,index_upp]=Gauss_upp(U,y);%輸出數(shù)值解xpause(0.3)%估計矩陣逆的無窮大范數(shù)normA=normAinv(L,U,u);%估計計算解的殘量rnorm=norm(b-a*x,inf);%計算右端項bnorm=norm(b,inf);%計算矩陣A的范數(shù)Anorm=norm(a,inf);%計算解的精度error=normA*Anorm*rnorm/bnorm;運(yùn)行結(jié)果：>>[x,error]=soluerr(a,b)x=0.8147236863931790.9057919370756190.1269868162935060.9133758561390200.632359246225410x=0.8147236863931790.9057919370756190.1269868162935060.9133758561390200.632359246225410error=5.215941218821592e-016X1與真實(shí)相對誤差做比較M文件：bijiao.mfunction[x,x1,error,error1]=bijiao(a,b,n)[a,x1]=matrix1(n)[x,error]=soluerr(a,b)error1=abs((x-x1)/x1)運(yùn)行結(jié)果：>>x=0.8147236863931790.9057919370756190.1269868162935060.9133758561390200.632359246225410x1=0.9578935505663481.6132601689718680.629241553693582-0.799716694069468-1.770467991515150error=5.215941218821592e-016error1=Columns1through4000000