北航数值分析大作业一Word下载.docx
《北航数值分析大作业一Word下载.docx》由会员分享,可在线阅读,更多相关《北航数值分析大作业一Word下载.docx(15页珍藏版)》请在冰豆网上搜索。
二.源程序
#include<
stdio.h>
iostream.h>
stdlib.h>
math.h>
float.h>
iomanip.h>
time.h>
#defineE1.0e-12/*定义全局变量相对误差限*/
intmax2(inta,intb)/*求两个整型数最大值的子程序*/
{
if(a>
b)
returna;
else
returnb;
}
intmin2(inta,intb)/*求两个整型数最小值的子程序*/
intmax3(inta,intb,intc)/*求三整型数最大值的子程序*/
{intt;
b)
t=a;
elset=b;
if(t<
c)t=c;
return(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;
k=500-j;
array[0][j]=-0.064;
array[4][k]=-0.064;
j++;
//第1,3行赋值
for(j=1;
array[1][j]=0.16;
array[3][k]=0.16;
//第2行赋值
{k=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++)
a=a+u[i]*u[i];
a=sqrt(a);
//求yk-1
i++)
y[i]=u[i]/a;
//求uk
u[i]=0;
for(j=max2(i-2,0);
=min2(i+2,500);
j++)
{
u[i]+=array[i-j+2][j]*y[j];
}
u[i]=u[i]-p*y[i];
/*引入平移量*/
//求βk
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));
a[i]=(10.0*rand()/RAND_MAX);
/*生成0~10的随机数*/
}
voidchuzhi2(doublea[],intj)/*令初始迭代向量为ei*/
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<
k++)
for(j=k;
=min2((k+2),500);
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;
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)
{/*带原点平移的反幂法*/
//引入平移量
array[2][i]-=p;
//先将矩阵Doolittle分解
LU(array);
//求ηk-1
//回带过程,求解uk
u[i]=y[i];
for(i=1;
for(j=max2(0,(i-2));
=(i-1);
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;
=min2((i+2),500);
u[i]=u[i]/array[2][i];
//求βk
=E)/*达到精度要求,迭代终止*/
return(p+(1/b));
/*直接返回距离原点P最接近的A的特征值*/
//主函数
main()
{inti;
doubled1,d501,ds,d,a;
doubleu[501];
doubleMatrixC[5][501];
printf("
《数值分析》计算实习题目第一题\n"
);
SY1103120朱舜杰\n"
//将矩阵A转存为MatrixC
assignment(MatrixC);
//用带原点平移的幂法求解λ1,λ501
chuzhi(u);
d=mifa(u,MatrixC,0);
a=mifa(u,MatrixC,d);
if(d<
0)
d1=d;
d501=a;
{
d501=d;
d1=a;
}
λ1=%.12e\n"
d1);
λ501=%.12e\n"
d501);
//用反幂法求λs
ds=fmifa(u,MatrixC,0);
λs=%.12e\n"
ds);
//用带原点平移的反幂法求λik
=39;
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));
A的(谱范数)条件数cond<
A>
2=%.12e\n"
d);
//求detA
assignment(MatrixC);
LU(MatrixC);
a=1;
a*=MatrixC[2][i];
行列式detA=%.12e\n"
a);
//测试不同迭代初始向量对λ1计算结果的影响。
改变迭代初始向量对λmax计算结果的测试如下:
\n"
chuzhi2(u,i);
d1=mifa(u,MatrixC,0);
u%03d,λmax=%+e"
i,d1);
if(((i+1)%3)==0)
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,则以此类推。