北航研究生数值分析编程大作业1.doc

数值分析大作业
一、 算法设计方案
1、 矩阵初始化
矩阵的下半带宽r=2,上半带宽s=2，设置矩阵，在矩阵C中检索矩阵A中的带内元素的方法是：。这样所需要的存储单元数大大减少，从而极大提高了运算效率。
2、 利用幂法求出
幂法迭代格式：

当时，迭代终止。
首先对于矩阵A利用幂法迭代求出一个，然后求出矩阵B，其中(为单位矩阵），对矩阵B进行幂法迭代，求出，之后令，比较，大者为，小者为。
3、 利用反幂法求出
反幂法迭代格式：

当时，迭代终止，。
每迭代一次都要求解一次线性方程组，求解过程为：
（1） 作分解
对于执行

（2） 求解（数组b先是存放原方程组右端向量，后来存放中间向量y)

使用反幂法，直接可以求得矩阵按模最小的特征值。
求与数最接近的特征值，对矩阵实行反幂法，即可求出对应的。
4、 求出A的条件数和行列式
根据，其中分子分母分别对应按模最大和最小的特征值。
的计算：由于,其中为下三角矩阵，且对角线元素为1，故，所以有，又为上三角矩阵，故为对其对角线上各元素的乘积，最后可得。

二、 程序源代码
（1）定义所需要的函数：
#include <stdio.h>
#include <conio.h>
#include <math.h>
#define N 501
#define R 2
#define S 2
int min(int a,int b); // 求最小值
int max(int a,int b,int c); // 求最大值
double Fan_two(double x[N]);//计算二范数
void FenjieLU(double (*C)[N]);//解线性方程组的LU分解过程
void Solve(double (*C)[N], double *b,double *x);//解线性方程组的求解过程
double PowerMethod(double C[][N],double u[N],double y[N],double bta,double D);//幂法
double InversePowerMethod(double C[][N],double u[N],double y[N],double bta,double D);//反幂法
};

（2）程序的主函数，Main.cpp代码如下：
void main()
{

double C[R+S+1][N];
double u[N];
double y[N];
double miu[39];
double C1[R+S+1][N];
double bta = 1.0;
double Namda1,Namda501,NamdaS;
double Namda[39];
double CondA2;
double detA = 1.0;
double D = 1.0e-12;
int i, j, k;
FILE * fp;
fp = fopen("Namda.txt","w");
//对数组进行初始化//
int i, j;
for (i = 0; i < N; i++)
{
u[i] = 1;
}
for (i = 0;i< R + S + 1;i++)
{
for (j = 0;j< N;j++)
{
if (i==0||i==4)
{
C[i][j]=-0.064;
}
else if (i==1||i==3)
{
C[i][j]=0.16;
}
else if (i==2)
{
C[i][j]=(1.64-0.024*(j+1))*sin(0.2*(j+1))
-0.64*exp(0.1/(j+1));
}
}
}

//幂法求Namda1//
Namda1 = PowerMethod(C, u, y, bta, D);
printf("\n================================================\n");
printf("Namda1 = %12.11e", Namda1);
printf("\n================================================\n");
//幂法求Namda501//
bta = 1.0;
for (i = 0; i < R + S + 1; i++)
{
for (j = 0; j < N; j++)
{
if (i == 2)
C1[i][j] = C[i][j] -Namda1;
else
C1[i][j] = C[i][j];
}
}
Namda501 = algorism.PowerMethod(C1, u, y, bta, D) +Namda1;
printf("\n================================================\n");
printf("Namda501 = %12.11e", Namda501);
printf("\n================================================\n");
//反幂法求NamdaS//
bta = 1.0;
NamdaS = InversePowerMethod(C, u, y, bta, D);
printf("\n================================================\n");
printf("NamdaS = %12.11e", NamdaS);
printf("\n================================================\n");
//反幂法求Namda[k]//
printf("\n================================================\n");
for (k = 0; k < 39; k++)
{
miu[k] = Namda1 + (k + 1) * (Namda501 - Namda1) / 40.0;
bta = 1.0;
for (i = 0; i < R + S + 1; i++)
{
for (j = 0; j < N; j++)
{
if (i == 2)
C1[i][j] = C[i][j] - miu[k];
else
C1[i][j] = C[i][j];
}
}
Namda[k] = InversePowerMethod(C1, u, y, bta, D) + miu[k];
fprintf(fp,"与%12.11e最接近的特征值为:%12.11e\n",miu[k],Namda[k]);
}
printf("求与miu[k]最接近的Namda[k]的计算结果已经输出到文件Namda.txt中");
printf("\n================================================\n");
//求A的谱范数//
printf("\n================================================\n");
printf("A的谱范数为:%12.11e", sqrt(Namda501));
printf("\n================================================\n");
//求A的条件数//
CondA2 = fabs( Namda1 / NamdaS);
printf("\n================================================\n");
printf("A的谱范数的条件数Cond(A)2为:%12.11e",CondA2);
printf("\n================================================\n");
//求det(A)2的值//
for (j = 0; j < N; j++)
detA *= C[2][j];
printf("\n================================================\n");
printf("行列式A的值为:%12.11e",detA);
printf("\n================================================\n");
fclose(fp);
_getch();
return;
}

