清华大学高等数值分析 第一次实验作业

10
0
CG法收敛曲线 (阶数 n=1002)
10
-1
10
-2
10
-3
||rk||/||b||
10
-4
10
-5
10
-6
10
-7
10
-8
0
20
40
60
80
100
120
140
160
180
迭代次数
图1 CG法求解良态问题的收敛曲线
（3）CG法求解Ax=b，A病态： 利用matlab编程实现CG算法。b = ones(1002,1)，x0 = zeros(1002, 1)。计算 每一步迭代的残差rk相对于初始残差的2范数。相对残差2范数的对数值与 迭代步数的关系曲线如图2所示：
高等数值分析实验作业一
高等数值分析实验作业一
1. 构造例子说明 CG 的数值性态。当步数 = 阶数时 CG 的解如何？当 A 的 最大特征值远大于第二个最大特征值 , 最小特征值远小于第二个最小特征值时 方法的收敛性如何? 解： （1）构造 1002 阶正定矩阵 A： A的特征值向量为，分为良态和病态情形; （2）CG法求解Ax=b，A良态 : 利用matlab编程实现CG算法。b = ones(1002,1)，x0 = zeros(1002, 1)。计算 每一步迭代的残差rk相对于初始残差的2范数。相对残差2范数的对数值与 迭代步数的关系曲线如图1所示：
10
0
Lanzcos 算法的收敛曲线 (阶数 n=1002)
m=10 m=50 m=100 m=400 m=800
10
-2
10
-4
||rk||/||b||
10
-6
10
-8
10
-10
10
-12
0
2
4
6
8
10
12
14
16
18
20
迭代次数
图6 对于b由m个特征向量线性表示时，Lanczos法求解Ax=b的收敛曲线
10
2
Lanzcos 算法的收敛曲线 (阶数 n=1002）
10
0
10
-2
||rk||/||b||
10
-4
10
-6
10
-8
10
-10
0
50
100
150
200
250
300
350
迭代次数
图7
10
4
m=10时，Lanczos法求解Ax=b的收敛曲线
Lanzcos 算法的收敛曲线 (阶数 n=1002）
10
10
4
Lanzcos 算法的收敛曲线 (阶数 n=1002)
10
2
10
0
||rk||/||b||
10
-2
10
-4
10
-6
10
-8
0
50
100
150
200
250
300
350
400
迭代次数
图4
Lanczos法求解病态问题的收敛曲线
结论：对于良态正定问题， Lanczos法和CG法收敛性一样，相对残差的2范数随 迭代步数的增加而减小；对于病态正定问题 Lanczos法较CG法收敛的要慢一点， 而且收敛曲线更加不平滑，震荡较严重，相对残差的2范数同样没有最优性。 事实上，Lanczos法主要是解对称不定问题，针对这类问题才显现出它的优 越性，这也从侧面表明，对称正定问题不适合用Lanzcos法求解。
Minres 算法的收敛曲线 (阶数 n=1002）
10
0
10
-2
||rk||/||b||
10
-4
10
-6
10
-8
10
-10
0
100
200
300
400
500
600
700
800
900
迭代次数
图14
m=100时，Minres法求解Ax=b的收敛曲线
高等数值分析实验作业一
10
2
Minres 算法的收敛曲线 (阶数 n=1002）
10
-5
||rk||/||b||
10
-10
10
-15
0
20
40
60
80
100
120
140
160
迭代次数
图5 对不同的m，Lanczos法求解Ax=b的收敛曲线
高等数值分析实验作业一
结论：如果 A 只有 m 个不同的特征值，则 Lanczos 方法至多 m 步就可以找到精 确解。实验中，在 m 较大的时候，算法收敛较快，迭代次数远小于 m。当 m 较 小时，可能需要接近于 m 步才能找到准确解。 4. 取初始值近似解为零向量， 右端项 b 仅由 A 的 m 个不同特征向量的线性组合 表示时，Lanczos 方法的收敛性如何？数值计算中方法的收敛性和 m 的大小关 系如何？ 解：设A的特征值向量为：1=[1002:-1:1]，分别构建m =10、50、100、400、800 时的b，即b分别由10、50、100、400、800个A的特征向量线性表示，分别求解 Ax=b，收敛曲线如图6所示：
10
0
Lanzcos 算法的收敛曲线 (阶数 n=1002)
10
-1
10
-2
10
-3
||rk||/||b||
10
-4
10
-5
10
-6
10
-7
10
-8
10
-9
0
20
40
60
80
100
120
140
160
180
迭代次数
图3
Lanczos法求解良态问题的收敛曲线
高等数值分析实验作业一
（3）Lanczos法求解Ax=b，A病态： 利用matlab编程实现CG算法。b = ones(1002,1)，x0 = zeros(1002, 1)。计算 每一步迭代的残差rk相对于初始残差的2范数。相对残差2范数的对数值与 迭代步数的关系曲线如图4所示：
高等数值分析实验作业一
end
Lanczos 法 lanczos.m
function [T,Q,x,k,Errs,tiao]=Lanczos(A,b,x0,Err) x=[]; tiao=[]; [m,n]=size(b); if m<n b=b'; end [m,n]=size(x0); if m<n m=n; x0=x0'; end size(b) size(A*x0) r=b-A*x0; r_zeros=r; q=r/norm(r); q0=0; beta0=0; T=[0]; k=1; Q=[q]; Errs=[norm(r,2)/norm(b,2)]; while 1 disp('k=');disp(k); T=[T,zeros(k,1);zeros(1,k),0]; r=A*q-beta0*q0; T(k,k)=q'*r; r=r-T(k,k)*q; T(k+1,k)=norm(r,2); beta0=T(k+1,k); T(k,k+1)=beta0; q0=q; q=r/beta0; Tk=T(1:k,1:k) ; if k==1 disp('Ìø¹ý´Ë²½') tiao=[tiao,k]; else [L,U]=LanczosLU(Tk); bk=norm(r_zeros,2)*[1;zeros(k-1,1)];
结论：负特征值越多，Lanczos 方法的收敛曲线的峰点个数越多，纹波越大，发 生近似中断的次数越来越多。
高等数值分析实验作业一
10
2
Minres 算法的收敛曲线 (阶数 n=1002）
10
0
10
-2
||rk||/||Βιβλιοθήκη Baidu||
10
-4
10
-6
10
-8
10
-10
0
50
100
150
200
250
300
350
结论：理论上，b 仅由 A 的 m 个不同特征向量的线性组合表示时，Lanczos 方法 必然 m 步收敛。但由于 A 的阶数为 1002，是比较大的，精度方面的限制导致计 算得到的 m 个特征向量并不都线性无关，所以，m = 800 时只需 20 步迭代。
高等数值分析实验作业一
5. 构造对称不定矩阵，验证 Lanczos 方法的近似中断，观察收敛曲线中的峰点 个数和特征值的分布关系； 观测当出现峰点时， MINRES 方法的收敛形态怎样。 解：分别构建负特征值个数为 m =10、50、100、200、500 的矩阵 A，分别计算 Ax = b 的解，Lanzcos 方法的收敛曲线如图 7-图 11 所示，Minres 的收敛曲线如 图 12-图 16 所示：
10
1
10
0
||rk||/||b||
10
-1
10
-2
10
-3
0
200
400
600
800
1000
1200
迭代次数
图15
10
2
m=200时，Minres法求解Ax=b的收敛曲线
Minres 算法的收敛曲线 (阶数 n=1002）
10
1
10
0
||rk||/||b||
10
-1
10
-2
10
-3
0
200
400
600
高等数值分析实验作业一
10
4
CG法收敛曲线 (阶数 n=1002)
10
2
10
0
||rk||/||b||
10
-2
10
-4
10
-6
10
-8
0
50
100
150
200
250
300
350
400
迭代次数
图2 CG法求解病态问题的收敛曲线
结论：对于良态问题，CG法收敛较快，相对残差的2范数随迭代步数的增加而减 小；对于病态问题CG法收敛的很慢，而且收敛曲线不平滑，震荡非常严重，相 对残差的2范数没有最优性。
迭代次数
图12 m=10时，Minres法求解Ax=b的收敛曲线
10
2
Minres 算法的收敛曲线 (阶数 n=1002）
10
0
10
-2
||rk||/||b||
10
-4
10
-6
10
-8
10
-10
0
100
200
300
400
500
600
700
迭代次数
图13
10
2
m=50时，Minres法求解Ax=b的收敛曲线
2
10
0
||rk||/||b||
10
-2
10
-4
10
-6
10
-8
10
-10
0
100
200
300
400
500
600
700
迭代次数
图8
10
2
m=50时，Lanczos法求解Ax=b的收敛曲线
Lanzcos 算法的收敛曲线 (阶数 n=1002）
10
0
10
-2
||rk||/||b||
10
-4
10
-6
10
-8
高等数值分析实验作业一
附件：主要算法代码
CG 法 CG.m
function [x,Error,i,flag]=CG(A,b,x,ErrSet,uplimit) [m,n]=size(b); if m<n b=b'; end [m,n]=size(x); if m<n x=x'; end r=b-A*x; p=r; i=1; temp_rkrkplus=r'*r; Error=sqrt(temp_rkrkplus)/norm(b,2); while 1 temp_AP=A*p; temp_rkrk=temp_rkrkplus; temp_pAP=p'*temp_AP; if abs(temp_pAP)<1e-12 disp('¶ñÐÔÖжϣ¡') break; end a=temp_rkrk/(temp_pAP); x=x+a*p; r=r-a*temp_AP; temp_rkrkplus=r'*r; beta=temp_rkrkplus/temp_rkrk; p=r+beta*p; Err=sqrt(temp_rkrkplus)/norm(b,2); %/(norm(b)+temp_AP); if Err<ErrSet disp('Method converge£¡') disp(i) flag=1; break; end Error=[Error,sqrt(temp_rkrkplus)/norm(b,2)]; if i>uplimit flag=0; break end end end
2. 对于同样的例子，比较 CG 和 Lanczos 的计算结果。 解： （1）构造 1002 阶正定矩阵 A： （2）Lanczos法求解Ax=b，A良态： 利用matlab编程实现CG算法。b = ones(1002,1)，x0 = zeros(1002, 1)。计算 每一步迭代的残差rk相对于初始残差的2范数。相对残差2范数的对数值与 迭代步数的关系曲线如图3所示：
800
1000
1200
迭代次数
图16
m=500时，Minres法求解Ax=b的收敛曲线
结论：对于 Lanczos 方法，随着负特征值的增多，收敛曲线的峰点个数增多，振 荡越来越严重，发生近似中断的次数越来越多；然而，对于相同的 A，Minres 方法的相对残差没有出现峰值， 随着迭代数增加而单调下降。 收敛性态比 Lanczos 好，但代价是计算时间的急剧增大。
400
600
800
1000
1200
迭代次数
图10
10
2
m=200时，Lanczos法求解Ax=b的收敛曲线
Lanzcos 算法的收敛曲线 (阶数 n=1002）
10
0
10
-2
||rk||/||b||
10
-4
10
-6
10
-8
10
-10
0
500
1000
1500
迭代次数
图11
m=500时，Lanczos法求解Ax=b的收敛曲线
10
-10
0
100
200
300
400
500
600
700
800
900
迭代次数
图9
m=100时，Lanczos法求解Ax=b的收敛曲线
高等数值分析实验作业一
10
4
Lanzcos 算法的收敛曲线 (阶数 n=1002）
10
2
10
0
||rk||/||b||
10
-2
10
-4
10
-6
10
-8
10
-10
0
200
3. 当 A 只有 m 个不同特征值时， 对于大的 m 和小的 m， 观察有限精度下 Lanczos 方法如何收敛。 解：分别构建 m = 10、50、100、400、800 五个矩阵 A，分别求解 Ax=b，收敛 曲线如图 5 所示：
10
0
Lanzcos 算法的收敛曲线 (阶数 n=1002)
m=10 m=50 m=100 m=400 m=800