北航数值分析大作业一.docx
- 文档编号:9013713
- 上传时间:2023-02-02
- 格式:DOCX
- 页数:15
- 大小:336.89KB
北航数值分析大作业一.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(t 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;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,则以此类推。
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 北航 数值 分析 作业