2016矩阵与数值分析上机作业满分

第一题考虑计算给定向量的范数：输入向量(),输出‖‖‖‖‖‖。请编制一个通用程序，并用你编制的程序计算如下向量的范数：()()对n=10，100，1000甚至更大的n计算其范数，你会发现什么结果?你能否修改你的程序使得计算结果相对精确呢？代码：functionProblem1()clcN=input('inputN=');X=zeros(1,N);Y=1:N;fori=1:NX(i)=1/i;endx1=0;x2=0;y1=0;y2=0;fori=1:Nx1=x1+abs(X(i));x2=x2+X(i)*X(i);y1=y1+abs(Y(i));y2=y2+Y(i)*Y(i);endx1;x2=sqrt(x2);x3=max(abs(X));y1;y2=sqrt(y2);y3=max(abs(Y));fprintf('x的1范数是%.2f\n',x1);fprintf('x的2范数是%.2f\n',x2);fprintf('x的无穷范数是%.2f\n',x3);fprintf('y的1范数是%.2f\n',y1);fprintf('y的2范数是%.2f\n',y2);fprintf('y的无穷范数是%.2f\n',y3);结果：inputN=10x的1范数是2.93x的2范数是1.24x的无穷范数是1.00y的1范数是55.00y的2范数是19.62y的无穷范数是10.00inputN=100x的1范数是5.19x的2范数是1.28x的无穷范数是1.00y的1范数是5050.00y的2范数是581.68y的无穷范数是100.00inputN=1000x的1范数是7.49x的2范数是1.28x的无穷范数是1.00y的1范数是500500.00y的2范数是18271.11y的无穷范数是1000.00第二题考虑()(),其中定义f(0)=1，此时f(x)是连续函数。用此公式计算当[]时的函数值，画出图像。另一方面，考虑下面算法：d=1+xifd=1theny=1elsey=lnd/(d-1)endif用此算法计算[]时的函数值，画出图像。比较一下发生了什么？代码：functionproblem2()clcN=1000;n=2*10^(-15)/N;%公式计算f=zeros(1,N+1);t=1;fori=-10^(-15):n:10^(-15)if(i~=0)f(t)=log(1+i)/i;elsef(t)=1;endt=t+1;end%算法计算a=zeros(1,N+1);t=1;fori=-10^(-15):n:10^(-15)d=1+i;if(d~=1)a(t)=log(d)/(d-1);elsea(t)=1;endt=t+1;end%画图，左边是公式算的，右边是算法算的subplot(1,2,1);plot(-10^(-15):n:10^(-15),f);holdonsubplot(1,2,2);plot(-10^(-15):n:10^(-15),a);holdon结果：第三题首先编制一个利用秦九韶算法计算一个多项式在给定点的函数值的通用程序，你的程序包括输入多项式的系数以及给定点，输出函数值。利用你编制的程序计算()()在x=2邻域附近的值。画出()在[]上的图像。代码：functionproblem3()%秦九韶算法clcn=9;c=[1,-18,144,-672,2016,-4032,5376,-4608,2304,-512];N=100;m=(2.05-1.95)/N;y=zeros(1,N+1);ydex=0;forx=1.95:m:2.05ydex=ydex+1;dex=1;temp=c(dex);whiledex=ndex=dex+1;temp=temp*x+c(dex);endy(ydex)=temp;endplot(1.95:m:2.05,y);结果：第四题编制计算给定矩阵A的LU分解和PLU分解的通用程序，然后用你编制的程序完成下面两个计算任务：(1)考虑[]自己取定，并计算b=Ax。然后用你编制的不选主元和列主元的Gauss消去法求解该方程组，记你计算出的解为̂。对n从5到30估计计算解的精度。(2)对n从5到30计算其逆矩阵。代码：%LU分解functionproblem4()clc%生成矩阵Aforn=5:30A=-ones(n);A(:,n)=ones(n,1);foria=1:nA(ia,ia)=1;forja=(ia+1):(n-1)A(ia,ja)=0;endendx=rand(n,1);b=A*x;%LU分解[n,n]=size(A);L=eye(n);M=eye(n);fori=1:n-1ifA(i,i)==0fprintf('Error:A(%d,%d)=0!\n',i,i);return;endforj=i+1:nL(j,i)=A(j,i)/A(i,i);endfork=i+1:nforl=i:nA(k,l)=A(k,l)-L(k,i)*A(i,l);endendendU=A;M=inv(U)*inv(L)*My1=inv(L)*b;x1=inv(U)*y1;error=norm(x-x1)endreturn;%PLU分解functionproblem4()clc%生成矩阵An=5;A=-ones(n);A(:,n)=ones(n,1);foria=1:nA(ia,ia)=1;forja=(ia+1):(n-1)A(ia,ja)=0;endendx=rand(n,1);b=A*x;%PLU分解[n,n]=size(A);Q=eye(n);R=eye(n);S=eye(n);L=zeros(n,n,n-1);foril=1:(n-1)L(:,:,il)=eye(n);endP=zeros(n,n,n-1);foril=1:(n-1)P(:,:,il)=eye(n);endfori=1:n-1m=i;a=abs(A(i,i));fork=i+1:nifabs(A(k,i))aa=abs(A(k,i));m=k;endendA([i,m],:)=A([m,i],:);P([i,m],:,i)=P([m,i],:,i);forj=i+1:nL(j,i,i)=-A(j,i)/A(i,i);forl=i:nA(j,l)=A(j,l)+L(j,i,i)*A(i,l);endendendU=A;Uforx=1:n-1Q=L(:,:,x);fory=x+1:n-1Q=P(:,:,y)*Q*P(:,:,y);endS=S*inv(Q);end%L矩阵Sforz=1:(n-1)R=P(:,:,z)*R;end%P矩阵Rb1=R*b;b2=inv(S)*b1;x2=inv(U)*b2;M=inv(U)*inv(S)*Rerror=norm(x-x2)结果：由于篇幅有限仅列出n=5时的精度及逆矩阵LU分解M=0.5000-0.2500-0.1250-0.0625-0.062500.5000-0.2500-0.1250-0.1250000.5000-0.2500-0.25000000.5000-0.50000.50000.25000.12500.06250.0625error=8.8818e-16PLU分解M=0.5000-0.2500-0.1250-0.0625-0.062500.5000-0.2500-0.1250-0.1250000.5000-0.2500-0.25000000.5000-0.50000.50000.25000.12500.06250.0625error=7.8516第五题编制计算对称正定阵的Cholesky分解的通用程序，并用你编制的程序计算Ax=b，其中()。b可以由你自己取定，对n从10到20验证程序的可靠性。代码：functionproblem5()%cholesky分解clc%生成矩阵AN=10;A=zeros(N);foria=1:Nforja=1:NA(ia,ja)=1/(ia+ja-1);endendL=chol(A);L结果：由于篇幅有限仅列出n=10时的L矩阵L=1.00000.50000.33330.25000.20000.16670.14290.12500.11110.100000.28870.28870.25980.23090.20620.18560.16840.15400.1417000.07450.11180.12780.13310.13310.13040.12650.12200000.01890.03780.05250.06300.07020.07480.077700000.00480.01190.01950.02650.03260.0378000000.00120.00360.00680.01030.01390000000.00030.00110.00220.003800000000.00010.00030.0007000000000.00000.00010000000000.0000第六题(1)编制程序House(x)，其作用是对输入的向量x，输出单位向量u使得()‖‖。(2)编制Householder变换阵乘以的程序HA，注意，你的程序并不显式的计算出H。(3)考虑矩阵(√√√)用你编制的程序计算H使得HA的第一列为的形式，并将HA的结果显示。代码：functionproblem6()clcx=[1;-1;-2;-10^(1/2);0];ifx(1)0sg=1;elsesg=-1;endNx=length(x);e1=zeros(Nx,1);e1(1)=1;sum=0;fori=1:Nxsum=sum+x(i)*x(i);endx2=sqrt(sum);clearsum;w=x-sg*x2*e1;u=w/sqrt(w'*w);u%生成矩阵AA=[1,2,3,4;-1,3,2^(1/2),3^(1/2);-2,-2,exp(1),pi;-10^(1/2),2,-3,7;0,2,7,5/2];H=eye(Nx)-2*u*u';H*A结果：u=-0.6124-0.2041-0.4082-0.64550ans=4.0000-0.83111.4090-6.53780.00002.05630.8839-1.78050.0000-3.88741.6576-3.88360.0000-0.9843-4.6770-4.107802.00007.00002.5000第七题用Jacobi和Gauss-Seidel迭代求解下面的方程组，输出迭代每一步的误差‖‖:{代码：functionproblem7()clcinter=10;%迭代次数rx=[69*7/61-8,69*3/(61*2)-7/2,69/61]';%根据题意生成A,bA=[5-1-3;-124;-3415];b=[-2110]';N=length(b);x=ones(N,1);%Jacobi迭代法L=zeros(N,N);D=L;U=L;foria=1:ND(ia,ia)=A(ia,ia);ifiaNforiu=(ia+1):NU(ia,iu)=A(ia,iu);endendifia1foril=1:(ia-1)L(ia,il)=-A(ia,il);endendendforintt=1:interx=inv(D)*(L+U)*x+inv(D)*b;wuca1=0;%求误差forix=1:Nwuca1=wuca1+(x(ix)-rx(ix))^2;endwuca1end%G-S迭代法forintt=1:interx=inv(D-L)*U*x+inv(D-L)*b;wuca2=0;%求误差forix=1:Nwuca2=wuca2+(x(ix)-rx(ix))^2;endwuca2end结果:wuca1=24.6036wuca1=12.3259wuca1=1.9276wuca1=5.4520wuca1=15.9787wuca1=10.0748wuca1=3.6731wuca1=6.2770