数值分析第一次大作业.docx
《数值分析第一次大作业.docx》由会员分享，可在线阅读，更多相关《数值分析第一次大作业.docx（12页珍藏版）》请在冰豆网上搜索。
数值分析第一次大作业
《数值分析》计算实习报告
第一题
院系：
机械工程及自动化学院_
学号：
_____
姓名：
_______
2017年11月7日
一、算法设计方案
1、求，和的值
1）利用幂法计算出矩阵按模最大的特征值，设其为。
2）令矩阵（I为单位矩阵），同样利用幂法计算出矩阵按模最大的特征值。
3）令。
由计算过程可知和分别为矩阵所有特征值按大小排序后，序列两端的值。
即，，。
4）利用反幂法计算。
其中，反幂法每迭代一次都要求解线性方程组，由于矩阵A为带状矩阵，故可用三角分解法解带状线性方程组的方法求解得到。
2、求的与数最接近的特征值
1）令矩阵，利用反幂法计算出矩阵按模最小的特征值,则。
3、求的（谱范数）条件数cond（）2和行列式det
1），前文已算出和，直接带入即可。
2）反幂法计算时，已经对矩阵A进行过Doolittle分解，得到A=LU。
而L为对角线上元素全为1的下三角矩阵，U为上三角矩阵，可知，，即有。
最后，为节省存储量，需对矩阵进行压缩，将中带内元素存储为数组。
二、源程序代码
#include
#include
#include
#include
usingnamespacestd;
#defineN501
#defineK39
#definer2
#defines2
#defineEPSI1.0e-12
//求两个整数中的最大值
intint_max2（inta,intb）
{
return（a>b?
a:
b）;
}
//求两个整数中的最小值
intint_min2（inta,intb）
{
return（aa:
b）;
}
//求三个整数中的最大值
intint_max3（inta,intb,intc）
{
intt;
if（a>b）
t=a;
elset=b;
if（treturn（t）;
}
//定义向量内积
doubledianji（doublex[],doubley[]）
{
doublesum=0;
for（inti=0;isum=sum+x[i]*y[i];
return（sum）;
}
//计算两个数之间的相对误差
doubleerro（doublelamd0,doublelamd1）
{
doublee,d,l;
e=fabs（lamd1-lamd0）;
d=fabs（lamd1）;
l=e/d;
return（l）;
}
//矩阵A的压缩存储初始化成C
voidinit_c（doublec[][N]）
{
inti,j;
for（i=0;ifor（j=0;jif（i==0||i==4）
c[i][j]=-0.064;
elseif（i==1||i==3）
c[i][j]=0.16;
else
c[i][j]=（1.64-0.024*（j+1））*sin（0.2*（j+1））-0.64*exp（0.1/（j+1））;
}
//矩阵复制
voidfuzhi_c（doublec_const[][N],doublec[][N]）
{
inti,j;
for（i=0;ifor（j=0;jc[i][j]=c_const[i][j];
}
//LU三角分解
voidLUDet_c（doublec_const[][N],doublec_LU[][N]）
{
doublesum;
intk,i,j;
fuzhi_c（c_const,c_LU）;
for（k=1;k<=N;k++）
{
for（j=k;j<=int_min2（k+s,N）;j++）
{
sum=0;
for（i=int_max3（1,k-r,j-s）;i<=k-1;i++）
sum+=c_LU[k-i+s][i-1]*c_LU[i-j+s][j-1];
c_LU[k-j+s][j-1]-=sum;
}
for（j=k+1;j<=int_min2（k+r,N）;j++）
{
sum=0;
for（i=int_max3（1,j-r,k-s）;i<=k-1;i++）
sum+=c_LU[j-i+s][i-1]*c_LU[i-k+s][k-1];
c_LU[j-k+s][k-1]=（c_LU[j-k+s][k-1]-sum）/c_LU[s][k-1];
}
}
}
//三角分解法解带状线性方程组
voidjiefc（doublec_const[][N],doubleb_const[],doublex[]）
{
inti,j;
doubleb[N],c_LU[r+s+1][N],sum;
for（i=0;ib[i]=b_const[i];
LUDet_c（c_const,c_LU）;
for（i=2;i<=N;i++）
{
sum=0;
for（j=int_max2（i-2,1）;j<=i-1;j++）
sum+=c_LU[i-j+2][j-1]*b[j-1];
b[i-1]-=sum;
}
x[N-1]=b[N-1]/c_LU[2][N-1];
for（i=N-1;i>=1;i--）
{
sum=0;
for（j=i+1;j<=int_min2（i+2,N）;j++）
sum+=c_LU[i-j+2][j-1]*x[j-1];
x[i-1]=（b[i-1]-sum）/c_LU[2][i-1];
}
}
//幂法求按模最大特征值
doublemifa_c（doublec_const[][N]）
{
doubleu[N],y[N];
doublesum,length_u,beta0,beta1;
inti,j;
for（i=0;iu[i]=0.5;
length_u=sqrt（dianji（u,u））;
for（i=0;iy[i]=u[i]/length_u;
for（i=1;i<=N;i++）
{
sum=0;
for（j=int_max2（i-2,1）;j<=int_min2（i+2,N）;j++）
sum=sum+c_const[i-j+2][j-1]*y[j-1];
u[i-1]=sum;
}
beta1=dianji（u,y）;
do
{
beta0=beta1;
length_u=sqrt（dianji（u,u））;
for（i=0;iy[i]=u[i]/length_u;
for（i=1;i<=N;i++）
{
sum=0;
for（j=int_max2（i-2,1）;j<=int_min2（i+2,N）;j++）
sum=sum+c_const[i-j+2][j-1]*y[j-1];
u[i-1]=sum;
}
beta1=dianji（u,y）;
}while（erro（beta0,beta1）>=EPSI）;
return（beta1）;
}
//反幂法求按模最小特征值
doublefmifa_c（doublec_const[][N]）
{
doubleu[N],y[N];
doublelength_u,beta0,beta1;
inti;
for（i=0;iu[i]=0.5;
length_u=sqrt（dianji（u,u））;
for（i=0;iy[i]=u[i]/length_u;
jiefc（c_const,y,u）;
beta1=dianji（y,u）;
do
{
beta0=beta1;
length_u=sqrt（dianji（u,u））;
for（i=0;iy[i]=u[i]/length_u;
jiefc（c_const,y,u）;
beta1=dianji（y,u）;
}while（erro（beta0,beta1）>=EPSI）;
beta1=1/beta1;
return（beta1）;
}
//计算lamd_1、lamd_501、lamd_s
voidcalculate1（doublec_const[][N],double&lamd_1,double&lamd_501,double&lamd_s）
{
inti;
doublelamd_mifa0,lamd_mifa1,c[r+s+1][N];
lamd_mifa0=mifa_c（c_const）;
fuzhi_c（c_const,c）;
for（i=0;ic[2][i]=c[2][i]-lamd_mifa0;
lamd_mifa1=mifa_c（c）+lamd_mifa0;
if（lamd_mifa0{
lamd_1=lamd_mifa0;
lamd_501=lamd_mifa1;
}
else
{
lamd_501=lamd_mifa0;
lamd_1=lamd_mifa1;
}
lamd_s=fmifa_c（c_const）;
}
//平移+反幂法求最接近u_k的特征值
voidcalculate2（doublec_const[][N],doublelamd_1,doublelamd_501,doublelamd_k[]）
{
i
升级会员 