北航数值分析A第一题.docx
- 文档编号:10777154
- 上传时间:2023-02-22
- 格式:DOCX
- 页数:14
- 大小:226.86KB
北航数值分析A第一题.docx
《北航数值分析A第一题.docx》由会员分享,可在线阅读,更多相关《北航数值分析A第一题.docx(14页珍藏版)》请在冰豆网上搜索。
北航数值分析A第一题
《数值分析A》计算实习第一题
一、设计方案
由于
为带状矩阵,且除对角线其他各带上的元素相同,为节约存贮空间,设置数组
存储
对角线上的元素,设置
,
两个变量存储其他带上的元素。
1.采用幂法求出矩阵
按模最大的特征值
,再进行原点平移采用幂法求
按模最大的特征值
,若
,则
,
;否则
,
。
采用反幂法求矩阵
按模最小的特征值
,涉及一系列系数矩阵相同的线性方程组求解,定义三角分解函数将
分解为单位下三角矩阵
和上三角矩阵
。
2.采用循环语句,每求出一个
值,就将
中的元素都减去
,采用反幂法求解,与求
相同,需调用三角分解函数求解线性方程组。
3.
的条件数
,行列式
的值为三角分解上三角阵
对角线元素的乘积。
二、源代码
#include
#include
#defineM5
#defineN501
#defineS2
#defineR2//定义各项矩阵参数
#defineK2000//最大迭代次数
voidA(double(*a)[N])//定义存储A中元素的转换矩阵,对计算不同平移量时进行初始化
{
inti;
for(i=0;i { a[0][i]=-0.064; a[1][i]=0.16; a[2][i]=(1.64-0.024*(i+1))*sin(0.2*(i+1))-0.64*exp(0.1/(i+1)); a[3][i]=0.16; a[4][i]=-0.064; } } intmax3(intx,inty,intz)//求3个数的最大值 { intt; if(x>y) t=x; else t=y; if(t t=z; returnt; } intmax(intx,inty) { return(x>y? x: y);//求2个数的最大值 } intmin(intx,inty)//求2个数的最小值 { return(x x: y); } voidCF(double*u,double(*a)[N],double*y)//矩阵A与向量y的乘法 { doubleb=0.16,c=-0.064; inti; u[0]=a[2][0]*y[0]+b*y[1]+c*y[2]; u[1]=a[2][1]*y[1]+b*(y[0]+y[2])+c*y[3]; for(i=2;i u[i]=a[2][i]*y[i]+b*(y[i-1]+y[i+1])+c*(y[i-2]+y[i+2]); u[N-2]=a[2][N-2]*y[N-2]+b*(y[N-3]+y[N-1])+c*y[N-4]; u[N-1]=a[2][N-1]*y[N-1]+b*y[N-2]+c*y[N-3]; } doubleMF(double(*a)[N])//幂法 { doubleu[N],y[N],s,t,beta; inti,k; for(i=0;i u[i]=1.0; for(k=0;k { beta=t; s=0.0; t=0.0; for(i=0;i s=s+u[i]*u[i]; s=sqrt(s); for(i=0;i y[i]=u[i]/s; CF(u,a,y); for(i=0;i t=t+y[i]*u[i]; if(fabs((t-beta)/t)<=1e-12) { //printf("迭代次数为%d",k); returnt; } } printf("超过最大迭代次数%d\n",K); printf("相对误差为%.12e\n",fabs((t-beta)/t)); returnt; } voidFJ(double(*a)[N])//Doolittle分解 { doubles; intk,j,t,i; for(k=1;k<=N;k++) { for(j=k;j<=min(k+S,N);j++) { s=0.0; for(t=max3(1,k-R,j-S);t<=k-1;t++) s=s+a[k-t+S][t-1]*a[t-j+S][j-1]; a[k-j+S][j-1]=a[k-j+S][j-1]-s; } for(i=k+1;i<=min(k+R,N)&&k { s=0.0; for(t=max3(1,i-R,k-S);t<=k-1;t++) s=s+a[i-t+S][t-1]*a[t-k+S][k-1]; a[i-k+S][k-1]=(a[i-k+S][k-1]-s)/a[2][k-1]; } } } voidJFC(double*u,double(*a)[N],double*y)//解同系数矩阵方程组 { doubleuu[N],s; inti,t; uu[0]=y[0]; for(i=2;i<=N;i++) { s=0.0; for(t=max(1,i-R);t<=i-1;t++) s=s+a[i-t+S][t-1]*uu[t-1]; uu[i-1]=y[i-1]-s; } u[N-1]=uu[N-1]/a[2][N-1]; for(i=N-1;i>=1;i--) { s=0.0; for(t=i+1;t<=min(i+S,N);t++) s=s+a[i-t+S][t-1]*u[t-1]; u[i-1]=(uu[i-1]-s)/a[S][i-1]; } } doubleFMF(double(*a)[N])//反幂法 {doubleu[N],y[N],s,t,beta; inti,k; FJ(a); for(i=0;i u[i]=1.0; for(k=0;k { beta=t; s=0.0; t=0.0; for(i=0;i s=s+u[i]*u[i]; s=sqrt(s); for(i=0;i y[i]=u[i]/s; JFC(u,a,y); for(i=0;i t=t+y[i]*u[i]; t=1/t; if(fabs((t-beta)/t)<=1e-12) { //printf("迭代次数为%d",k); returnt; } } printf("超过最大迭代次数%d\n",K); printf("相对误差为%.12e\n",fabs((t-beta)/t)); returnt; } voidmain() { doublea[M][N],b[M][N],m1,m2,m3,uk,det=1.0; inti,k; A(a); printf("第1题\n"); m1=MF(a);//幂法计算m1 for(i=0;i b[S][i]=a[S][i]-m1;//原点平移 m2=MF(b);//幂法计算m2 printf("λ1=%.12e\nλ501=%.12e\n",(m1 m1: (m1+m2)),(m1>m1+m2? m1: (m1+m2))); m3=FMF(a);//反幂法计算按模最小的特征值 printf("λs=%.12e\n",m3); printf("第2题\n"); for(k=1;k<40;k++) { uk=(m1 m1: (m1+m2))+(double)k*fabs(m2)/40; A(a); for(i=0;i a[S][i]=a[S][i]-uk;//原点平移 printf("λi%d=%.12e\t",k,FMF(a)+uk);//反幂法计算特征值 if(k%2==0) printf("\n"); } printf("\n第3题\n"); printf("cond(A)=%.12e\n",fabs(m1/m3));//计算条件数 A(a); FJ(a);//将A分解为L与U相乘 for(i=0;i det=det*a[S][i];//计算行列式的值 printf("detA=%.12e\n",det); } 三、计算结果 四、计算结果分析 迭代初始向量 ,采用幂法计算得到矩阵 按模最大的特征值 如下: 结果表明,当 取不同值时,可以得到不同的特征值,当 取某些值时才能得到正确的结果。 这是因为以 的 个线性无关的特征向量为一组基,将初始向量线性表出时,如果按摸最大的特征值 对应的特征向量 的系数 ,则无法求出特征值 。
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 北航 数值 分析 第一
![提示](https://static.bdocx.com/images/bang_tan.gif)