北航数值分析作业第一题

一、算法设计方案

利用带状Dollittle分解，将A[501][501]转存到数组C[5][501]，以节省存储空间

1、计算λ1和λ501 

    首先使用幂法求出矩阵的按模最大的特征值λ0：如果λ0>0，则其必为按模最大值，因此λ501=λ0，然后采用原点平移法，平移量为λ501，使用幂法迭代求出矩阵A-λ501I的按模最大的特征值，其特征值按从小到大排列应为λ1-λ501、λ2-λ501、……、。因此A-λ501I的按模最大的特征值应为λ1-λ501。又因为λ501的值已求得，由此可直接求出λ1。

2、计算λS

    λS为矩阵A按模最小的特征值，可以通过反幂法直接求出。

3、计算λik

    λik是对矩阵A进行λik平移后，再用反幂法求出按模最小的特征值λmin，λik=λik+λmin。

4、计算矩阵A的条件数计算cond（A）2和行列式det(A)

    矩阵A的条件数为，其中λ1和λn分别是矩阵A的模最大和最小特征值，直接利用上面求得的结果直接计算。

    矩阵A的行列式可先对矩阵A进行LU分解后，det(A)等于U所有对角线上元素的乘积。

二、源程序：

#include

#include

#include

#include

#define s 2

#define r 2

int Max(int v1,int v2);   

int Min(int v1,int v2); 

int maxt(int v1,int v2,int v3); 

void storage(double C[5][501],double b,double c); 

double mifa(double C[5][501]);

void LU(double C[5][501]);

double fmifa(double C[5][501]);

int Max(int v1,int v2)  //求两个数的最大值

{    return((v1>v2)?v1:v2);

}   

int Min(int v1,int v2)      //求两个数最小值

