Matlab数值计算

“高等数值算法与应用”课程作业(2015秋)


课本：《Matlab数值计算》
通过网络学堂提交，包含文档与源程序（若有）。文档答题，可不抄题；源程序为.m文件，注意文件命名规范，例如ex1_3.m.
>> 请将文档与源程序文件放在同一个目录中，不带目录压缩为一个.zip/.rar/.gzip压缩包，命名时包含自己的姓名。谢谢！
>> 对有些作业题，在答题文档中最好用贴图的方式直观说明所编程序的运行效果。谢谢！




第2次作业（deadline: 2015.10.29）


一. 证明: 非零元分布对称的稀疏矩阵A，对它做不选主元的LU分解，填入元的位置也是对称的。


二. 下面两个证明题中任选一个做
1. 在共轭梯度法中（见第3讲课件pp. 28~29），假设循环体执行了k次，得到近似解序列{x_0,⋯,x_k }, 搜索方向序列{p_0,⋯,p_k }，残差序列{r_0,⋯,r_k }，且{p_0,⋯,p_(k-1) }这k个向量线性无关。试证明（不考虑舍入误差）：
(1) r_i∈K^k (A,b), p_i∈K^k (A,b), 其中0≤i≤k-1，且x_k∈K^k (A,b).
(2) x_k是K^k （A,b）空间中使φ(x)达到最小的那个x，即设x=∑_(i=0)^(k-1)▒α_i  p_i，证明只有当{α_i}的取值如课件所述时φ(x)取最小值。
提示：采用数学归纳法证明。
(3) r_k⊥K^k (A,b)，即r_k⊥{p_0,⋯,p_(k-1) }这k个向量.
提示：反证法。假设r_k^T p_i=θ≠0, 考虑x_k^'=x_k+cp_i，可证明取适当的c值可以使φ(x_k^' )<φ(x_k )。
(4) 若p_k与{p_0,⋯,p_(k-1) }线性相关，则一定有r_k为零向量（发生良性中断）。
(5) 搜索方向满足A正交，即p_i^T Ap_k=0,(i<k).
提示：数学归纳法证明，利用公式p_k=r_k+β_(k-1) p_(k-1)，以及r_k的性质。


2. 见第10章习题10.4，证明该矩阵A的特征值公式为(n为矩阵A的阶数):
λ_i=-2 cos⁡(iπ/(n+1)),i=1,2,…,n.


三．第5章习题5.5（提示：从Householder矩阵的几何意义考虑，(b)问就不难了）。


四. 矩阵A=I-uv^T称为初等阵（elementary matrix），其中I为单位阵，u,v为两个向量。
(1)证明：当v^T u≠1时，初等阵A非奇异。
(2) LU分解中的消去阵M、行交换阵P（见第2讲课件）、Householder变换矩阵H、Givens变换矩阵G都是初等阵吗？若是，请说明原因；若不是，又为什么？


五. 编程实现基于Givens旋转的稀疏矩阵QR分解算法。测试下面两个矩阵：
A= sprand(50, 40, 0.1);
B= sprand(1000, 800, 0.01);
测试程序正确性，即得到的Q满足正交矩阵要求、R是上三角阵、A=QR.
观察填入现象，即看R中非零元数目比原始A中的多了多少。
实现一种重排序的QR分解，函数接口为[Q, R, E]= givensQR(A), 其中AE=QR, E为排列阵。排序策略是在未消元的列中找非零元数目最少的，将其与当前列交换，然后用Givens旋转消元。测试算法正确性，并看看是否缓解了R矩阵的填入问题。
提示：用find命令确定非零元位置，例如find(A(:,1))返回A矩阵第1列哪些行有非零元。


六. 第10章习题10.2(b).
提示：我们在第一堂课上已经讨论过10.2(a)，以此为启发进行思考。


七. 见第10章习题10.9，研究三种算矩阵奇异值的方法
(1) 说明它们为什么都能算矩阵A的奇异值。
(2) 通过随机矩阵A, 测试三种算法的计算速度区别，并尝试分析原因。（注意：用阶数足够大的矩阵测试，tic/toc命令测出的时间也有意义。）
(3) 从准确度角度看，上述三种方法的优劣如何？


八. 基于qrtrid.m程序，编写计算实对称矩阵特征值的QR算法[D, it]= myqr(A)，其中D为包含特征值的一维数组，it为QR迭代次数（参考第5讲课件pp.19）。
(1) 考虑稠密矩阵A，实现基于Householder变换的正交三对角化；
(2) 修改qrtrid函数，不要使用qr命令（程序第10行），而采用Givens旋转变换实现一步QR迭代；
(3) 在qrtrid函数中，实现隐式位移技术。
(4) 用命令A= rand(n); A=A+A’; 得到的矩阵进行测试，并与eig命令的结果比较。