北航数值分析大作业一.doc
文本预览下载声明
北 京 航 空 航 天 大 学
数值分析大作业一
学院名称 自动化
专业方向 控制工程
学 号 ZY1403140
学生姓名 许阳
教 师 孙玉泉
日 期 2014 年 11月26 日
设有的实对称矩阵A，
其中，。矩阵A的特征值为，并且有
1.求，和的值。
2.求A的与数最接近的特征值。
3.求A的(谱范数)条件数和行列式detA。
一 方案设计
1 求，和的值。
为按模最小特征值，。可使用反幂法求得。
，分别为最大特征值及最小特征值。可使用幂法求出按模最大特征值，如结果为正，即为，结果为负，则为。使用位移的方式求得另一特征值即可。
2 求A的与数最接近的特征值。
题目可看成求以为偏移量后，按模最小的特征值。即以为偏移量做位移，使用反幂法求出按模最小特征值后，加上，即为所求。
3 求A的(谱范数)条件数和行列式detA。
矩阵A为非奇异对称矩阵，可知，
(1-1)
其中为按模最大特征值，为按模最小特征值。
detA可由LU分解得到。因LU均为三角阵，则其主对角线乘积即为A的行列式。
二 算法实现
1 幂法
使用如下迭代格式：
(2-1)
终止迭代的控制理论使用，
实际使用
(2-2)
由于不保存A矩阵中的零元素，只保存主对角元素a[501]及b,c值。则上式中简化为：
(2-3)
2 反幂法
使用如下迭代格式：
(2-4)
其中，解方程求出。
求解过程中使用LU分解，由于A为5对角矩阵，选择追赶法求取LU分解。求解过程如下：
追赶法求LU分解的实现：
(2-5)
由上式推出分解公式如下：
(2-6)
推导出回代求解公式如下：
(2-7)
(2-8)
3 及A行列式求解
(2-9)
由式(2-5)可得：
(2-10)
三 源程序
#include stdio.h
#include math.h
double ep=1e-12,b=0.16,c=-0.064;
int j=0;
double power(double a[501]); //幂法
double inv_power(double a[501]); //反幂法
double det(double a[501]) ; //求det
int main() //主程序
{
int i,k;
double A[501],B[501],beta_1,beta_501,beta_s,beta_k;
double mu;
for(i=0;i501;i++)
A[i]=(1.64-0.024*(i+1))*sin(0.2*(i+1))-0.64*exp(0.1/(i+1));
beta_1=power(A); //第一问
printf(λ1\t= %.12e\t迭代次数：%d\n,beta_1,j);
for(i=0;i501;i++) //位移
B[i]=A[i]-beta_1;
beta_501=power(B)+beta_1;
printf(λ501\t= %.12e\t迭代次数：%d\n,beta_501,j);
beta_s=inv_power(A);
printf(λs\t= %.12e\t迭代次数：%d\n,beta_s,j);
for(k=1;k=39;k++) //第二问
{
mu=beta_1+k*(beta_501-beta_1)/40;
for(i=0;i501;i++)
B[i]=A[i]-mu;
beta_k=inv_power(B)+mu;
printf(λi%d\t= %.12e\t迭代次数：%d\n,k,beta_k,j);
}
printf(cond(A)2= %.12e\n,beta_1/beta_s); //第三问
printf(detA\t= %.12e\n,det(A));
}
double power(double a[501]) //幂法
{
int i=0,N=5000;
double b=0.16,c=-0.064;
double u[501],y[501];
double m=1,beta;
for(i=0;i501;i++)
显示全部