贾哥高等数值第一次实验解析.docVIP

高等数值分析第一次实验 第一题： 构造例子说明CG的数值形态。当步数 = 阶数时CG的解如何？当A的最大特征值远大于第二个最大特征值，最小特征值远小于第二个最小特征值时，方法的收敛性如何？ 解： 用Housholder变换和对角阵构造1000阶正定对称矩阵A： 构造对角阵D = diag( linspace(1, 1000, 1000) )； 构造Householder阵H。取单位向量w=[1,0,0,.....0]T，I为1000阶单位矩阵，H = I – wTw。 构造对称正定矩阵A。A = HTDH。由于D是对角阵，H是对称的，所以A对称；且A与D具有相同的特征值linspace(1, 1000, 1000) 0，因此A对阵正定。 b=rand(1000,1);取初始解x0=zeros(1000,1); 1.计算Ax = b 利用matlab编程实现CG算法。由于实际计算存在机器误差，因此迭代1000步后的残差不等于0，因此不能用rk=0作为停机准则，否则matlab会无休止地计算下去。本例采用停机准则为：迭代步数=1000步。 当D = diag( linspace(1, 1000, 1000) )时，条件数k=1000； 当D = diag( linspace(1, 100, 1000) )时，条件数k=100； 当D = diag( linspace(1, 10, 1000) )时，条件数k=10； 分别计算上述三种条件数下的CG算法，得到迭代步数与残差的曲线图。 图1：log(‖rk‖)与步数关系曲线。横坐标是迭代步数，纵坐标是残差的对数值。 图 1 clc %条件数k=1000 D=diag(linspace(1,1000,1000)); w=eye(1,1000); I=eye(1000); H=I-w*w; A=H*D*H;%生成1000阶的对称矩阵 b=rand(1000,1); x=zeros(1000,1);%初始解 r=b-A*x;%初始残量 p=r;%初始搜索方向 k=0; semilogy(0,norm(r),r-); hold on; while k1000 alpha = r*p/(p*A*p); x = x+alpha*p; rold = r; r = rold-alpha*A*p; beta = r*r/(rold*rold); p = r+beta*p; semilogy(k,norm(r),r-); hold on; k=k+1; end %条件数k=100 clear all D=diag(linspace(1,100,1000)); w=eye(1,1000); I=eye(1000); H=I-w*w; A=H*D*H;%生成1000阶的对称矩阵 b=rand(1000,1); x=zeros(1000,1);%初始解 r=b-A*x;%初始残量 p=r;%初始搜索方向 k=0; semilogy(0,norm(r),b-); hold on; while k1000 alpha = r*p/(p*A*p); x = x+alpha*p; rold = r; r = rold-alpha*A*p; beta = r*r/(rold*rold); p = r+beta*p; semilogy(k,norm(r),b-); hold on; k=k+1; end %条件数k=10 clear all D=diag(linspace(1,10,1000)); w=eye(1,1000); I=eye(1000); H=I-w*w; A=H*D*H;%生成1000阶的对称矩阵 b=rand(1000,1); x=zeros(1000,1);%初始解 r=b-A*x;%初始残量 p=r;%初始搜索方向 k=0; semilogy(0,norm(r),black-); hold on; while k1000 alpha = r*p/(p*A*p); x = x+alpha*p; rold = r; r = rold-alpha*A*p; beta = r*r/(rold*rold); p = r+beta*p; semilogy(k,norm(r),black-); hold on; k=k+1; end title(条件数的大小对CG法收敛特性的影响); xlabel(迭代步数) ylabel(残差对数log(||rk||)) 2.构造特殊特征值分布 构造对称