数值分析第二题.docx
《数值分析第二题.docx》由会员分享，可在线阅读，更多相关《数值分析第二题.docx（12页珍藏版）》请在冰点文库上搜索。
数值分析第二题
一、算法的设计方案
（一）求矩阵特征值
（1）使用拟上三角化的算法把矩阵A化为拟上三角化矩阵
给定精度水平和迭代最大次数L。
（2）记
，令k=1，m=n。
（3）如果
，则得到A的一个特征值
，置
（降阶），转（4）；否则转（5）。
（4）如果m=1，则得到A的一个特征值
，转（10）；如果m=0，则直接转（10）；如果m>1，则转
（2）。
（5）求二阶子阵
的两个特征值
和
，即计算二次方程的两个根。
（6）如果m=2，则得到A的两个特征值
和
，转（11）；否则转（7）。
（7）如果
，则得到A的两个特征值
和
，置
（降阶），转（4）；否则转（8）。
（8）如果k=L，则终止计算，未得到A的全部特征值；否则转（9）。
（9）记
（
），计算
（I是m阶单位阵）
（10）置
，转（4）。
（11）A的全部特征值已计算完毕，停止计算。
（二）求矩阵特征向量
特征向量为方程组
的解向量，采用列主元素Gauss消去法求解方程组的算法：
1、消元过程
对于k=1，2，…，n-1执行
（1）选行号
，使
。
（2）交换
与
以及
与
所含的数值。
（3）对于
计算
2、回代过程
二、
全部源程序
#include
#include
voidCreatA（doublea[10][10]）/*创建矩阵*/
{inti,j;
for（i=0;i<=9;i++）
{for（j=0;j<=9;j++）
{if（i!
=j）
a[i][j]=sin（0.5*（i+1）+0.2*（j+1））;
else
a[i][j]=1.5*cos（（i+1）+1.2*（j+1））;
}
}
}
voidHessenberg（double（*p）[10],intn）/*拟上三角化*/
{intw,i,r,j;doublesum,tr,dr,cr,hr;doubleur[10],pr[10],qr[10],wr[10];
for（r=0;r<=n-3;r++）
{w=0;
for（i=r+2;i<=n-1;i++）
{if（（*（*（p+i）+r））!
=0）{w=1;break;}}
if（w==0）continue;
sum=0;
for（i=r+1;i<=n-1;i++）
{sum=sum+（*（*（p+i）+r））*（*（*（p+i）+r））;}
dr=sqrt（sum）;
if（（*（*（p+r+1）+r））>0）{cr=-dr;}
else{cr=dr;}
hr=cr*cr-cr*（*（*（p+r+1）+r））;
for（i=0;i<=r;i++）
{ur[i]=0;}
ur[r+1]=（*（*（p+r+1）+r））-cr;
for（i=r+2;i<=n-1;i++）
{ur[i]=（*（*（p+i）+r））;}
tr=0;
for（i=0;i<=n-1;i++）
{pr[i]=0;qr[i]=0;
for（j=0;j<=n-1;j++）
{pr[i]=（*（*（p+j）+i）*ur[j]）/hr+pr[i];qr[i]=（*（*（p+i）+j）*ur[j]）/hr+qr[i];}
tr=pr[i]*ur[i]/hr+tr;
}
for（i=0;i<=n-1;i++）
{wr[i]=qr[i]-tr*ur[i];}
for（i=0;i<=n-1;i++）
for（j=0;j<=n-1;j++）
{*（*（p+i）+j）=*（*（p+i）+j）-wr[i]*ur[j]-ur[i]*pr[j];}
}
for（i=0;i<=9;i++）
{for（j=0;j<=9;j++）
printf（"%e",*（*（p+i）+j））;
printf（"\n"）;
}
}
voidMkMAT（doublea[10][10],doubleMk[10][10],doubles,doublet,intm）
{inti,j,k;
for（i=0;i<=m-1;i++）
for（j=0;j<=m-1;j++）
{Mk[i][j]=0;
for（k=0;k<=m-1;k++）
Mk[i][j]=a[i][k]*a[k][j]+Mk[i][j];
}
for（i=0;i<=m-1;i++）
for（j=0;j<=m-1;j++）
{if（i==j）
{Mk[i][j]=Mk[i][j]-s*a[i][j]+t;}
else
{Mk[i][j]=Mk[i][j]-s*a[i][j];}
}
}
voidQRMAT（double（*b）[10],double（*c）[10],intm）
{intflag,i,r,j;doublex,tr,dr,cr,hr;doubleur[10],pr[10],qr[10],wr[10],vr[10];
for（r=0;r<=m-2;r++）
{flag=0;
for（i=r+1;i<=m-1;i++）
{if（（*（*（b+i）+r））!
=0）{flag=1;break;}}
if（flag==0）continue;
x=0;
for（i=r;i<=m-1;i++）
{x=x+（*（*（b+i）+r））*（*（*（b+i）+r））;}
dr=sqrt（x）;
if（（*（*（b+r）+r））>0）{cr=-dr;}
else{cr=dr;}
hr=cr*cr-cr*（*（*（b+r）+r））;
for（i=0;i<=r-1;i++）
{ur[i]=0;}
ur[r]=（*（*（b+r）+r））-cr;
for（i=r+1;i<=m-1;i++）
{ur[i]=（*（*（b+i）+r））;}
tr=0;
for（i=0;i<=m-1;i++）
{pr[i]=0;qr[i]=0;vr[i]=0;
for（j=0;j<=m-1;j++）
{pr[i]=（*（*（c+j）+i）*ur[j]）/hr+pr[i];qr[i]=（*（*（c+i）+j）*ur[j]）/hr+qr[i];vr[i]=（*（*（b+j）+i）*ur[j]）/hr+vr[i];}
tr=pr[i]*ur[i]/hr+tr;
}
for（i=0;i<=m-1;i++）
{wr[i]=qr[i]-tr*ur[i];}
for（i=0;i<=m-1;i++）
for（j=0;j<=m-1;j++）
{*（*（c+i）+j）=*（*（c+i）+j）-wr[i]*ur[j]-ur[i]*pr[j];*（*（b+i）+j）=*（*（b+i）+j）-ur[i]*vr[j];}
}
}
voidDL（doublea[10][10],doubleb[10],doublex[10]）/*杜立特尔分解*/
{
inti,j,k,r;
doublel[10][10],u[10][10],y[10];
for（i=0;i<10;i++）
{
for（j=0;j<10;j++）
{
if（i==j）
l[i][i]=1.0;
else
{
l[i][j]=0.0;
if（i>j）
u[i][j]=0.0;
}
}
}
for（k=0;k<10;k++）
{
for（j=k;j<10;j++）
{
u[k][j]=a[k][j];
for（r=0;r<=k-1;r++）
u[k][j]-=l[k][r]*u[r][j];
}
for（i=k+1;i<10;i++）
{
l[i][k]=a[i][k];
for（r=0;r<=k-1;r++）
l[i][k]-=l[i][r]*u[r][k];
l[i][k]/=u[k][k];
}
}
for（i=0;i<10;i++）
{
y[i]=b[i];
for（j=0;j<=i-1;j++）
y[i]-=l[i][j]*y[j];
}
for（i=9;i>=0;i--）
{
x[i]=y[i];
for（j=i+1;j<10;j++）
y[i]-=u[i][j]*x[j];
x[i]=y[i]/u[i][i];
}
}
voidVector（doubleA[10][10],doublea）/*求特征向量*/
{doubleu[10]={1,2,3,4,5,6,7,8,9,0},b,sum,x,y[10],B[10][10];inti,j,k;
for（i=0;i<=9;i++）
for（j=0;j<=9;j++）
{if（i==j）B[i][j]=A[i][j]-（float）a;
elseB[i][j]=A[i][j];
}
for（k=1;k<=20;k++）
{sum=0;
for（i=0;i<=9;i++）
{sum=sum+u[i]*u[i];}
x=sqrt（sum）;
for（i=0;i<=9;i++）
y[i]=u[i]/x;
DL（B,y,u）;
b=0;
for（i=0;i<=9;i++）
{b=b+y[i]*u[i];}
if（b>1e+12）break;
}
printf（"\n"）;
for（i=0;i<=9;i++）
{printf（"%.13.7e",y[i]）;}
printf（"\n"）;
}
voidmain（）
{doubleA[10][10],Mk[10][10],precision=1e-12,s,t,S1,S2,Re[10],Im[10],a;intL=200,k=1,m=10,i,j;
CreatA（A）;
Hessenberg（A,10）;
while（m>=1）
{if（fabs（A[m-1][m-2]）<=precision）{Re[m-1]=A[m-1][m-1];Im[m-1]=0;m--;}
else
{s=A[m-2][m-2]+A[m-1][m-1];
t=A[m-2][m-2]*A[m-1][m-1]-A[m-1][m-2]*A[m-2][m-1];
S1=s/2;S2=sqrt（fabs（s*s-4*t））/2;
if（m==2）
{if（（s*s-4*t）<0）{Re[m-1]=S1;Im[m-1]=-S2;Re[m-2]=S1;Im[m-2]=S2;}
else{Re[m-1]=S1+S2;Im[m-1]=0;Re[m-2]=S1-S2;Im[m-2]=0;}
break;}
else
{if（fabs（A[m-2][m-3]）<=precision）
{if（（s*s-4*t）<0）{Re[m-1]=S1;Im[m-1]=-S2;Re[m-2]=S1;Im[m-2]=S2;}
else{Re[m-1]=S1+S2;Im[m-1]=0;Re[m-2]=S1-S2;Im[m-2]=0;};m=m-2;}
else{if（k==L）{printf（"fail!
"）;break;}
else{MkMAT（A,Mk,s,t,m）;QRMAT（Mk,A,m）;k=k+1;continue;}
}
}
}
if（m==1）{Re[m-1]=A[m-1][m-1];Im[m-1]=0;m--;}
}
for（m=0;m<=9;m++）
{printf（"%13.7e，%13.7ei\n",Re[m],Im[m]）;
}
for（i=0;i<=9;i++）
{for（j=0;j<=9;j++）
printf（"%e\t",A[i][j]）;
printf（"\n"）;}
CreatA（A）;
for（i=0;i<=9;i++）
{if（Im[i]==0）{a=Re[i];Vector（A,a）;}}
}
三、
计算结果
（1）拟上三角化矩阵及特征值
（2）QR分解后矩阵
（3）实特征值的特征向量
升级会员 