北航数值分析大作业一.docx
《北航数值分析大作业一.docx》由会员分享，可在线阅读，更多相关《北航数值分析大作业一.docx（15页珍藏版）》请在冰豆网上搜索。
北航数值分析大作业一
《数值分析B》大作业一
SY1103120朱舜杰
一．算法设计方案：
1.矩阵A的存储与检索
将带状线性矩阵A[501][501]转存为一个矩阵MatrixC[5][501].
由于C语言中数组角标都是从0开始的，所以在数组MatrixC[5][501]中检索A的带内元素aij的方法是：
A的带内元素aij=C中的元素ci-j+2,j
2.求解λ1，λ501，λs
1首先分别使用幂法和反幂法迭代求出矩阵按摸最大和最小的特征值λmax和λmin。
λmin即为λs；
如果λmax>0,则λ501=λmax；如果λmax<0,则λ1=λmax。
②使用带原点平移的幂法（mifa（）函数），令平移量p=λmax，求出对应的按摸最大的特征值λ，max，
如果λmax>0,则λ1=λ，max+p；如果λmax<0,则λ501=λ，max+p。
3.求解A的与数μk=λ1+k（λ501-λ1）/40的最接近的特征值λik（k=1,2，…，39）。
使用带原点平移的反幂法，令平移量p=μk，即可求出与μk最接近的特征值λik。
4.求解A的（谱范数）条件数cond（A）2和行列式detA。
①cond（A）2=|λ1/λn|，其中λ1和λn分别是矩阵A的模最大和最小特征值。
②矩阵A的行列式可先对矩阵A进行LU分解后，detA等于U所有对角线上元素的乘积。
二．源程序
#include
#include
#include
#include
#include
#include
#include
#defineE1.0e-12/*定义全局变量相对误差限*/
intmax2（inta,intb）/*求两个整型数最大值的子程序*/
{
if（a>b）
returna;
else
returnb;
}
intmin2（inta,intb）/*求两个整型数最小值的子程序*/
{
if（a>b）
returnb;
else
returna;
}
intmax3（inta,intb,intc）/*求三整型数最大值的子程序*/
{intt;
if（a>b）
t=a;
elset=b;
if（treturn（t）;
}
voidassignment（doublearray[5][501]）/*将矩阵A转存为数组C[5][501]*/
{
inti,j,k;
//所有元素归零
for（i=0;i<=4;）
{
for（j=0;j<=500;）
{
array[i][j]=0;
j++;
}
i++;
}
//第0,4行赋值
for（j=2;j<=500;）
{
k=500-j;
array[0][j]=-0.064;
array[4][k]=-0.064;
j++;
}
//第1,3行赋值
for（j=1;j<=500;）
{
k=500-j;
array[1][j]=0.16;
array[3][k]=0.16;
j++;
}
//第2行赋值
for（j=0;j<=500;）
{k=j;
j++;
array[2][k]=（1.64-0.024*j）*sin（（double）（0.2*j））-0.64*exp（（double）（0.1/j））;
}
}
doublemifa（doubleu[501],doublearray[5][501],doublep）/*带原点平移的幂法*/
{
inti,j;/*u[501]为初始迭代向量*/
doublea,b,c=0;/*array[5][501]为矩阵A的转存矩阵*/
doubley[501];/*p为平移量*/
for（;;）
{
a=0;
b=0;
/*选用第一种迭代格式*/
//求ηk-1
for（i=0;i<=500;i++）
{
a=a+u[i]*u[i];
}
a=sqrt（a）;
//求yk-1
for（i=0;i<=500;i++）
{
y[i]=u[i]/a;
}
//求uk
for（i=0;i<=500;i++）
{
u[i]=0;
for（j=max2（i-2,0）;j<=min2（i+2,500）;j++）
{
u[i]+=array[i-j+2][j]*y[j];
}
u[i]=u[i]-p*y[i];/*引入平移量*/
}
//求βk
for（i=0;i<=500;i++）
{
b+=y[i]*u[i];
}
if（fabs（（b-c）/b）<=E）/*达到精度水平，迭代终止*/
break;
c=b;
}
return（b+p）;/*直接返回A的特征值*/
}
voidchuzhi（doublea[]）/*用随机数为初始迭代向量赋值*/
{
inti;
srand（（int）time（0））;
for（i=0;i<=500;i++）
{
a[i]=（10.0*rand（）/RAND_MAX）;/*生成0~10的随机数*/
}
}
voidchuzhi2（doublea[],intj）/*令初始迭代向量为ei*/
{
inti;
for（i=0;i<=500;i++）
{
a[i]=0;
}
a[j]=1;
}
voidLU（doublearray[5][501]）/*对矩阵A进行Doolittle分解*/
{/*矩阵A转存在C[5][501]中*/
intj,k,t;/*分解结果L，U分别存在C[5][501]的上半部与下半部*/
for（k=0;k<=500;k++）
{
for（j=k;j<=min2（（k+2）,500）;j++）
{
for（t=max3（0,k-2,j-2）;t<=（k-1）;t++）
{
array[k-j+2][j]-=array[k-t+2][t]*array[t-j+2][j];
}
}
if（k<500）
for（j=k+1;j<=min2（（k+2）,500）;j++）
{
for（t=max3（0,k-2,j-2）;t<=（k-1）;t++）
{
array[j-k+2][k]-=array[j-t+2][t]*array[t-k+2][k];
}
array[j-k+2][k]=array[j-k+2][k]/array[2][k];
}
}
}
doublefmifa（doubleu[501],doublearray[5][501],doublep）
{/*带原点平移的反幂法*/
inti,j;
doublea,b,c=0;
doubley[501];
//引入平移量
for（i=0;i<=500;i++）
{
array[2][i]-=p;
}
//先将矩阵Doolittle分解
LU（array）;
for（;;）
{
a=0;
b=0;
//求ηk-1
for（i=0;i<=500;i++）
{
a=a+u[i]*u[i];
}
a=sqrt（a）;
//求yk-1
for（i=0;i<=500;i++）
{
y[i]=u[i]/a;
}
//回带过程，求解uk
for（i=0;i<=500;i++）
{
u[i]=y[i];
}
for（i=1;i<=500;i++）
{
for（j=max2（0,（i-2））;j<=（i-1）;j++）
{
u[i]-=array[i-j+2][j]*u[j];
}
}
u[500]=u[500]/array[2][500];
for（i=499;i>=0;i--）
{
for（j=i+1;j<=min2（（i+2）,500）;j++）
{
u[i]-=array[i-j+2][j]*u[j];
}
u[i]=u[i]/array[2][i];
}
//求βk
for（i=0;i<=500;i++）
{
b+=y[i]*u[i];
}
if（fabs（（b-c）/b）<=E）/*达到精度要求，迭代终止*/
break;
c=b;
}
return（p+（1/b））;/*直接返回距离原点P最接近的A的特征值*/
}
//主函数
main（）
{inti;
doubled1,d501,ds,d,a;
doubleu[501];
doubleMatrixC[5][501];
printf（"《数值分析》计算实习题目第一题\n"）;
printf（"SY1103120朱舜杰\n"）;
//将矩阵A转存为MatrixC
assignment（MatrixC）;
//用带原点平移的幂法求解λ1，λ501
chuzhi（u）;
d=mifa（u,MatrixC,0）;
chuzhi（u）;
a=mifa（u,MatrixC,d）;
if（d<0）
{
d1=d;
d501=a;
}
else
{
d501=d;
d1=a;
}
printf（"λ1=%.12e\n",d1）;
printf（"λ501=%.12e\n",d501）;
//用反幂法求λs
chuzhi（u）;
ds=fmifa（u,MatrixC,0）;
printf（"λs=%.12e\n",ds）;
//用带原点平移的反幂法求λik
for（i=1;i<=39;i++）
{
a=d1+（i*（d501-d1））/40;
assignment（MatrixC）;
chuzhi（u）;
d=fmifa（u,MatrixC,a）;
printf（"与μ%02d=%+.12e最接近的特征值λi%02d=%+.12e\n",i,a,i,d）;
}
//求A的条件数
d=fabs（（d1/ds））;
printf（"A的（谱范数）条件数cond2=%.12e\n",d）;
//求detA
assignment（MatrixC）;
LU（MatrixC）;
a=1;
for（i=0;i<=500;i++）
{
a*=MatrixC[2][i];
}
printf（"行列式detA=%.12e\n",a）;
//测试不同迭代初始向量对λ1计算结果的影响。
printf（"改变迭代初始向量对λmax计算结果的测试如下：
\n"）;
assignment（MatrixC）;
for（i=0;i<=500;i++）
{
chuzhi2（u,i）;
d1=mifa（u,MatrixC,0）;
printf（"u%03d,λmax=%+e",i,d1）;
if（（（i+1）%3）==0）
printf（"\n"）;
}
printf（"Pressanykeytocontinue\n"）;
getchar（）;
}
三.程序结果：
四．分析初始向量选择对计算结果的影响
矩阵的初始向量选择，对结果的影响很大，选择不同的初始向量可能会得到的特征值。
以幂法为例（反幂法原理相同），选取初始迭代向量ui=ei（i=0,1，…500）；
即u[j]=0,j≠i；u[j]=1,j=i。
测试结果如下：
试验结果发现只有当i取特定的一些值时才能得到正确的结果，即得到按摸最大的特征值。
i取不同值时，即取不同的初始向量时，可能得到不同的特征值。
这是因为以A的n个线性无关的特征向量为一组基，将初始向量线性表出时，按摸最大的特征值λ1对应的特征向量x1的系数α1如果为0，就无法求出特征值λ1。
如果按摸第二大的特征值λ2对应的特征向量x2的系数α2不为0，则求出该特征值。
若为0，则以此类推。
升级会员 