北航数值分析大作业一

北 京 航 空 航 天 大 学

数值分析大作业一

学院名称 自动化 专业方向 控制工程 学 号 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)

for(i=0;i<501;i++) B[i]=A[i]-mu;

beta_k=inv_power(B)+mu;

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

printf(\ //第三问

printf(\ }

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;i<501;i++) u[i]=1; j=0;

while(j

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)

if(fabs(beta-m)/fabs(beta)

return beta; }

double inv_power(double a[501]) //反幂法 {

double p[501],r[501],t[501],q[501],u[501],y[501]; double beta,m=1; int i,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=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)

if(fabs(beta-m)/fabs(beta)

return 1/beta; }

double det(double a[501]) //求det {

double det_A=1;

double p[501],r[501],t[501],q[501]; int i;

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]; return det_A; }

四 程序结果

五 计算过程中的现象

使用 |?k??k?1|/|?k|??作为终止迭代条件时，出现迭代无法终止的情况，通过调试发现按模最大特征值为负时，当k充分大后，迭代向量uk各分量不断变号，使得?k与?k?1异号，判别式|?k??k?1|/|?k|不收敛。

因此将终止迭代条件修改为||?k|?|?k?1||/|?k|??，程序实现如下： if(beta<0)

if(fabs(fabs(beta)-fabs(m))/fabs(beta)

if(fabs(beta-m)/fabs(beta)

从迭代次数可以看出?1与?501收敛较慢，由按模最大特征值与按模次大特征值的比值越小，收敛速度越慢，可知存在与?1和?501的模相近的特征值。

本文来源：https://www.bwwdw.com/article/y34w.html