北航研究生数值分析作业第二题：
一、算法设计方案
1. 按照题目给出的矩阵定义对矩阵A赋初值：对应的函数为a_init()； 2. 对矩阵A进行householder变换，使其拟上三角化：对应的函数为
householder()；
3. 输出拟上三角化后的A：对应的函数为aout(int)；
4. 对拟上三角化后的矩阵A使用带双步位移的QR分解法逐次迭代（最大迭代
次数L=500），逐个求出其特征值，对应的函数为eigen_a()；中间包含两个子程序：calc_mk()和qr_analyze()，分别用来计算矩阵Mk和对Mk进行QR分解并得到Ak+1;
5. 输出QR分解过程完毕后的A及求得的特征向量：对应的函数为aout()和
eigenvalout()；
6. 对于在第三步中求得的每个实特征值，使用带原点平移的反幂法求出其对应
的特征向量，对应的函数为eigenvec()；其中包含一个解方程(A-μI)=yk-1的程序段。这部分也用迭代完成，仍然将最大迭代次数L设置为500； 7. 输出矩阵A的特征向量，结束计算：对应的函数为eigenvecout()。 算法编译环境：vlsual c++6.0 二、源程序如下：
#include
#define N 10 //矩阵阶数；
#define EPSL 1.0e-12 //迭代的精度水平； #define L 500 //迭代最大次数；
#define OUTPUTMODE 1 //输出格式：0--输出至屏幕，1--输出至文件
double a[N][N], a2[N][N], eigen[N][N]; //声明矩阵A；
double sa_re[N] = {0}, sa_im[N] = {0}; //声明矩阵的特征值数组； double u_init[N] = {2,1,2,1,2,1,2,1,2,1}; //定义反幂法中使用的初始向量u；
//主程序开始； int main()
{
FILE *p;
void a_init(); void householder();
void equal_zero(double matrix[N][N], int); void eigenvec(); int eigen_a(); void aout(int);
void eigenvalout(int); void eigenvecout(int);
if(OUTPUTMODE) {
p = fopen(\ fprintf(p, \计算结果：\\n\ fclose(p); }
a_init(); //对矩阵A进行初始化；
householder(); //对矩阵A进行拟上三角化；
equal_zero(a, N); //对矩阵A的元素进行归零处理，消除误差； aout(OUTPUTMODE); //输出A；
if(eigen_a()) printf(\迭代超过最大次数，特征值求解结果可能不正确。\\n\//求矩阵A的特征值；
equal_zero(a, N); //对矩阵A的元素进行归零处理，消除误差； aout(OUTPUTMODE); //输出A；
eigenvalout(OUTPUTMODE); //输出矩阵的特征值； eigenvec(); //求矩阵A的特征向量；
eigenvecout(OUTPUTMODE); //输出矩阵A的特征向量； getchar(); }
void a_init() {
int i,j;
for(i = 1; i <= N; i++) {
for(j = 1; j <= N; j++)
a2[i - 1][j - 1] = a[i - 1][j - 1] = sin(0.5 * i + 0.2 * j); }
for(i = 1; i <= N; i++)
a2[i - 1][i - 1] = a[i - 1][i - 1] = cos(i + 1.2 * i) * 1.5; //这里使用if语句反而更慢，所以赋值赋了两次。
}
void
householder()
//对矩阵进行拟上三角化； {
int r,i,j;
double tmp_ir, tmp_dr, tmp_pr, tmp_qr, tmp_tr, dr, cr, hr, tr; double ur[N] = {0}, pr[N] = {0}, qr[N] = {0}, wr[N] = {0}; int sgn2(double);
for(r = 0; r < N-2; r++) {
//第一步；
tmp_ir = 0;
for(i = r + 2; i < N; i++) tmp_ir = tmp_ir + (a[i][r]==0); if (tmp_ir == N-2 - r) continue; else {
//第二步；
tmp_dr = 0;
for(i = r+1; i < N; i++) tmp_dr = tmp_dr + a[i][r] * a[i][r];
dr = sqrt(tmp_dr);
cr = -1.0 * sgn2(a[r + 1][r]) * dr; hr = cr * cr - cr * a[r + 1][r];
//第三步；
for(i = 0; i < r + 1; i++) ur[i] = 0;
for(i = r + 2; i < N; i++) ur[i] = a[i][r]; ur[r + 1] = a[r + 1][r] - cr;
//第四步；
for(i = 0; i < N; i++) {
tmp_pr = 0; tmp_qr = 0;
for(j = 0; j < N; j++) {
tmp_pr = tmp_pr + a[j][i] * ur[j]; tmp_qr = tmp_qr + a[i][j] * ur[j]; }
pr[i] = tmp_pr / hr; qr[i] = tmp_qr / hr; }
tmp_tr = 0;
for(i = 0; i < N; i++) tmp_tr = tmp_tr + pr[i] * ur[i]; tr = tmp_tr / hr;
for(i = 0; i < N; i++) wr[i] = qr[i] - tr * ur[i]; for(i = 0; i < N; i++) {
for(j = 0; j < N; j++) {
a[i][j] = a[i][j] - wr[i] * ur[j] - ur[i] * pr[j];
} }
} //第五步：(继续) } } int
sgn2(double a)
//求cr时用到的sgn子程序 {
if(a >= 0) return 1; else return -1; }
void equal_zero(double
matrix[N][N], int rank)
//对矩阵进行归零处理； {
int i,j;
for(i = 0; i < rank; i++) {
for(j = 0; j < rank; j++) if (-EPSL < a[i][j] && a[i][j] < EPSL) matrix[i][j] = 0; } } int eigen_a() //计算A的特征值； {
int snum = N-1, m = N-1, flag = 0, flag_7to4 = 0, step4_cont = 0, k=1;
double mk[N][N], det_dk, stmpb, s1_re, s1_im, s2_re, s2_im, ms, mt, b24ac;
void calc_mk(double mk[N][N], int, double, double); void qr_analyze(double mk[N][N], int);
while (1)
{
//第三步；
if ((-EPSL < a[m][m - 1] && a[m][m - 1] < EPSL) && !flag_7to4) {
sa_re[snum] = a[m][m], snum--, m--, flag = 1; } else {}
//第四步；
flag_7to4 = 0; if (flag)
{ flag = 0; if (m==0) {
sa_re[snum] = a[0][0]; return 0; }
else if (m==-1) return 0; else step4_cont = 1; }
//第五步； else {
det_dk - 1];
stmpb = -1.0 * (a[m - 1][m - 1] + a[m][m]); b24ac = stmpb * stmpb - 4.0 * det_dk; if(b24ac < 0) {
b24ac = -b24ac;
s1_re = -0.5 * stmpb;
s1_im = 0.5 * sqrt(b24ac); s2_re = -0.5 * stmpb;
s2_im = -0.5 * sqrt(b24ac); } else {
s1_re = (-1.0 * stmpb + sqrt(b24ac)) / 2; s1_im = 0;
s2_re = (-1.0 * stmpb - sqrt(b24ac)) / 2; s2_im = 0; } }
if (step4_cont)
= a[m - 1][m - 1] * a[m][m] - a[m - 1][m] * a[m][m
百度搜索“77cn”或“免费范文网”即可找到本站免费阅读全部范文。收藏本站方便下次阅读，免费范文网，提供经典小说综合文库北航研究生数值分析作业第二题在线全文阅读。
相关推荐：