北航数值分析大作业一.docx

资源描述

北航数值分析大作业一.docx

《北航数值分析大作业一.docx》由会员分享，可在线阅读，更多相关《北航数值分析大作业一.docx（19页珍藏版）》请在冰点文库上搜索。

北航数值分析大作业一.docx

北航数值分析大作业一

《数值分析》大作业一

1、算法设计方案：

矩阵A的存储与检索：

将矩阵A[501][501],转存为新矩阵A[5][501]。

在存入时，按每一行依次存入，存入时原矩阵与新矩阵的列号保持不变，其中原矩阵的主对角线上的元素存为新矩阵的第三行，最后所有元素存入新矩阵中。

求λ1，λ501，λs

λs：

λs表示矩阵的按模最小特征值，为求得λs直接对待求矩阵A应用反幂法即可。

λ1、λ501：

已知矩阵A的特征值满足关系λ1＜…＜…λ501，要求λ1、及λ501时，可按如下方法求解：

a．对矩阵A用幂法，求得按模最大的特征值λn1。

b．按平移量λn1对矩阵A进行原点平移得矩阵，对矩阵B=A+λn1I用反幂法求得B的按模最小特征值λn2。

c．λn3=λn2一λn1

d．则：

λ1=min（λn1，λn3），λ501=max（λn1，λn3）

求A中与数μk=λ1+k（λ501-λ1）/40最接近的特征值λik（k=1，…39）：

先求矩阵B=A-

I对应的按模最小特征值

，则

即为矩阵A与

最接近的特征值。

重复以上过程39次即可求得

（k=0，1，…39）的值。

求A的（谱范数）条件数

和行列式

：

在

（1）中用反幂法求矩阵A的按模最小特征值时，要用到Doolittle分解方法，在Doolittle分解完成后得到的两个矩阵分别为L和U，detA等于U所有对角线上元素的乘积。

，λmax和λmin分别为模最大特征值与模最小特征值。

2、程序源代码：

#include

#defineN501/*定义列数*/

#defineM5/*定义行数*/

#definee1.0e-12/*误差限*/

doubleA[M][N];/*初始矩阵*/

doubleu[N];/*初始向量*/

doubley[N],yy[N];

doublemaximum,TZ1,TZ2,TZ_1,TZ_N,TZ_s,TZ_abs_max;

intmax_sign,max_position;

/*对矩阵A进行初始化程序，存为A[5][501]*/

voidCSH_A（）

{

doubleb=0.16,c=-0.064;

inti;

for（i=2;i

A[0][i]=c;

for（i=1;i

A[1][i]=b;

for（i=0;i

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

for（i=0;i

A[3][i]=b;

for（i=0;i

A[4][i]=c;

}

/*初始化迭代向量,且随机赋值*/

voidCSH_u（）

{

inti;

for（i=0;i

u[i]=double（15*rand（）/32767）;

}

/*幂法中对迭代向量进行单位化程序*/

voidGet_y（）

{

inti;

for（i=0;i

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

}

/*幂法中求解迭代向量的无穷范数*/

voidGet_max（）

