北航数值分析大作业第一题幂法与反幂法

《数值分析》计算实习题目

第一题：

1. 算法设计方案

（1）1λ，501λ和s λ的值。

1)首先通过幂法求出按模最大的特征值λt1，然后根据λt1进行原点平移求出另一特征值λt2，比较两值大小，数值小的为所求最小特征值λ1，数值大的为是所求最大特征值λ501。

2)使用反幂法求λs ，其中需要解线性方程组。因为A 为带状线性方程组，此处采用LU 分解法解带状方程组。

（2）与140k λλμλ-5011=+k 最接近的特征值λik 。

通过带有原点平移的反幂法求出与数k μ最接近的特征值 λik 。

（3）2cond(A)和det A 。

1）1=n

λλ2cond(A)，其中1λ和n λ分别是按模最大和最小特征值。 2）利用步骤（1）中分解矩阵A 得出的LU 矩阵，L 为单位下三角阵,U 为上三角阵，其中U 矩阵的主对角线元素之积即为det A 。

由于A 的元素零元素较多，为节省储存量，将A 的元素存为6×501的数组中，程序中采用get_an_element()函数来从小数组中取出A 中的元素。

2.全部源程序

#include <stdio.h>

#include <math.h>

void init_a();//初始化A

double get_an_element(int,int);//取A 中的元素函数

double powermethod(double);//原点平移的幂法

double inversepowermethod(double);//原点平移的反幂法

int presolve(double);//三角LU 分解

int solve(double [],double []);//解方程组

int max(int,int);

int min(int,int);

double (*u)[502]=new double[502][502];//上三角U 数组

double (*l)[502]=new double[502][502];//单位下三角L 数组

double a[6][502];//矩阵A

int main()

{

int i,k;

double lambdat1,lambdat2,lambda1,lambda501,lambdas,mu[40],det;

init_a();//初始化A

lambdat1=powermethod(0);

lambdat2=powermethod(lambdat1);

lambda1=lambdat1<lambdat2?lambdat1:lambdat2;

lambda501=lambdat1>lambdat2?lambdat1:lambdat2;

presolve(0);

lambdas=inversepowermethod(0);

det=1;

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

det=det*u[i][i];

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

{

mu[k]=lambda1+k*(lambda501-lambda1)/40;

presolve(mu[k]);

lambda[k]=inversepowermethod(mu[k]);

}

printf("------------所有特征值如下------------\n");

printf("λ=%1.11e λ=%1.11e\n",lambda1,lambda501);

printf("λs=%1.11e\n",lambdas);

printf("cond(A)=%1.11e\n",fabs(lambdat1/lambdas));

printf("detA=%1.11e \n",det);

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

{printf("λi%d=%1.11e ",k,lambda[k]);if(k % 3==0) printf("\n");} delete []u;

delete []l;//释放堆内存

return 0;

}

void init_a()//初始化A

{

int i;

for (i=3;i<=501;i++) a[1][i]=a[5][502-i]=-0.064;

for (i=2;i<=501;i++) a[2][i]=a[4][502-i]=0.16;

for (i=1;i<=501;i++) a[3][i]=(1.64-0.024*i)*sin(0.2*i)-0.64*exp(0.1/i); }

double get_an_element(int i,int j)//从A中节省存储量的提取元素方法{

if (fabs(i-j)<=2) return a[i-j+3][j];else return 0;

}

double powermethod(double offset)//幂法

{

int i,x1;

double beta=0,prebeta=-1000,yita=0;

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

u[i]=1,y[i]=0;//设置初始向量u[]

for (int k=1;k<=10000;k++)

{

yita=0;

for (i=1;i<=501;i++) yita=sqrt(yita*yita+u[i]*u[i]);

for (i=1;i<=501;i++) y[i]=u[i]/yita;

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

{u[x1]=0;for (int x2=1;x2<=501;x2++)

u[x1]=u[x1]+((x1==x2)?(get_an_element(x1,x2)-offset):get_an_element(x1,x2))*y[x2];} prebeta=beta;beta=0;

for (i=1;i<=501;i++) beta=beta+ y[i]*u[i];

if (fabs((prebeta-beta)/beta)<=1e-12) {printf("offset=%f lambda=%f err=%e k=%d\n",offset,(beta+offset),fabs((prebeta-beta)/beta),k);break;};

//输出中间过程,包括偏移量,误差,迭代次数

}

return (beta+offset);

}

double inversepowermethod(double offset)//反幂法

{

int i;

double u[502],y[502];

double beta=0,prebeta=0,yita=0;

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

u[i]=1,y[i]=0; //设置初始向量u[]

for (int k=1;k<=10000;k++)

{

yita=0;

for (i=1;i<=501;i++) yita=sqrt(yita*yita+u[i]*u[i]);

for (i=1;i<=501;i++) y[i]=u[i]/yita;

solve(u,y);

prebeta=beta;beta=0;

for (i=1;i<=501;i++) beta=beta+ y[i]*u[i];

beta=1/beta;

if (fabs((prebeta-beta)/beta)<=1e-12) {printf("offset=%f lambda=%f err=%e k=%d\n",offset,(beta+offset),fabs((prebeta-beta)/beta),k);break;};

//输出中间过程,包括偏移量,误差,迭代次数

}

return (beta+offset);

}

int presolve(double offset)//三角LU分解