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

《数值分析》计算实习题目第一题：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)和detA。1）1=n2cond(A)，其中1和n分别是按模最大和最小特征值。2）利用步骤（1）中分解矩阵A得出的LU矩阵，L为单位下三角阵,U为上三角阵，其中U矩阵的主对角线元素之积即为detA。由于A的元素零元素较多，为节省储存量，将A的元素存为6×501的数组中，程序中采用get_an_element()函数来从小数组中取出A中的元素。2.全部源程序#includestdio.h#includemath.hvoidinit_a();//初始化Adoubleget_an_element(int,int);//取A中的元素函数doublepowermethod(double);//原点平移的幂法doubleinversepowermethod(double);//原点平移的反幂法intpresolve(double);//三角LU分解intsolve(double[],double[]);//解方程组intmax(int,int);intmin(int,int);double(*u)[502]=newdouble[502][502];//上三角U数组double(*l)[502]=newdouble[502][502];//单位下三角L数组doublea[6][502];//矩阵Aintmain(){inti,k;doublelambdat1,lambdat2,lambda1,lambda501,lambdas,mu[40],det;doublelambda[40];init_a();//初始化Alambdat1=powermethod(0);lambdat2=powermethod(lambdat1);lambda1=lambdat1lambdat2?lambdat1:lambdat2;lambda501=lambdat1lambdat2?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;//释放堆内存return0;}voidinit_a()//初始化A{inti;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);}doubleget_an_element(inti,intj)//从A中节省存储量的提取元素方法{if(fabs(i-j)=2)returna[i-j+3][j];elsereturn0;}doublepowermethod(doubleoffset)//幂法{inti,x1;doubleu[502],y[502];doublebeta=0,prebeta=-1000,yita=0;for(i=1;i=501;i++)u[i]=1,y[i]=0;//设置初始向量u[]for(intk=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(intx2=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=%flambda=%ferr=%ek=%d\n,offset,(beta+offset),fabs((prebeta-beta)/beta),k);break;};//输出中间过程,包括偏移量,误差,迭代次数}return(beta+offset);}doubleinversepowermethod(doubleoffset)//反幂法{inti;doubleu[502],y[502];doublebeta=0,prebeta=0,yita=0;for(i=1;i=501;i++)u[i]=1,y[i]=0;//设置初始向量u[]for(intk=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=%flambda=%ferr=%ek=%d\n,offset,(beta+offset),fabs((prebeta-beta)/beta),k);break;};//输出中间过程,包括偏移量,误差,迭代次数}return(beta+offset);}intpresolve(doubleoffset)//三角LU分解{inti,k,j,t;doublesum;for(k=1;k=501;k++)for(j=1;j=501;j++){u[k][j]=l[k][j]=0;if(k==j)l[k][j]=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+l[k][t]*u[t][j];u[k][j]=((k==j)?(get_an_element(k,j)-offset):get_an_element(k,j))-sum;}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+l[i][t]*u[t][k];l[i][k]=(((i==k)?(get_an_element(i,k)-offset):get_an_element(i,k))-sum)/u[k][k];}}return0;}intsolve(doublex[],doubleb[])//解方程组{inti,t;doubley[502];doublesum;y[1]=b[1];for(i=2;i=501;i++){sum=0;for(t=max(1,i-2);t=i-1;t++)sum=sum+l[i][t]*y[t];y[i]=b[i]-sum;}x[501]=y[501]/u[501][501];for(i=500;i=1;i--){sum=0;for(t=i+1;t=min(i+2,501);t++)sum=sum+u[i][t]*x[t];x[i]=(y[i]-sum)/u[i][i];}return0;}intmax(intx,inty){return(xy?x:y);}intmin(intx,inty){return(xy?x:y);}3.计算结果结果如下图所示:部分中间结果：给出了偏移量(offset)，误差(err)，迭代次数(k)4.讨论迭代初始向量的选取对计算结果的影响,并说明原因使用u[i]=1（i=1,2,...,501）作为初始向量进行迭代，可得出以上结果。经过Mathematica计算验证结果正确。现修改初始向量u[1]=1，u[i]=0，(i=2,3,...,501)。得出结果此结果与正确结果相差较多。令初始向量u[m]=1，u[n]=0，(m=1,2,...,250n=251,252,....,501)，得出结果：此结果也与正确结果相差较多。但与上次结果相比，更加靠近准确值一些。再增加初始向量u[]中等于1的元素个数，可以发现其结果更加靠近准确值。经验证，只有当不为0元素的个数达到比较高的一个值时，才能得到精确的收敛结果，与元素的绝对值大小无关。分析算法，设1为按模最大特征值，观察式子01122...nnuxxx和21112211kkknknnuxxx，如果选取的0u使得10，由于计算的精度较高，可能无法求出按模最大的特征值，而求出次大的特征值。实际上，当初始向量0u中的0元素较多时，可能0i的情况较为普遍，许多都有可能等于0，此时计算出的结果便与最大特征值差距较大。特征向量12,,,nxxx是已确定的，因此迭代初始向量的不合理选取完全有可能产生这样的i。