北航数值分析大作业一

《数值分析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

① 首先分别使用幂法和反幂法迭代求出矩阵按摸最大和最小的特征值λ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

#define E 1.0e-12 /*定义全局变量相对误差限*/ int max2(int a,int b) /*求两个整型数最大值的子程序*/ { if(a>b) return a; else return b; } int min2(int a,int b) /*求两个整型数最小值的子程序*/ { if(a>b) return b; else return a; } int max3(int a,int b,int c) /*求三整型数最大值的子程序*/ { int t; if(a>b) t=a; else t=b; if(treturn(t); }

void assignment(double array[5][501]) /*将矩阵A转存为数组C[5][501]*/ { int i,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)); }

} double mifa(double u[501],double array[5][501],double p) /*带原点平移的幂法*/ { int i,j; /* u[501]为初始迭代向量*/ double a,b,c=0; /* array[5][501]为矩阵A的转存矩阵*/ double y[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的特征值*/

} void chuzhi(double a[]) /*用随机数为初始迭代向量赋值*/ { int i; srand((int)time(0)); for(i=0;i<=500;i++) { a[i]=(10.0*rand()/RAND_MAX); /*生成0~10的随机数*/ } } void chuzhi2(double a[],int j) /*令初始迭代向量为ei*/ { int i; for(i=0;i<=500;i++) { a[i]=0; } a[j]=1; }

void LU(double array[5][501]) /*对矩阵A进行Doolittle分解*/ { /*矩阵A转存在C[5][501]中*/ int j,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]; }

} }

double fmifa(double u[501],double array[5][501],double p) { /*带原点平移的反幂法*/ int i,j; double a,b,c=0; double y[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() { int i; double d1,d501,ds,d,a; double u[501]; double MatrixC[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;