北航数值分析报告大作业一.docx

北航数值分析报告大作业一.docx

《北航数值分析报告大作业一.docx》由会员分享，可在线阅读，更多相关《北航数值分析报告大作业一.docx（14页珍藏版）》请在冰豆网上搜索。

北航数值分析报告大作业一

北京航空航天大学

数值分析大作业一

学院名称自动化

专业方向控制工程

学号ZY*******

学生姓名许阳

教师孙玉泉

日期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

#include

doubleep=1e-12,b=0.16,c=-0.064;

intj=0;

doublepower（doublea[501]）;//幂法

doubleinv_power（doublea[501]）;//反幂法

doubledet（doublea[501]）;//求det

intmain（）//主程序

{

inti,k;

doubleA[501],B[501],beta_1,beta_501,beta_s,beta_k;

doublemu;

for（i=0;i<501;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;i<501;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;i<501;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））;

}

doublepower（doublea[501]）//幂法

{

inti=0,N=5000;

doubleb=0.16,c=-0.064;

doubleu[501],y[501];

doublem=1,beta;

for（i=0;i<501;i++）

u[i]=1;

j=0;

while（j

{

for（i=0;i<501;i++）

{

y[i]=u[i]/fabs（m）;

}

u[0]=a[0]*y[0]+b*y[1]+c*y[2];

u[1]=b*y[0]+a[1]*y[1]+b*y[2]+c*y[3];

u[499]=c*y[497]+b*y[498]+a[499]*y[499]+b*y[500];

u[500]=c*y[498]+b*y[499]+a[500]*y[500];

for（i=2;i<499;i++）

{u[i]=c*y[i-2]+b*y[i-1]+a[i]*y[i]+b*y[i+1]+c*y[i+2];}

beta=0;

for（i=0;i<501;i++）

{

if（fabs（u[i]）>=fabs（beta））

beta=u[i];

}

if（beta<0）

if（fabs（fabs（beta）-fabs（m））/fabs（beta）

break;

if（fabs（beta-m）/fabs（beta）

break;

m=beta;j++;

}

returnbeta;

}

doubleinv_power（doublea[501]）//反幂法

{

doublep[501],r[501],t[501],q[501],u[501],y[501];

doublebeta,m=1;

inti,N=1000;

p[0]=a[0];t[0]=b/p[0];r[1]=b;

p[1]=a[1]-r[1]*t[0];

q[0]=c/p[0];q[1]=c/p[1];

t[1]=（b-r[1]*q[0]）/p[1];

for（i=2;i<501;i++）

{

r[i]=b-c*t[i-2];

p[i]=a[i]-c*q[i-2]-r[i]*t[i-1];

q[i]=c/p[i];

t[i]=（b-r[i]*q[i-1]）/p[i];

}

for（i=0;i<501;i++）

u[i]=1;

j=0;

while（j

{

for（i=0;i<501;i++）

{

y[i]=u[i]/fabs（m）;

}

u[0]=y[0]/p[0];

u[1]=（y[1]-r[1]*u[0]）/p[1];

for（i=2;i<501;i++）

u[i]=（y[i]-c*u[i-2]-r[i]*u[i-1]）/p[i];

u[499]=u[499]-t[499]*u[500];

for（i=498;i>=0;i--）

u[i]=u[i]-t[i]*u[i+1]-q[i]*u[i+2];

beta=0;

for（i=0;i<501;i++）

{

if（fabs（u[i]）>=fabs（beta））

beta=u[i];

}

if（beta<0）

if（fabs（fabs（beta）-fabs（m））/fabs（beta）

break;

if（fabs（beta-m）/fabs（beta）

break;

m=beta;j++;

}

return1/beta;

}

doubledet（doublea[501]）//求det

{

doubledet_A=1;

doublep[501],r[501],t[501],q[501];

inti;

p[0]=a[0];t[0]=b/p[0];r[1]=b;

p[1]=a[1]-r[1]*t[0];

q[0]=c/p[0];q[1]=c/p[1];

t[1]=（b-r[1]*q[0]）/p[1];

for（i=2;i<501;i++）

{

r[i]=b-c*t[i-2];

p[i]=a[i]-c*q[i-2]-r[i]*t[i-1];

q[i]=c/p[i];

t[i]=（b-r[i]*q[i-1]）/p[i];

}

for（i=0;i<501;i++）

det_A=det_A*p[i];

returndet_A;

}

四程序结果

五计算过程中的现象

使用

作为终止迭代条件时，出现迭代无法终止的情况，通过调试发现按模最大特征值为负时，当k充分大后，迭代向量

各分量不断变号，使得

与

异号，判别式

不收敛。

因此将终止迭代条件修改为

，程序实现如下：

if（beta<0）

if（fabs（fabs（beta）-fabs（m））/fabs（beta）

break;

if（fabs（beta-m）/fabs（beta）

break;

从迭代次数可以看出

与

收敛较慢，由按模最大特征值与按模次大特征值的比值越小，收敛速度越慢，可知存在与

和

的模相近的特征值。