北航数值分析大作业一

北京航空航天大学

数值分析大作业一

学院名称自动化

专业方向控制工程

学号

学生姓名________

教师_________

日期2021年11月26日

设有501x501的实对称矩阵A,

0.1

其中，%=(1.64-0.024i)sin(0.20-0.64e~(i=1,2,--,501),^=0.16,c=-0.064。矩阵A的特征

值为4a=1,2,•••,501),并且有

儿，2501和4的值。

A的与数4=4+4最接近的特征值4,(k=l,2,-.,39)o

A的(谱范数)条件数co〃d(A)2和行列式detAo

一方案设计

1求4，和4的值。

人为按模最小特征值，田=酶田。可使用反客法求得。

4，分别为最大特征值及最小特征值。可使用事法求出按模最大特征值，如结果

为正，即为结果为负，那么为乙。使用位移的方式求得另一特征值即可。

2求A的与数4=4+&%■*最接近的特征值4(k=1,2,...,39)。

40

题目可看成求以出为偏移量后，按模最小的特征值。即以外为偏移量做位移，使用

反鬲法求出按模最小特征值后，加上外，即为所求。

3求A的(谱范数)条件数co〃d(A)2和行列式detA。

矩阵A为非奇异对称矩阵，可知，

cond(A)^=|—

-4in(I」)

其中Anax为按模最大特征值，2min为按模最小特征值。

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

二算法实现

1幕法

使用如下迭代格式：

任取非零向羸0…,〃俨)7

MT/maxIki(?1)

以=A%

A=sgn(max||)max|uk_x\

终止迭代的控制理论使用I乩-4马/14区£,

实际使用

IIAI-IA-1(2-2)

由于不保存A矩阵中的零元素，只保存主对角元素a[501]及b,c值。那么上式中

uk=A%7简化为：

W(l)=。⑴MD+by(2)+cv(3)

«(2)=与，⑴十。(2))，(2)十byQ)十①(4)

4(500)=.(498)+勿(499)+67(500)y(500)+by(50\)

“(501)=cy(499)+办,(500)+。(50l)y(501)(")

u(i)=cy(i-2)+by(i-1)+a(i)y(i)+by(i+1)+cy(i+2)

"=3,…,499)

2反幕法

使用如下迭代格式：

任取非零向勤。=(瞰)，…,町〉)T

)'i=/T/max|uk_{\(24)

14k=A/工

Pk=sgn(max|uk_x|)-max|uk_x\

其中uk=AIk=>Auk=yk_i，解方程求出uko

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

程如下:

追赶法求LU分解的实现：

%

b

c\=LU

A=c

cb%0l

」「「1(2-5)

P1•/・..

&

二z：，・•«・・,•见99

,,•，’500

Z501OP5O111_

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

Pi=%$=b/a]

6=b,p2=a2-r2t}

q{-c!pi,i=1,...499

y=S-WT)/Pi，i=2,…50()(2-6)

Zj=cj=3,

r=/?-cf_9,z=3,...501

Pi=%_c%_2-aT，i=3,...501

推导出回代求解公式如下：

玉二y/巧

，々=()‘2-")/〃2(2-7)

=(X-z/"2-Um)/2,j=3,...5°1

M501=X501

,"500=/00—，5(X)〃501Q-8)

%=xi-tiui+i-qixi+29i=499,...J

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

czmd(A)2=q(2-9)

1A1

由式(2-5)可得：

501

detA=Pi(2-10)

1=1

三源程序

#include<sldio.h>

#include<math.h>

doubleep=le-12,b=0.16,c=-0.064;

intj=0;

doublepower(doublea[501]);〃暴法

doubleinv_powcr(doublca[501J);〃反器法

doubledet(doublea[501]);〃求det

intmain()〃主程序

(

inti,k;

doubleA[5()1],B[5()1],beta_l,beta_501,beta_s,beta_k;

doublemu;

for(i=0;i<501;i++)

A[i]=(1.64-0.024*(i+l))*sin(0.2*(i+1))-0.64*exp(0.1/(i+1));

beta_l=power(A);〃第一问

printf("Xl\t=%.12e\t迭代次数：%d\n\beta_lj);

for(i=0;i<501;i++)//位移

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

beta_501=power(B)+beta_l;

printf(MX501\t=%.12e\t迭代次数：%d\n'\bcta_501,j);

beta_s=inv_power(A);

printf(MXs\t=%.12e\t迭代次数:%d\nH,beta_s,j);

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

{

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

for(i=0;i<501;i++)

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

beta_k=inv_power(B)+mu;

printf("Xi%d\t=%.12e\t迭代次数:%d\nH,k,beta_kj);

)

printf("cond(A)2=%.12e^nu,beta_l/beta_s);〃第三问

printf("detA\t=%.l2e\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]=l;

j=0;

while(j<N)

(

for(i=0;i<501;i++)

y[i]=u[i]/fabs(m);

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

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

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

u[500]=c：!：yl498]+b*y[499]+al500]*y[5001;

for(i=2;i<499;i++)

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

beta=O;

for(i=0;i<501;i++)

(

if(fabs(u[i])>=fabs(bcta))

beta=u[i];

}

if(beta<0)

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

break;

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

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[l]=b;

p[l]=a[l]-r[l]*t[O];

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

t[l]=(b-r[l]*q[O])/p[l];

for(i=2;i<501;i++)

(

rfil=b-c*t[i-21;

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

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

t[i]=(b-r[i]*q[i-l])/p[i];

)

for(i=0;i<501;i++)

u[i]=l;

j=0;

while(j<N)

(

for(i=0;i<501;i++)

y[i]=u[i]/fabs(m);

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

u[l]=(y[l]-r[l]*u[O])/p[l];

for(i=2;i<501;i++)

u[i]=(y[i]-c*u[i-2]-r[i]*u[i-l])/p[i];

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

for(i=498;i>=0;i-)

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

beta=O;

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)<ep)

break;

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

break;

m=beia;j++;

}

return1/beta;

)

doubledet(doublea[5()l])〃求det

(

doubledet_A=1;

doublep[5Ol],r[5Ol],t[5Ol],q[5Ol];

inti;

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