北航硕士研究生数值分析大作业一.docx

北航硕士研究生数值分析大作业一.docx

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

北航硕士研究生数值分析大作业一

数值分析

—计算实习作业一

学院：

17系

专业：

精密仪器及机械

姓名：

张大军

学号：

DY1417114

2014-11-11

数值分析计算实现第一题报告

1、算法方案

算法方案如图1所示。

（此算法设计实现完全由本人独立完成）

2、全部源程序

全部源程序如下所示

#include

#include

#include

intmain（）

{

doublea[501];

doublevv[5][501];

doubled=0;

doubler[3];

doubleuu;

inti,k;

doublemifayunsuan（double*a,doubleweiyi）;

doublefanmifayunsuan（double*a,doubleweiyi）;

voidyasuo（double*A,double（*C）[501]）;

voidLUfenjie（double（*C）[501]）;

//赋值语句

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

{

a[i-1]=（1.64-0.024*i）*sin（0.2*i）-0.64*exp（0.1/i）;

}

//程序一：

使用幂方法求绝对值最大的特征值

r[0]=mifayunsuan（a,d）;

//程序二：

使用幂方法求求平移λ[0]后绝对值最大的λ，得到原矩阵中与最大特征值相距最远的特征值

d=r[0];

r[1]=mifayunsuan（a,d）;

//比较λ与λ-λ[0]的大小，由已知得

if（r[0]>r[1]）

{

d=r[0];

r[0]=r[1];

r[1]=d;

}

//程序三：

使用反幂法求λ

r[2]=fanmifayunsuan（a,0）;

cout<

:

right）;

cout<<"λ["<<1<<"]="<

:

scientific）<

cout<<"λ["<<501<<"]="<

:

scientific）<

cout<<"λ[s]="<

:

scientific）<

//程序四：

求A的与数u最接近的特征值

for（k=1;k<40;k++）

{

uu=r[0]+k*（r[1]-r[0]）/40;

cout<<"最接近u["<

<

:

scientific）<

}

//程序五：

谱范数的条件数是绝对值最大的特征值除以绝对值最小的特征值的绝对值

cout<<"cond（A）2="<

//程序六：

A的行列式的值就是A分解成LU之U的对角线的乘积

yasuo（a,vv）;

LUfenjie（vv）;

uu=1;

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

{

uu=uu*vv[2][i];

}

cout<<"Det（A）="<

return1;

}

doublemifayunsuan（double*a,doubleweiyi）

{

inti,k;

doubleb=0.16;

doublec=-0.064;

doubleee,w,v1,v2,mm,sum;

doubleu[501];

doubley[505]={0};

for（i=0;i<501;i++）u[i]=1;

//给u赋初值

if（weiyi!

=0）

{

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

a[i]-=weiyi;

}

ee=1;k=0;

//使得初始计算时进入循环语句

while（ee>1e-12）

{

mm=0;

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

{

mm=mm+u[i]*u[i];

}

w=sqrt（mm）;

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

{

y[i+2]=u[i]/w;

//注意此处编程与书上不同，之后会解释它的巧妙之处1

}

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

{

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

//1显然巧妙之处凸显出来

}

sum=0;

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

{

sum+=y[i+2]*u[i];

}

v1=v2;

v2=sum;

//去除特殊情况，减少漏洞

if（k==0）

{

k++;

}

else

{

ee=fabs（v2-v1）/fabs（v2）;

}

}

if（weiyi!

=0）

{

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

a[i]+=weiyi;

}//还原A矩阵

return（v2+weiyi）;

}

doublefanmifayunsuan（double*a,doubleweiyi）

{

inti,k;

doubleb=0.16;

doublec=-0.064;

doubleee,w,v1,v2,mm,sum;

doubleu[501];

doubley[501];

doubleC[5][501];

voidyasuo（double*A,double（*C）[501]）;

voidLUfenjie（double（*C）[501]）;

voidqiuU（double（*C）[501],double*y,double*u）;

//把A阵压缩到C阵中

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

u[i]=1;//给u赋初值

if（weiyi!

=0）

{

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

a[i]-=weiyi;

}

yasuo（a,C）;

LUfenjie（C）;

ee=1;k=0;//使得初始计算时进入循环语句

while（ee>1e-12）

{

mm=0;

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

{

mm=mm+u[i]*u[i];

}

w=sqrt（mm）;

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

{

y[i]=u[i]/w;

}

qiuU（C,y,u）;

sum=0;

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

{

sum+=y[i]*u[i];

}

v1=v2;

v2=sum;

//去除特殊情况，减少漏洞

if（k==0）

{

k++;

}

else

{

ee=fabs（1/v2-1/v1）/fabs（1/v2）;

}

}

if（weiyi!

=0）

{

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

a[i]+=weiyi;

}//还原A矩阵

return（1/v2+weiyi）;

}

voidyasuo（double*A,double（*C）[501]）

