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

1、数值分析计算实习题目第一题：1. 算法设计方案（1），和的值。1)首先通过幂法求出按模最大的特征值t1，然后根据t1进行原点平移求出另一特征值t2，比较两值大小，数值小的为所求最小特征值1，数值大的为是所求最大特征值501。2)使用反幂法求s，其中需要解线性方程组。因为A为带状线性方程组，此处采用LU分解法解带状方程组。（2）与最接近的特征值。通过带有原点平移的反幂法求出与数最接近的特征值 。（3）和。1），其中和分别是按模最大和最小特征值。2）利用步骤（1）中分解矩阵A得出的LU矩阵，L为单位下三角阵,U为上三角阵，其中U矩阵的主对角线元素之积即为。由于A的元素零元素较多，为节省储存量，将A

2、的元素存为6×501的数组中，程序中采用get_an_element()函数来从小数组中取出A中的元素。2.全部源程序#include <stdio.h>#include <math.h>void init_a();/初始化Adouble get_an_element(int,int);/取A中的元素函数double powermethod(double);/原点平移的幂法double inversepowermethod(double);/原点平移的反幂法int presolve(double);/三角LU分解int solve(double ,double

3、 );/解方程组int max(int,int);int min(int,int);double (*u)502=new double502502;/上三角U数组double (*l)502=new double502502;/单位下三角L数组double a6502;/矩阵Aint main() int i,k; double lambdat1,lambdat2,lambda1,lambda501,lambdas,mu40,det;double lambda40; init_a();/初始化A lambdat1=powermethod(0);lambdat2=powermethod(lamb

4、dat1);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*uii;for (k=1;k<=39;k+)muk=lambda1+k*(lambda501-lambda1)/40;presolve(muk);lambdak=inversepowermethod(muk);printf("

5、-所有特征值如下-n");printf("=%1.11e =%1.11en",lambda1,lambda501);printf("s=%1.11en",lambdas);printf("cond(A)=%1.11en",fabs(lambdat1/lambdas);printf("detA=%1.11e n",det);for (k=1;k<=39;k+)printf("i%d=%1.11e ",k,lambdak);if(k % 3=0) printf("n&quo

6、t;);delete u;delete l;/释放堆内存return 0;void init_a()/初始化Aint i;for (i=3;i<=501;i+) a1i=a5502-i=-0.064;for (i=2;i<=501;i+) a2i=a4502-i=0.16;for (i=1;i<=501;i+) a3i=(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 ai-j+3j;el

7、se return 0;double powermethod(double offset)/幂法int i,x1;double u502,y502;double beta=0,prebeta=-1000,yita=0;for (i=1;i<=501;i+)ui=1,yi=0;/设置初始向量u for (int k=1;k<=10000;k+)yita=0;for (i=1;i<=501;i+) yita=sqrt(yita*yita+ui*ui);for (i=1;i<=501;i+) yi=ui/yita;for (x1=1;x1<=501;x1+)ux1=0;

8、for (int x2=1;x2<=501;x2+)ux1=ux1+(x1=x2)?(get_an_element(x1,x2)-offset):get_an_element(x1,x2)*yx2;prebeta=beta;beta=0;for (i=1;i<=501;i+) beta=beta+ yi*ui;if (fabs(prebeta-beta)/beta)<=1e-12) printf("offset=%f lambda=%f err=%e k=%dn",offset,(beta+offset),fabs(prebeta-beta)/beta),

9、k);break;/输出中间过程,包括偏移量,误差,迭代次数return (beta+offset);double inversepowermethod(double offset)/反幂法int i;double u502,y502;double beta=0,prebeta=0,yita=0;for (i=1;i<=501;i+)ui=1,yi=0; /设置初始向量u for (int k=1;k<=10000;k+)yita=0;for (i=1;i<=501;i+) yita=sqrt(yita*yita+ui*ui);for (i=1;i<=501;i+) y

10、i=ui/yita;solve(u,y);prebeta=beta;beta=0;for (i=1;i<=501;i+) beta=beta+ yi*ui;beta=1/beta;if (fabs(prebeta-beta)/beta)<=1e-12) printf("offset=%f lambda=%f err=%e k=%dn",offset,(beta+offset),fabs(prebeta-beta)/beta),k);break; /输出中间过程,包括偏移量,误差,迭代次数 return (beta+offset);int presolve(dou

11、ble offset)/三角LU分解int i,k,j,t;double sum;for (k=1;k<=501;k+)for (j=1;j<=501;j+)ukj=lkj=0;if (k=j) lkj=1; /初始化LU矩阵for (k=1;k<=501;k+)for (j=k;j<=min(k+2,501);j+)sum=0;for (t=max(1,max(k-2,j-2) ; t<=(k-1) ; t+)sum=sum+lkt*utj;ukj=(k=j)?(get_an_element(k,j)-offset):get_an_element(k,j)-su

12、m;if (k=501) continue;for (i=k+1;i<=min(k+2,501);i+)sum=0;for (t=max(1,max(i-2,k-2);t<=(k-1);t+)sum=sum+lit*utk;lik=(i=k)?(get_an_element(i,k)-offset):get_an_element(i,k)-sum)/ukk;return 0;int solve(double x,double b)/解方程组int i,t;double y502;double sum;y1=b1;for (i=2;i<=501;i+)sum=0;for (t=

13、max(1,i-2);t<=i-1;t+)sum=sum+lit*yt;yi=bi-sum;x501=y501/u501501;for (i=500;i>=1;i-)sum=0;for (t=i+1;t<=min(i+2,501);t+)sum=sum+uit*xt;xi=(yi-sum)/uii;return 0;int max(int x,int y)return (x>y?x:y);int min(int x,int y)return (x<y?x:y);3.计算结果结果如下图所示:部分中间结果：给出了偏移量(offset)，误差(err)，迭代次数(k)4.讨论迭代初始向量的选取对计算结果的影响,并说明原因使用ui=1（i=1,2,.,501）作为初始向量进行迭代，可得出以上结果。经过Mathematica计算验证结果正确。现修改初始向量u1=1，ui=0，(i=2,3,.,501)。得出结果此结果与正确结果相差较多。令初始向量um=1，un=0，(m=1,2,.,250 n=251,252,.,501)，得出结果：此结果也与正确结果相差较多。但与上次结果相比，更加靠近准确值一些。再增加初始向量u中等于1的元素个数，可以发现其