贾哥高等数值分析第一次实验

=10、50、100、1000四个矩阵A,设置停机准则为:norm(rk)<sqrt(eps),分别计算Ax=b的解,如图6所示。图6各个m值达到收敛时迭代步数如表2所示。表2m5008009509901000迭代步数162856119242可见,相同的特征值个数越多,Lanczos收敛速度越快,而且只要m值稍小于矩阵阶数,收敛速度就明显提高,如图中m=990比m=1000的收敛速度快了一倍。事实上,由于所构造的矩阵条件数都不坏(最大为1000),因此算法的收敛步数远小于m值。附matlab程序3:clearallclcfori=1:5n=1000;M=[500,800,950,990,1000];m=M(i);%m=500、800、950、990、1000D=diag([(1001-m)*ones(1,n-m),linspace(1001-m,1000,m)]);w=eye(1,n);I=eye(n);H=I-w*w';A=H'*D*H;%生成1000阶的对称矩阵b=rand(n,1);X0=zeros(n,1);%初始解x=X0;r0=b-A*X0;R0=r0;%======Lanczos解法========ticy=norm(r0);F(1)=norm(r0);Q(:,1)=r0/norm(r0);r0=A*Q(:,1);alpha(1)=Q(:,1)'*r0;r0=r0-alpha(1)*Q(:,1);bita(1)=norm(r0);Q(:,2)=r0/bita(1);%给各变量赋初始值forj=2:nr0=A*Q(:,j)-bita(j-1)*Q(:,j-1);alpha(j)=Q(:,j)'*r0;r0=r0-alpha(j)*Q(:,j);ifnorm(r0)>1bita(j)=norm(r0);Q(:,j+1)=r0/bita(j);endfork=1:j