{

doubleb=0.16;

doublec=-0.064;

inti;

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

{

C[0][i]=c;

C[1][i]=b;

C[2][i]=A[i];

C[3][i]=b;

C[4][i]=c;

}

}

voidLUfenjie（double（*C）[501]）

{

intk,t,j;

intr=2,s=2;

doublesum;

intminn（int,int）;

intmaxx（int,int）;

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

{

for（j=k;j<=minn（k+s,501-1）;j++）

{

if（k==0）sum=0;

else

{

sum=0;

for（t=maxx（k-r,j-s）;t

{

sum=sum+C[k-t+s][t]*C[t-j+s][j];

}

}

C[k-j+s][j]=C[k-j+s][j]-sum;

}

for（j=k+1;j<=minn（k+r,501-1）;j++）

{

if（k<501-1）

{

if（k==0）sum=0;

else

{

sum=0;

for（t=maxx（j-r,k-s）;t

{

sum=sum+C[j-t+s][t]*C[t-k+s][k];

}

}

C[j-k+s][k]=（C[j-k+s][k]-sum）/C[s][k];

}

}

}

}

voidqiuU（double（*C）[501],double*y,double*u）

{

inti,t;

doubleb[501];

doublesum;

intr=2,s=2;

intminn（int,int）;

intmaxx（int,int）;

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

{

b[i]=y[i];

}

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

{

sum=0;

for（t=maxx（0,i-r）;t

{

sum=sum+C[i-t+s][t]*b[t];

}

b[i]=b[i]-sum;

}

u[500]=b[500]/C[s][500];

for（i=501-2;i>=0;i--）

{

sum=0;

for（t=i+1;t<=minn（i+s,500）;t++）

{

sum=sum+C[i-t+s][t]*u[t];

}

u[i]=（b[i]-sum）/C[s][i];

}

}

intminn（intx,inty）

{

intmin;

if（x>y）

min=y;

else

min=x;

returnmin;

}

intmaxx（intb,intc）

{

intmax;

if（b>c）

{

if（b>0）

max=b;

else

max=0;

}

else

{

if（c>0）

max=c;

else

max=0;

}

returnmax;

}

三、特征值以及的值

λ[1]=-1.070011361502e+001λ[501]=9.724634098777e+000

λ[s]=-5.557910794230e-003

最接近u[1]的特征值为-1.018293403315e+001

最接近u[2]的特征值为-9.585707425068e+000

最接近u[3]的特征值为-9.172672423928e+000

最接近u[4]的特征值为-8.652284007898e+000

最接近u[5]的特征值为-8.093483808675e+000

最接近u[6]的特征值为-7.659405407692e+000

最接近u[7]的特征值为-7.119684648691e+000

最接近u[8]的特征值为-6.611764339397e+000

最接近u[9]的特征值为-6.066103226595e+000

最接近u[10]的特征值为-5.585101052628e+000

最接近u[11]的特征值为-5.114083529812e+000

最接近u[12]的特征值为-4.578872176865e+000

最接近u[13]的特征值为-4.096470926260e+000

最接近u[14]的特征值为-3.554211215751e+000

最接近u[15]的特征值为-3.041090018133e+000

最接近u[16]的特征值为-2.533970311130e+000

最接近u[17]的特征值为-2.003230769563e+000

最接近u[18]的特征值为-1.503557611227e+000

最接近u[19]的特征值为-9.935586060075e-001

最接近u[20]的特征值为-4.870426738850e-001

最接近u[21]的特征值为2.231736249575e-002

最接近u[22]的特征值为5.324174742069e-001

最接近u[23]的特征值为1.052898962693e+000

最接近u[24]的特征值为1.589445881881e+000

最接近u[25]的特征值为2.060330460274e+000

最接近u[26]的特征值为2.558075597073e+000

最接近u[27]的特征值为3.080240509307e+000

最接近u[28]的特征值为3.613620867692e+000

最接近u[29]的特征值为4.091378510451e+000

最接近u[30]的特征值为4.603035378279e+000

最接近u[31]的特征值为5.132924283898e+000

最接近u[32]的特征值为5.594906348083e+000

最接近u[33]的特征值为6.080933857027e+000

最接近u[34]的特征值为6.680354092112e+000

最接近u[35]的特征值为7.293877448127e+000

最接近u[36]的特征值为7.717111714236e+000

最接近u[37]的特征值为8.225220014050e+000

最接近u[38]的特征值为8.648666065193e+000

最接近u[39]的特征值为9.254200344575e+000

cond（A）2=1.925204273902e+003Det（A）=2.772786141752e+118

四、现象讨论

在大作业的程序设计过程当中，初始向量的赋值我顺其自然的设为第一个分量为1，其它分量为0的向量，计算结果与参考答案存在很大差别，计算结果对比如下图2所示（左侧为正确结果，右侧为错误结果），导致了我花了很多的时间去检查程序算法。

通过总结和分析原因，本作业的矩阵是个带状的矩阵，当只赋予u某些分量值时，可能在给定精度条件下不能将u完全转移到特征向量的情况，由此使得计算结果出现较大误差。