数值线性代数大作业报告

包头师范学院数学科学学院2011级信息与计算科学专业数值线性代数实验大报告指导老师：赵国忠姓名：1108300001刘帅1108300004王敏1108300032郭蒙实验名称：16题P75上机习题实验目的：编制通用的子程序，完成习题的计算任务实验内容与要求：P75上机习题先用熟悉的计算机语言将算法2.5.1编制成通用的子程序，然后再用所编制的子程序完成下面两个计算任务：估计5到20阶Hilbert矩阵的无穷范数条件数。设A=先随机地选取x∈R，并计算出b=Ax;然后再用列主元Gauss消去法求解该方程组，假定计算解为.试对n从5到30估计计算解的精度，并且与真实的相对误差作比较。实验原理：（1）矩阵范数（martixnorm）是数学上向量范数对矩阵的一个自然推广。利用for循环和cond（a）Hilbert求解Hilbert矩阵的无穷范数，再利用norm(a,inf)求矩阵的无穷范数条件数。（2）本题分为4步来求解。先运用rand随机选取x∈R，输入A矩阵，编制一个M文件计算出b。第二步用列主元高斯消去法求解出方程的解X2。第三步建立M文件:soluerr.m估计计算解的精度。第四步,建立M文件:bijiao.m,与真实相对误差作比较。实验过程：(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运行结果：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运行结果：>>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运行结果：>>LZYgauss(a,b)ans=0.8147236863931790.9057919370756190.1269868162935060.9133758561390200.632359246225410估计计算解的精度：M文件:soluerr.mfunction[x,error]=soluerr(a,b)formatlong%估计计算解的精度%算法：列主元Gauss消去法，其中%A---系数矩阵%b-右端项%index---index=0表示计算成败；index=1表示计算成功%输出结果：error--本算法给出的计算解的估计%normA--逆矩阵无穷范数估计%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);%输出数值解xpause(0.3)%估计矩阵逆的无穷大范数normA=normAinv(L,U,u);%估计计算解的残量rnorm=norm(b-a*x,inf);%计算右端项bnorm=norm(b,inf);%计算矩阵A的范数Anorm=norm(a,inf);%计算解的精度error=normA*Anorm*rnorm/bnorm;运行结果：>>[x,error]=soluerr(a,b)x=0.8147236863931790.9057919370756190.1269868162935060.9133758561390200.632359246225410x=0.8147236863931790.9057919370756190.1269868162935060.9133758561390200.632359246225410error=5.215941218821592e-016X1与真实相对误差做比较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)运行结果：>>x=0.8147236863931790.9057919370756190.1269868162935060.9133758561390200.632359246225410x1=0.9578935505663481.6132601689718680.629241553693582-0.799716694069468-1.770467991515150error=5.215941218821592e-016error1=Columns1through4000000