（3）成员函数的实现
int min(int a,int b)
{
return a < b ? a : b;
}
int max(int a,int b,int c)
{
int temp;
temp = a > b ? a : b;
return temp > c ? temp : c;
}
double Fan_two(double x[N])
{
double sum = 0.0;
int i;
for (i = 0; i < N; i++)
{
sum += pow(x[i],2);
}
return sqrt(sum);
}

void FenjieLU(double (*C)[N])
{
double sum = 0;
int i, j, k,t;
for (k = 0; k < N; k++)
{
j = k;
i = k + 1;
while (1)
{
if (j == min(k + S + 1, N))
break;
for (t = max(0, k - R, j - S); t <= k - 1; t++)
{
sum += C[k-t+S][t] * C[t-j+S][j];
}
C[k-j+S][j] = C[k-j+S][j] - sum;
sum = 0.0;
j++;
if (k == N-1)
break;
if (i == min(k + R + 1, N))
break;
for (t = max(0, i - R,k - S); t <= k - 1; t++)
{
sum += C[i-t+S][t] * C[t-k+S][k];
}
C[i-k+S][k] = (C[i-k+S][k] - sum) / C[S][k];
sum = 0;
i++;
}
}
}
void Solve(double (*C)[N], double *b,double *x)
{
double sum = 0;
int i, t;
sum = 0;
for (i = 1; i < N; i++)
{
for (t = max(0, i - R); t <= i - 1; t++)
{
sum += C[i-t+S][t] * b[t];
}
b[i] = b[i] - sum;
sum = 0;
}
x[N-1] = b[N-1] / C[S][N-1];
for (i = N - 2; i >= 0; i--)
{
for (t = i+1; t <= min(i + S, N - 1); t++)
{
sum += C[i-t+S][t] * x[t];
}
x[i] = (b[i] - sum) / C[S][i];
sum = 0;
}
}
double PowerMethod(double C[][N],double u[N],double y[N],double bta,double D)
{
double ita;
double sum = 0;
double temp = 0.0;
int i,j,k = 0;
while (fabs(bta - temp) / fabs(bta) > D)
{
temp = bta;
ita = Fan_two(u);
for (i = 0; i < N; i++)
{
y[i] = u[i] / ita;
}
for (i = 0; i < N; i++)
{
for (j = max(0,i - R); j < min(i + S + 1,N); j++)
{
sum += C[i - j + S][j] * y[j];
}
u[i] = sum;
sum = 0;
}
for (i = 0; i < N; i++)
{
sum += y[i] * u[i];
}
bta = sum;
sum = 0;
k++;
}
return bta;
}
double InversePowerMethod(double C[][N],double u[N],double y[N],double bta,double D)
{
double TC[R+S+1][N];
double ty[N];
double ita;
double sum = 0;
double temp = 0.0;
int i,j,k = 0;
FenjieLU(C);
while (abs(1/bta - 1/temp) / abs(1/bta) > D)
{
temp = bta;
ita = Fan_two(u);
for (i = 0; i < N; i++)
{
y[i] = u[i] / ita;
}
//用到临时存储数组TC[][]和ty[][]是因为函数Solve执行过程中会改变A[][]和y[][]
for (i = 0; i < R + S + 1; i++)
{
for (j = 0; j < N; j++)
TC[i][j] = C[i][j];
}
for (i = 0; i < N; i++)
ty[i] = y[i];
Solve(C, y, u);
for (i = 0; i < R+S+1; i++)
{
for (j = 0; j < N; j++)
C[i][j] = TC[i][j];
}
for (i = 0; i < N; i++)
y[i] = ty[i];
for (i = 0; i < N; i++)
{
sum += y[i] * u[i];
}
bta = sum;
sum = 0;
k++;
}
bta = 1.0 / bta;
return bta;
}

三、 程序运行结果
下图为主程序运行结果

其中的结果输出在Namda.txt文件中，结果如下：

四、分析迭代初始向量对计算结果的影响
选择不同的初始向量可能会得到不同的特征值。
选取时，运行结果如下：

选取时，运行结果如下：

选取时（i<N/2时为1，i>=int(N/2)时为0），运行结果如下：

选取时（i<N/2时为0，i>=int(N/2)时为1），运行结果如下：

通过以上类似的实验可以大致看出这样的规律：
的值趋近于有两种情况：
（1） 当的元素中，1的个数较多时;
（2） 在1的个数相同的条件下，1的分布越靠中后段，
观察对应的特征向量可以发现：
（1）随着i的增加，特征向量元素的绝对值不断增大，即绝对值较大的数集中于中后位置。因此，如果初始向量的非零元素集中在中后段，该初始向量会更容易逼近对应的特征向量，得到的结果也越准确。
对于，初始向量的非零元素集中在前半段的情况进行实验，会发现当算法中不考虑给定的精度水平，强制性执行足够高次数（大约在300多次以上）的迭代，运算结果也会趋近于。这就说明，程序之前没有得到准确结果的原因，是因为迭代次数不够。当迭代次数在100到200次左右时，每一次迭代所造成的相对误差小于给定的精度水平，因此，如果由精度水平来控制循环迭代的次数，程序将错误地判断已经收敛，但实际上，当继续迭代到300次以上时，运算结果会突然变化，直至最终稳定在。
由此，可以得出结论，当迭代次数足够高（300次以上）时，得到的结果会趋于稳定，不同的初始向量和选定的精度水平，决定着程序是否出现以及何时出现假收敛。当所选取初始向量的非零元素越多，以及非零元素的位置越靠后时，收敛会更加迅速、准确。

13