北航数值分析大作业一

北 京 航 空 航 天 大 学

数值分析大作业一

学院名称 自动化 专业方向 控制工程 学 号 ZY1403140 学生姓名 许阳 教

师 孙玉泉

日 期 2014 年 11月26 日

设有501?501的实对称矩阵A，

?a1bc??b??????A??c???c?

?????b???cba501???其中，ai?(1.64?0.024i)sin(0.2i)?0.64e(i?1,2,???,501),b?0.16,c??0.064。矩阵A的特征值为?i(i?1,2,???,501)，并且有

0.1i?1??2??????501,|?s|?min|?i|

1?i?5011.求?1，?501和?s的值。 2.求A的与数?k??1?k?501??140最接近的特征值?ik(k?1,2,???,39)。

3.求A的(谱范数)条件数cond(A)2和行列式detA。 一 方案设计

1 求?1，?501和?s的值。

?s为按模最小特征值，|?s|?min|?i|。可使用反幂法求得。

1?i?501 ?1，?501分别为最大特征值及最小特征值。可使用幂法求出按模最大特征值，如结果为正，即为?501，结果为负，则为?1。使用位移的方式求得另一特征值即可。

2 求A的与数?k??1?k?501??140最接近的特征值?ik(k?1,2,...,39)。

题目可看成求以?k为偏移量后，按模最小的特征值。即以?k为偏移量做位移，使用反幂法求出按模最小特征值后，加上?k，即为所求。

3 求A的(谱范数)条件数cond(A)2和行列式detA。 矩阵A为非奇异对称矩阵，可知，

cond(A)2?|?max|?min

(1-1)

其中?max为按模最大特征值，?min为按模最小特征值。

detA可由LU分解得到。因LU均为三角阵，则其主对角线乘积即为A的行列式。

二 算法实现

1 幂法

使用如下迭代格式：

(0)T?任取非零向量u0?(u1(0),???,un)??yk?1?uk?1/max|uk?1| ?u?Ayk?1?k???k?sgn(max|uk?1|)?max|uk?1| (2-1)

终止迭代的控制理论使用|?k??k?1|/|?k|??， 实际使用

||?k|?|?k?1||/|?k|??

(2-2)

由于不保存A矩阵中的零元素，只保存主对角元素a[501]及b,c值。则上式中uk?Ayk?1简化为：

?u(1)?a(1)y(1)?by(2)?cy(3)?u(2)?by(1)?a(2)y(2)?by(3)?cy(4)??)?u(500)?cy(498)?by(499)?a(500)y(500)?by(501 ?)?cy(499)?by(500)?a(501)y(501)?u(501?u(i)?cy(i?2)?by(i?1)?a(i)y(i)?by(i?1)?cy(i?2)???(i?3,???,499)(2-3)

2 反幂法

使用如下迭代格式：

(0)T?任取非零向量u0?(h1(0),???,hn)??yk?1?uk?1/max|uk?1| ?-1u?Ayk?1?k???k?sgn(max|uk?1|)?max|uk?1| (2-4)

其中uk?A?1yk?1?Auk?yk?1，解方程求出uk。

求解过程中使用LU分解，由于A为5对角矩阵，选择追赶法求取LU分解。求解过程如下：

LUuk?yk?1?Lxk?yk?1?Uuk?xk?uk

追赶法求LU分解的实现：

??a1bc?b????A?????c???c??LU????b????cba501????p1??r2???????1t

1q1??????z3????????q499???????????t?500??z501r501p501????1??

由上式推出分解公式如下：

??p1?a1,t1?b/a1?r2?b,p2?a2?r2t1??qi?c/pi,i?1,...,499?ti?(b?riqi?1)/pi,i?2,...,500 ??zi?c,i?3,...,501??ri?b?cti?2,i?3,...,501?pi?ai?cqi?2?riti?1,i?3,...,501推导出回代求解公式如下：

??x1?y1/p1?x2?(y ?2?r2x1)/p2?xi?(yi?zixi?2?rixi?1)/pi,i?3,...,501(2-5)

(2-6)

(2-7)

?u501?x501? ?u500?x500?t500u501?u?x?tu?qx,i?499,...,1iii?1ii?2?i (2-8)

3 cond(A)2及A行列式求解

cond(A)2?|?1| (2-9)

|?s| 由式(2-5)可得：

501detA??pi

i?1

三 源程序

#include #include

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;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(\\\t= %.12e\\t迭代次数：%d\\n\ for(i=0;i<501;i++) //位移 B[i]=A[i]-beta_1;

beta_501=power(B)+beta_1;

printf(\\\t= %.12e\\t迭代次数：%d\\n\ beta_s=inv_power(A);

printf(\\\t= %.12e\\t迭代次数：%d\\n\

for(k=1;k<=39;k++) //第二问 {

mu=beta_1+k*(beta_501-beta_1)/40;

(2-10)