北京航空航天大学 2014年研究生
《数值分析》作业
第二题
院系：XXXXXXXXX学院 学号：XXXXXXXXX 姓名：XXXX
Email:XXXXXXXXXXXXXXX
2014年11月
算法设计方案
1、矩阵A的输入与存储
根据题目所给的条件输入需要求解的10阶矩阵A，即为源程序中的void assignment()子程序。
2、将矩阵A拟上三角化
利用矩阵的拟上三角化的算法把矩阵A?Rn?n化为拟上三角矩阵A(n?1)，由子程序void Hessenberg()完成。
矩阵A的拟上三角化算法：
(r)a记A(1)?A，并记A(r)??ij???n?n。对r=1，2，?，n-2执行
（1）若ai(,rr)(i?r?2,r?3,?,n)全为零，则令A(r?1)?A(r)，转（5）；否则转（2）。
（2）计算dr?i?r?1?(an(r)2ir)(r)则取cr?dr），)，cr??sign(ar(r?1,r)dr（若ar?1,r?0，
hr?cr2?cra(r)r?1,r
)(r)(r)Tn（3）令ur?(0,?,0,ar(r?1,r?cr,ar?2,r,?,an,r)?R
T（4）计算pr?A(r)Turhr，qr?A(r)urhr，tr?prurhr，?r?qr?trur，TT A(r?1)?A(r)??rur?urpr（5）继续
3、对拟上三角化后的矩阵A(n?1)进行QR分解
为了了解普通的QR分解过程及结果，下述程序中用void QRfenjie()子程序来对拟上三角化过后的A阵进行QR分解，并输出Q阵R阵RQ阵。其中QRfenjie()既用于基本的QR分解和输出Q阵、R阵、RQ阵，又用于带双步位移的QR方法。m为维数，m<=10。当n=0时，子程序进行基本QR分解；n≠0时，子程序进行带双步位移的QR分解。
- 1 -
4、对拟上三角化后的矩阵A(n?1)进行带双步位移的QR分解。
子程序void QRfa()实现对拟上三角化后的A阵进行带双步位移的QR分解，得出特征值并输出，并用子程序void vector()对其中的实数特征值进行求解，得出对应的特征向量。
带双步位移的QR方法具体算法如下：
（1）使用矩阵的拟上三角化的算法把矩阵A?Rn?n化为拟上三角矩阵A(n?1)；给定精度水平??0和迭代最大次数L。
(1)（2）记A1?A(n?1)???aij??n?n，令k=1，m=n。
(k)(k)（3）如果am,m?1??，则得到A的一个特征值am,m，置m=n-1，转（4）；
否则转（5）。
(k)（4）如果m=1，则得到A的一个特征值a1,1，转（11）；如果m=0，则直
接转（11）；如果m?1，则转（3）。
(k)(k)?am?a?1,m?1m?1,m（5）求二阶子阵Dk??(k)?的两个特征值s1和s2，即计算二次方程(k)??am,m?1am,m??(k)(k)?2?(am?1,m?1?am,m)??detDk?0的两个根s1和s2。
（6）如果m=2，则得到A的两个特征值s1和s2，转（11）；否则转（7）。
(k)（7）如果am则得到A的两个特征值s1和s2，置m=m-2，转（4）；?1,m?2??，
否则转（8）。
（8）如果k=L，则计算终止，未得到A的全部特征值，否则转（9）。
(k)a（9）记Ak???ij??m?m(1?i,j?m)，计算：
(k)(k)s?am?a?1,m?1m,m
(k)(k)(k)(k)t?am?1,m?1am,m?am,m?1am?1,m
Mk?Ak2?sAk?tI（I为m阶单位矩阵） Mk?QkRk（对Mk作QR分解）
TAk?1?QkAkQk
- 2 -
（10）置k=k+1，转（3）。
（11）A的全部特征值以计算完毕，停止计算。
注：对Mk作QR分解与Ak?1的计算算法为：
(1)(r)??bB?b记B1?Mk??，ijrij??m?m???m?m，C1?Ak。对于r=1，2，?，m-1执行
)（1）若bi(,r,r?2,?,m)全为零，则令Br?1?Br,Cr?1?Cr，转（5）；r(i?r?1否则转（2）。
（2）计算dr??(bi?rm(r)2ir)，cr??sign(br(,rr))dr（若br(,rr)?0，则取cr?dr），
hr?cr2?crbr(,rr)
r)(r)Tm（3）令ur?(0,?,0,br(,rr)?cr,br(?1,r,?,bm,r)?R
T（4）计算vr?BrTurhr，Br?1?Br?urvr，pr?CrTurhr，qr?Crurhr，TTT tr?prurhr，?r?qr?trur，Cr?1?Cr??rur?urpr（5）继续
- 3 -
源程序
#include
#define N 10 /*定义矩阵的阶数*/ #define E 1.0e-12 /*定义全局变量相对误差限*/ #define L 20 /*迭代次数*/ FILE *fp;
struct C /*定义结构体*/ { };
void shuchu(double a[N][N]) /*输出一个10*10的矩阵*/ {
int i,j; double b=0; for(i=0;i for(j=0;j if(fabs(a[i][j]) fprintf(fp,\ - 4 - double R; double I;