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

一、算法设计方案

1、矩阵初始化

    矩阵的下半带宽r=2,上半带宽s=2，设置矩阵，在矩阵C中检索矩阵A中的带内元素的方法是：。这样所需要的存储单元数大大减少，从而极大提高了运算效率。

2、利用幂法求出

幂法迭代格式：

当时，迭代终止。

   首先对于矩阵A利用幂法迭代求出一个，然后求出矩阵B，其中(为单位矩阵），对矩阵B进行幂法迭代，求出，之后令，比较，大者为，小者为。

3、利用反幂法求出

反幂法迭代格式：

当时，迭代终止，。

    每迭代一次都要求解一次线性方程组，求解过程为：

（1）作分解

对于执行

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

    使用反幂法，直接可以求得矩阵按模最小的特征值。

    求与数最接近的特征值，对矩阵实行反幂法，即可求出对应的。

4、求出A的条件数和行列式

    根据，其中分子分母分别对应按模最大和最小的特征值。

的计算：由于,其中为下三角矩阵，且对角线元素为1，故，所以有，又为上三角矩阵，故为对其对角线上各元素的乘积，最后可得。

二、程序源代码

（1）定义所需要的函数：

#include

#include

#include

#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

   //对数组进行初始化//

    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.0;

            }            

            else if (i==1||i==3)

            {

                C[i][j]=0.16;

            }    

            else if (i==2)

            {

                C[i][j]=(1.-0.024*(j+1))*sin(0.2*(j+1))

                    -0.*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=int(N/2)时为0），运行结果如下：

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

通过以上类似的实验可以大致看出这样的规律：

的值趋近于有两种情况：

（1）当的元素中，1的个数较多时;

（2）在1的个数相同的条件下，1的分布越靠中后段，

观察对应的特征向量可以发现：

（1）随着i的增加，特征向量元素的绝对值不断增大，即绝对值较大的数集中于中后位置。因此，如果初始向量的非零元素集中在中后段，该初始向量会更容易逼近对应的特征向量，得到的结果也越准确。

对于，初始向量的非零元素集中在前半段的情况进行实验，会发现当算法中不考虑给定的精度水平，强制性执行足够高次数（大约在300多次以上）的迭代，运算结果也会趋近于。这就说明，程序之前没有得到准确结果的原因，是因为迭代次数不够。当迭代次数在100到200次左右时，每一次迭代所造成的相对误差小于给定的精度水平，因此，如果由精度水平来控制循环迭代的次数，程序将错误地判断已经收敛，但实际上，当继续迭代到300次以上时，运算结果会突然变化，直至最终稳定在。

由此，可以得出结论，当迭代次数足够高（300次以上）时，得到的结果会趋于稳定，不同的初始向量和选定的精度水平，决定着程序是否出现以及何时出现假收敛。当所选取初始向量的非零元素越多，以及非零元素的位置越靠后时，收敛会更加迅速、准确。