{

inti;

max_position=0;

maximum=fabs（u[0]）;

for（i=1;i

{

if（maximum

{

max_position=i;

maximum=fabs（u[i]）;

}

if（u[max_position]<0）

max_sign=-1;

else

max_sign=1;

}

/*获得新迭代向量*/

voidGet_u（）

{

inti;

u[0]=A[2][0]*y[0]+A[1][1]*y[1]+A[0][2]*y[2];

u[1]=A[3][0]*y[0]+A[2][1]*y[1]+A[1][2]*y[2]+A[0][3]*y[3];

u[N-2]=A[4][N-4]*y[N-4]+A[3][N-3]*y[N-3]+A[2][N-2]*y[N-2]+A[1][N-1]*y[N-1];

u[N-1]=A[4][N-3]*y[N-3]+A[3][N-2]*y[N-2]+A[2][N-1]*y[N-1];

for（i=2;i

u[i]=A[4][i-2]*y[i-2]+A[3][i-1]*y[i-1]+A[2][i]*y[i]+A[1][i+1]*y[i+1]+A[0][i+2]*y[i+2];

}

/*获得迭代后特征值*/

voidGet_TZ（）

{

TZ2=TZ1;

TZ1=max_sign*u[max_position];

}

/*求解迭代向量无穷范数的幂法*/

voidCheck_TZ（）

{

inti;

CSH_u（）;

Get_max（）;

Get_y（）;

Get_u（）;

Get_TZ（）;

for（i=0;;i++）

{

Get_max（）;

Get_y（）;

Get_u（）;

Get_TZ（）;

if（fabs（（TZ2-TZ1）/TZ1）

break;

}

/*获取绝对值最大的特征值λ_501*/

voidThe_TZ（）

{

Check_TZ（）;

TZ_abs_max=TZ1;

}

/*获取特征值λ_1*/

voidThe_Other_TZ（）

{

inti;

doubleTZ_temp=TZ1;

for（i=0;i

{

A[2][i]-=TZ_temp;

}

Check_TZ（）;

TZ1+=TZ_temp;

if（TZ1

{

TZ_1=TZ1;

TZ_N=TZ_temp;

}

else

{

TZ_N=TZ1;

TZ_1=TZ_temp;

}

/*两值中取最小*/

intmin（inta,intb）

{

if（a

returna;

else

returnb;

}

/*两值中取最大*/

intmax（inta,intb）

{

if（a

returnb;

else

returna;

}

/*把分解矩阵为LU*/

voidResolve_LU（）

{

intk,i,j,t;

doubletemp;

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

{

for（j=k;j<=min（k+2,N）;j++）

{

temp=0;

for（t=max（max（1,k-2）,j-2）;t<=k-1;t++）

temp+=A[k-t+2][t-1]*A[t-j+2][j-1];

A[k-j+2][j-1]=A[k-j+2][j-1]-temp;

}

for（i=k+1;i<=min（k+2,N）;i++）

{

temp=0;

for（t=max（max（1,i-2）,k-2）;t<=k-1;t++）

temp+=A[i-t+2][t-1]*A[t-k+2][k-1];

A[i-k+2][k-1]=（A[i-k+2][k-1]-temp）/A[2][k-1];

}

/*方程组回代过程*/

voidGO_Back_substitution（）

{

inti,t;

doubletemp=0;

for（i=2;i

{

for（t=max（1,i-2）;t

y[i-1]-=A[i-t+2][t-1]*y[t-1];

}

u[N-1]=y[N-1]/A[2][N-1];

for（i=N-1;i>0;i--）

{

temp=0;

for（t=i+1;t<=min（i+2,N）;t++）

temp+=A[i-t+2][t-1]*u[t-1];

u[i-1]=（y[i-1]-temp）/A[2][i-1];

}

/*求矩阵行列式值*/

doubleDet_matrix（）

{

inti;

doubledet=1;

CSH_A（）;

Resolve_LU（）;

for（i=0;i

det=det*A[2][i];

returndet;

}

/*反幂法中获得迭代向量2一范数*/

doubleGet_norm（）

{

inti;

doublenormal=0;

for（i=0;i

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

normal=sqrt（normal）;

returnnormal;

}

/*反幂法中对迭代向量单位化*/

voidGet_yy（doublenormal）

{

inti;

for（i=0;i

{

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

yy[i]=y[i];

}

/*获得绝对值最小的特征值*/

voidGet_TZ_s（）

{

inti;

TZ2=TZ1;

TZ1=0;

for（i=0;i

TZ1+=yy[i]*u[i];

TZ1=1/TZ1;

}

/*反幂法求绝对值最小的特征值*/

voidTZ_min（）

{

doublenorm=0;

intcount=0;

TZ1=0,TZ2=0;

CSH_u（）;

norm=Get_norm（）;

Get_yy（norm）;

GO_Back_substitution（）;

Get_TZ_s（）;

while（count<10000）

{

count++;

norm=Get_norm（）;

Get_yy（norm）;

GO_Back_substitution（）;

Get_TZ_s（）;

if（fabs（（TZ2-TZ1）/TZ1）

break;

}

TZ_s=TZ1;

}

/*求矩阵条件数*/

doubleGet_cond_A（）/*求矩阵条件数*/

{

doublecond1;

cond1=fabs（TZ_abs_max/TZ_s）;

returncond1;

}

voidTZ_trans_min（）/*偏移条件下反幂法求特征值*/

{

inti,k;

doubletr,uu[40];

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

{

tr=TZ_1+k*（TZ_N-TZ_1）/40;

CSH_A（）;

for（i=0;i

A[2][i]-=tr;

Resolve_LU（）;

TZ_min（）;

TZ_s+=tr;

uu[k]=TZ_s+k*（（TZ_N-TZ_1））/k;

printf（"与数μ%2d=%.12e最近的特征值λ%d=%.12e\n",k,uu[k],k,TZ_s）;

printf（"\n"）;

}

/*主程序*/

voidmain（）

{

doublecond;

doubleTZ_det;

CSH_A（）;/*初始化矩阵A*/

The_TZ（）;/*获取绝对值最大的特征值λ_501*/

The_Other_TZ（）;/*获取特征值λ_1*/

printf（"输出结果\n"）;