数值分析大作业1

0.1 i
, b 0.16, c 0.064 ， 矩阵 A 的特征值
i (i 1, 2, , 501) 满足 1 2 …… 501，|s|=min|i|；
试求： （1） 1 ， 501 和 s 的值； （2）A 的与数 k 1 k
501 1
第三章 全部源程序
// shuzhi1.cpp : 定义控制台应用程序的入口点。
6
北京航空航天大学《数值分析 A》计算实习报告
// #include "stdafx.h" #include "math.h" #define e 1.0e-12 #define N 501 const double b=0.16; const double c=-0.064; //设定精度水平
2
北京航空航天大学《数值分析 A》计算实习报告
第二章 算法的设计方案 2.1 题目分析
由题目所求问题出发，结合题目已知条件方阵 A 为实对称矩阵，并且
1 2 501 ，根据已有知识：对非奇异的实对称矩阵 A 有 cond ( A) 1 2 n
（ 1 和 n 分别是矩阵 A 的模为最大和模为最小的特征值） ； 由定义可知，S 相当 于为矩阵 A 的按模最小的特征值； 数 k 1 k
第三章 全部源程序 ..................................................................................................... 6 第四章 特征值、条件数和行列式值的输出 .............................................. 15
1
'
k ；
1 （ 1 和 n 分别是 n
矩阵 A 的模为最大和模为最小的特征值） ；由已经求解出的 1 ， 501 和 S 的 数值（ 1 10.7001 ， 501 9.7246 ， S 0.0055579 ； ） ，可以判断得出按模 最大和按模最小的特征值分别为 1 10.7001 ， S 0.0055579 将其带入
yk 1 uk 1 / k 1
uk Ayk 1
T k yk 1u k
K=1,2，…
k k1 / k
Y
N
出口
图1
幂法程序流程图
3）反幂法的算法说明和程序流程图 （1）反幂法算法说明 反幂法的流程与幂法的迭代流程基本相同，但是在反幂法的迭代中， k 的求解由幂法中的 uk Ayk 1 变为通过解线性方程 Ak yk 1 得到，并且迭代
北京航空航天大学《数值分析 A》计算实习报告
第一章 概述
1.1《数值分析 A》计算实习概述
1）目的： 训练运用计算机进行科学与工程计算的能力。 2）要求： （1）独立进行算法设计、程序设计和上机运算，并得出正确的结果； （2）编制程序时全部采用双精度，要求按题目的要求设计输出，并执行 打印。
1.2 题目概述
5）矩阵 A 中与数 k 1 k 求解 根据已经求解得出 1 和 501 的数值，代入式 k 1 k
501 1
40
, k 1, 2, ,39 最接近的特征值 ik 的
501 1
40
中，利用
VC++ 中 的 for 循 环 语 句 求 解 出 k , k 1, 2, ,39 的 数 值 ， 矩 阵
//定义常量b、c
/************************************************************************/ /* 函数： 输入稀疏矩阵A的主对角线元素， 将A存储于二维数组C中(为了节约存储空间)*/ /************************************************************************/ void Structure(double Diagonal_A[] , double Store_C[][N] ) { int i; Store_C[0][0]=Store_C[0][1]=Store_C[1][0]=Store_C[4][N-1]=Store_C[4][N-2]=St ore_C[3][N-1]=0; for (i=2;i<N;i++) { Store_C[0][i]=c; } for (i=1;i<N;i++) { Store_C[1][i]=b; } for (i=0;i<N;i++) { Store_C[2][i]=Diagonal_A[i]; } for (i=0;i<N-1;i++) { Store_C[3][i]=b; } for (i=0;i<N-2;i++) { Store_C[4][i]=c; } }
4.1 结果输出 ............................................................................................................. 15
第五章 关于迭代初始向量的选取讨论 ......................................................... 16
cond ( A)2
1 10.7001 1925.204 ，即为所求矩阵 A 的条件数。行列式 n 0.0055579
detA 的求解，对矩阵 A 先进行 LU 分解，然后利用 VC++的 for 循环语句， 求取 U 的对角线元素乘积， 所得数值即为即为矩阵 A 的行列式值。 求解得出 的 detA=2.772786141752E+118；
501 ，数值小的即为 1 ；
5
北京航空航天大学《数值分析 A》计算实习报告
Ⅳ S 的求解：
S 为矩阵 A 的按模最小的特征值，采用反幂法求解。
求解结果：本题利用上述求解步骤调用幂法函数和反幂法函数求解得出
1 ， 501 和的数值分别为
1 10.700113615 ， 501 9.7246341， S 0.0055579 ；
501 三角分解法对矩阵 A 进行 LU 分解， 则行列式值 det A 1 矩 U 主对角线元素；
阵 A 中与 k 最接近的特征值 ik 可以用反幂法求解矩阵 （ A k I ） 按模最小的特 征值 '' ，则 ik
1
''
uk 。
2.2 算法设计
1）方阵 A 的存储 方阵 A 为 501 阶实对称方阵， 且上下半带宽均为 2， 属于大型带状矩阵， 为了节省存储量，方阵 A 的带外元素不给存储，仅存 A 的带内元素。为此， 设置一个二维数组 Store_C（m ，n） ，用于存放方阵 A 的带内元素，其中 m=2+2+1=5，n=501。数组 C 的第 j 列存放方阵 A 的第 j 列带内元素，并使 方阵 A 的主对角线元素存放在数组 C 的第 3 行中。数组 C 存放完方阵 A 的 带内元素后，多出的单元取零，这些零元共有（1+2）+（1+2）=6 个。在数 组 C 中检索方阵 A 的带内元素的方法是： A 的带内元素 aii =C 中的元素 ci j 2, j ；
1.1《数值分析 A》计算实习概述 .......................................................................... 1 1.2 题目概述 ................................................................................................................ 1
第二章 算法的设计方案 ........................................................................................... 3
2.1 题目分析 ................................................................................................................ 3 2.2 算法设计 ................................................................................................................ 3
1 0 5 1
40Байду номын сангаас
直接以 1 和 501 为系数因
子； 可以发现，1 和 501 的求解是本题的求解关键， 可以选择幂法求解 1 和 501 的 数值，用反幂法求解方阵 A 的按模最小的特征值 S 。由于反幂法中要求解线性 方程 Auk yk 1 并且题目要求求解行列式的值 detA，考虑到精度要求问题，选择
4
北京航空航天大学《数值分析 A》计算实习报告
k1 k11 k k 1 精度判断由幂法的 ，最后得到的特征值取倒 变为 k k1
数即为所要求的按模最小的特征值的数值。 （2）反幂法的程序流程图
图2 4） 1 、 501 和 S 的求解
反幂法程序流程图
5.1 迭代初始向量的选取对结果的影响 ............................................................ 16 5.2 讨论........................................................................................................................ 16
《数值分析 A》计算实习 题目：
幂法和反幂法求矩阵特征值
院（系）名 专 班 学 姓 业 名
称 称 级 号 名
机械工程及自动化学院 航空宇航制造工程 SY12073 班 SY1207415 陈俊
2012 年 10 月
目录 第一章 概述 ..................................................................................................................... 1
已知：设有 501×501 的矩阵
b c a1 b a2 b c b a3 A501501＝ …… c
其中：
c b b c
c a499 b c b a500 b
c b a501
ai (1.64 0.024i)sin(0.2i) 0.64e
40
最接近的特征值 ik (k 1,2 39) ；
（3）A 的（谱范数）条件数 cond(A)2 和行列式 detA； 说明： （1）在所有的算法中，凡是要给出精度水平的ξ，都取 1012 ；
1
北京航空航天大学《数值分析 A》计算实习报告
（2）选择算法时，应使 A 的所有零元素都不存储； （3）打印以下内容： ·算法的设计方案； ·全部源程序（要求注明主程序和每个子程序的功能） ·特征值和（k=1,2,…,39） ； ·讨论迭代初始向量的选取对计算结果的影响，并说明原因； （4）采用 e 型输出所有计算结果，并且至少显示 12 为有效数字。
3
北京航空航天大学《数值分析 A》计算实习报告
2）幂法的算法说明及程序流程图 （1）幂法的算法说明 A. 取初始向量u[i] B. 初始向量的2——范数 C. 将初始向量单位化，得到初始迭代向量 D. 算法迭代 E. 求出特征值 （2）幂法程序流程图
入口
任选初始迭代向量U 0
k 1 u(Tk 1)uk 1
因为方阵 A 为实对称矩阵，并且 1 2 501 ，所以矩阵 A 的按模最 大的特征值必为 1 或 501 ，采用如下步骤求解 1 和 501 ： Ⅰ 采用幂法求解矩阵 A 中按模最大的特征值 M 1 ； Ⅱ 矩阵 B A M I ，采用幂法求解矩阵 B 中按模最大的特征值 M 2 ； Ⅲ M 3 M 1 M 2 ；比较 M 1 和 M 3 二者的数值大小，则数值大的即为
S t o r e _A D k I ，即有 Store_D[2][ i ] Store_D[2][ i ]k ，对其采用反幂法
求解按模最小的特征值 ' ，则与 k 最接近的对应的特征值 ik 6）矩阵 A 的条件数 cond(A)2 和行列式 detA 的求解 因为 A 是非奇异的实对称矩阵，所以有 cond ( A) 2