{ return ((v1} 

int maxt(int v1,int v2,int v3)   //求三个数最大值

{ int t;

if(v1>v2) t=v1;

     else t=v2;

if(t     return(t);

}

/***将矩阵值转存在一个数组里，以节省存储空间***/    

void storage(double C[5][501],double b,double c)  

{    int i=0,j=0;

    C[i][j]=0,C[i][j+1]=0;

for(j=2;j<=500;j++)

        C[i][j]=c;

    i++;

    j=0;

    C[i][j]=0;

for(j=1;j<=500;j++)

        C[i][j]=b;

    i++;

for(j=0;j<=500;j++)

        C[i][j]=(1.-0.024*(j+1))*sin(0.2*(j+1))-0.*exp(0.1/(j+1));

    i++;

for(j=0;j<=499;j++)

        C[i][j]=b;

    C[i][j]=0;

    i++;

for(j=0;j<=498;j++)

        C[i][j]=c;

    C[i][j]=0,C[i][j+1]=0;

}                                     

//用于求解最大的特征值，幂法

double mifa(double C[5][501])

{    int m=0,i,j;

    double b2,b1=0,sum;

    double u[501],y[501];

for (i=0;i<501;i++)

{   u[i] = 1.0;

}             

    do

    {    sum=0;

        if(m!=0)b1=b2;

        m++;

     for(i=0;i<=500;i++)

            sum+=u[i]*u[i];

     for(i=0;i<=500;i++)

            y[i]=u[i]/sqrt(sum);

     for(i=0;i<=500;i++)

        {   u[i]=0;

         for(j=Max(i-r,0);j<=Min(i+s,500);j++)

                u[i]=u[i]+C[i-j+s][j]*y[j];

        }

        b2=0;

     for(i=0;i<=500;i++)

            b2=b2+y[i]*u[i];

    }

while(fabs(b2-b1)/fabs(b2)>=1.0e-12);

    return b2;

}

/*****行列式LU分解*****/

void LU(double C[5][501]) 

{    double sum;

     int k,i,j;

for(k=1;k<=501;k++)

     { for(j=k;j<=Min(k+s,501);j++)

            { sum=0;

         for(i=maxt(1,k-r,j-s);i<=k-1;i++)

                  sum+=C[k-i+s][i-1]*C[i-j+s][j-1];

                C[k-j+s][j-1]-=sum;

            }

         for(j=k+1;j<=Min(k+r,501);j++)

            {   sum=0;

         for(i=maxt(1,j-r,k-s);i<=k-1;i++)

                sum+=C[j-i+s][i-1]*C[i-k+s][k-1];

                C[j-k+s][k-1]=(C[j-k+s][k-1]-sum)/C[s][k-1];

            }

        }

}

/***带状DOOLITE分解，并且求解出方程组的解***/

void solve(double C[5][501],double x[501],double b[501])

{    int i,j,k,t;

    double B[5][501],c[501];

for(i=0;i<=4;i++)

    { for(j=0;j<=500;j++)

        B[i][j]=C[i][j];

    }

for(i=0;i<=500;i++)

        c[i]=b[i];

for(k=0;k<=500;k++)

    { for(j=k;j<=Min(k+s,500);j++)

        { for(t=Max(0,Max(k-r,j-s));t<=k-1;t++)

                B[k-j+s][j]=B[k-j+s][j]-B[k-t+s][t]*B[t-j+s][j];

        }

     for(i=k+1;i<=Min(k+r,500);i++)

        {  for(t=Max(0,Max(i-r,k-s));t<=k-1;t++)

                B[i-k+s][k]=B[i-k+s][k]-B[i-t+s][t]*B[t-k+s][k];

            B[i-k+s][k]=B[i-k+s][k]/B[s][k];

        }

    }

for(i=1;i<=500;i++)

     for(t=Max(0,i-r);t<=i-1;t++)

            c[i]=c[i]-B[i-t+s][t]*c[t];

        x[500]=c[500]/B[s][500];

     for(i=499;i>=0;i--)

        {    x[i]=c[i];

         for(t=i+1;t<=Min(i+s,500);t++)

                x[i]=x[i]-B[i-t+s][t]*x[t];

            x[i]=x[i]/B[s][i];

        }

}

//用于求解模最大的特征值，反幂法

double fmifa(double C[5][501])

{    int m=0,i;

    double b2,b1=0,sum=0,u[501],y[501];

for (i=0;i<=500;i++)

    { [i] = 1.0;

    } 

    do

    {    if(m!=0)b1=b2;

        m++;

        sum=0;

     for(i=0;i<=500;i++)

            sum+=u[i]*u[i];

     for(i=0;i<=500;i++)

            y[i]=u[i]/sqrt(sum);

        solve(C,u,y);

        b2=0;

     for(i=0;i<=500;i++)

            b2+=y[i]*u[i];

}while(fabs(b2-b1)/fabs(b2)>=1.0e-12);

    return 1/b2;

}

/***主程序***/

void main()

{    double b=0.16,c=-0.0,det=1.0;

    int i;

    double C[5][501],cond;

    storage(C,b,c);        //进行C的赋值

    cout.precision(12);   //定义输出精度 

    double k1=mifa(C);    //利用幂法计算矩阵的最大特征值和最小特征值

if(k1<0)

      printf("λ1=%.12e\\n",k1);

else if(k1>=0)

      printf("λ501=%.12e\\n",k1);

for(i=0;i<501;i++)

         C[2][i]=C[2][i]-k1;

    double k2=mifa(C)+k1;

if(k2<0)

         printf("λ1=%.12e\\n",k2);

else if(k2>=0)

         printf("λ501=%.12e\\n",k2);

    storage(C,b,c); 

    double k3=fmifa(C);   //利用反幂法计算矩阵A的按模最小特征值

    printf("λs=%.12e\\n",k3);

    storage(C,b,c);           //计算最接近特征值

    double u[39]={0};

for(i=0;i<39;i++)

    {    u[i]=k1+(i+1)*(k2-k1)/40;

        C[2][i]=C[2][i]-u[i];

        u[i]=fmifa(C)+u[i];

        printf("与数u%d 最接近的特征值λ%d:  %.12e\\n",i+1,i+1,u[i]);

    }

if(k1>0) //计算矩阵A的条件数，取2范数

        cond=fabs(k1/k3);

else if(k1<0)

          cond=fabs(k2/k3);

    storage(C,b,c); 

    LU(C);            //利用LU分解计算矩阵A的行列式

for(i=0;i<501;i++)

     det*=C[2][i];

    printf("\\ncond(A)=%.12e\\n",cond);

    printf("\\ndet(A)=%.12e\\n",det);

}

三、计算结果：

四、结果分析

迭代初始向量的选择对果有一定的影响，选择不同的初始向量可能会得到不同阶的特征值。

若改初始向量u[1]=1，u[i]=0，(i=2,3,...,501)。则得出结果

此结果与正确结果相差较多。

当初始向量u[m]=1，u[n]=0，(m=1，2，……，250；n=251，252，……501)，得出结果：

此结果也与正确结果相差较多。但与上次结果相比，更加靠近准确值一些。

再增加初始向量u[i]中等于1的元素个数，可以发现其结果更加靠近准确值。经验证，只有当不为0元素的个数达到比较高的一个值时，才能得到精确的收敛结果，与元素的绝对值大小无关。