北航梳理统计大作业

2024-12-28

北航梳理统计大作业（共1篇）

1.北航梳理统计大作业篇一

北京航空航天大学

数值分析大作业一

学院名称

自动化

专业方向

控制工程

学

号

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，501和s的值。

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

1i501 1，501分别为最大特征值及最小特征值。可使用幂法求出按模最大特征值，如结果为正，即为501，结果为负，则为1。使用位移的方式求得另一特征值即可。求A的与数k1k501140最接近的特征值ik(k1,2,...,39)。

题目可看成求以k为偏移量后，按模最小的特征值。即以k为偏移量做位移，使用反幂法求出按模最小特征值后，加上k，即为所求。求A的(谱范数)条件数cond(A)2和行列式detA。

矩阵A为非奇异对称矩阵，可知，cond(A)2|max|min

(1-1)其中max为按模最大特征值，min为按模最小特征值。

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

二算法实现幂法

使用如下迭代格式：

(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)反幂法

使用如下迭代格式：

(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)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(“λ1t= %.12et迭代次数：%dn”,beta_1,j);for(i=0;i<501;i++)

//位移

B[i]=A[i]-beta_1;

beta_501=power(B)+beta_1;

printf(“λ501t= %.12et迭代次数：%dn”,beta_501,j);

beta_s=inv_power(A);

printf(“λst= %.12et迭代次数：%dn”,beta_s,j);

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(“λi%dt= %.12et迭代次数：%dn”,k,beta_k,j);

}

printf(“cond(A)2= %.12en”,beta_1/beta_s);

//第三问

printf(“detAt= %.12en”,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;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++;

}

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

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++;

}

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)