资源描述
数值分析第一次上机作业 班级：ZY16131学号：ZY16131姓名：一、算法方案设计（1）存储矩阵A（参考课本25页）:矩阵A501501为大型带状矩阵，上半带宽S=2，下半带宽R=2，参照课本，可将其用循环语句存储在一个5行501列的二维数组A5501中，使得矩阵A的第j列存放数组A 的第j列带状元素，并使得矩阵的主对角元素存放在数组的第三行。检索矩阵的元素时只要按照公式：aij=数组ai-j+3j即可。（2）求1，501，s：第一步，用幂法求出矩阵A按模最大特征值，再对其判断正负，如果是正数，则该特征值为501，如果是负数，则该特征值为1。幂法实现过程具体为：第二步，用反幂法求矩阵A按模最小特征向量S。第三步，由于12.501，可采用原点平移法对矩阵A平移1，得到矩阵（-1I+A），记为矩阵B，再用用幂法求出B的模最大特征值B，则501=B+1。（3）距离k最近的特征值：还是用原点平移法，将矩阵A平移k个单位，再用反幂法求出平移后矩阵模最小特征值k，矩阵A与k最接近的特征值为：ik = k + k =k + （4）A的谱范数条件数与detA：a、求矩阵A的行列式值：在用反幂法时需要对矩阵A进行Doolittle三角分解，A=LU，根据det(A)=det（LU）=det（L）*det（U），其中det（L）=1，det(A)即为U矩阵的对角线元素的乘积。b、求A的（谱范数）条件数cond（A）2:由于A是非奇异实对称阵，从而cond（A）2 =1/s。2、 运行结果分析（1） 取初始向量为 u_0501=1,11，计算结果截图如下：（2）讨论初始迭代向量取值对计算结果的影响在编程实现算法的过程中遇到了很多问题，多次尝试才得以解决。例如，对各个函数的定义后需要对矩阵A进行初始化；为简化程序，方便改变变量值，应该尽可能地将某些多级变量写成函数的形式，只需要对初级变量赋值即可。猜想：幂法迭代的精度跟初始向量 u_0501的取值有关。取初始向量为1,1,0T（共50个1），计算结果如下：对比初始向量为 1,1,1T 计算结果可发现，最小特征值和按模最小特征值均发生较大变化。取初始向量为1,1,0T（共100个1），计算结果如下：计算结果与初始向量为 1,1,1T 时的计算结果基本一样。取初始向量为1,1,0T（共300个1），计算结果如下：对比初始向量为 1,1,1T 计算结果可发现，最大、最小特征值以及按模最小特征值都发生变化，其中最小特征值和按模最小特征值相较初始向量中含0较多的情况更接近初始向量为 1,1,1T 的计算结果，但最大特征值变化比较明显。取初始向量为10,10,10T（元素全为10），计算结果如下：计算结果跟初始向量为 1,1,1T 时的计算结果一样，没有变化。取初始向量为0,0,0T（元素全为0），计算结果如下：计算出错。总结：通过以上五组计算结果对比，可以初步发现：采用幂法迭代时，选取的初始向量中0元素的多少对于计算结果和运算次数有很大的影响。基本规律是：0元素越少，结果越准确。由于对初始向量进行了单位化，所以初始向量的模的大小对于运算结果没有太大影响。这是因为，如果所选的u0使得1=0，由于计算过程中的舍入误差影响，必然会使得迭代过程中的某一步出现x1方向分量不为零的情况，相当于重新进行迭代，增大计算量和误差。（参考课本49页）附：C程序代码#include#include #define N 502 /*数组记法：从a1,1到a501,501*/ #define accuracy 1.0e-12 #define r 2 #define s 2doublec6N; /*定义c矩阵存储压缩后的带状矩阵*/ double fuzhi(); /*赋值函数*/void Doolittle(); /*Doolittle三角分解程序*/int max(int a,int b); /*求两个数较大值函数*/int min(int a,int b); /*求两个数较小值函数*/double mifa(); /*幂法计算*/double fanmifa(); /*反幂法计算*/double fuzhi() /*定义赋值程序，将带状函数逆时针旋转45后，按行赋值，行从1到5，列从1到501，未赋值区域值为0*/int i,j;i=1;for(j=3;jN;j+)cij=-0.064;i=2;for(j=2;jN;j+)cij=0.16;i=3;for(j=1;jN;j+)cij=(1.64-0.024*j)*sin(0.2*j)-0.64*exp(0.1/j);i=4;for(j=1;jN-1;j+)cij=0.16;i=5;for(j=1;jb)?a:b);int min(int a,int b) /*定义求最小值函数*/ return(ab)?a:b);void Doolittle() /*L-U分解程序，采用的是带状矩阵压缩存储后的LU分解法*/ double temp; int i,j,k,t; for(k=1;kN;k+) for(j=k;j=min(k+s,N-1);j+) temp=0; for(t=max(1,max(k-r,j-s);t=(k-1);t+) temp=temp+ck-t+s+1t*ct-j+s+1j; ck-j+s+1j=ck-j+s+1j-temp; for(i=k+1;i=min(k+r,N-1);i+) temp=0; for(t=max(1,max(i-r,k-s);t=(k-1);t+) temp=temp+ci-t+s+1t*ct-k+s+1k; ci-k+s+1k=(ci-k+s+1k-temp)/cs+1k; double mifa() /*幂法计算程序*/ double u0N,u1N;double temp,yita,beta=0,beta0; /*beta是当前的beta值，beta0是上一次的beta值*/int i,j; for(i=1;iN;i+) /*选取迭代初始向量*/u0i=1; do beta0=beta; temp=0; for(i=1;iN;i+) temp=temp+u0i*u0i; yita=sqrt(temp); for(i=1;iN;i+) u1i=u0i/yita; for(i=1;iN;i+) temp=0; for(j=max(i-1,1);j=min(i+2,N-1);j+) temp=temp+ci-j+s+1j*u1j; /* u=A*y */u0i=temp; /*新的u0*/ temp=0; for(i=1;i=accuracy); /*判断精度水平是否满足，不满足则继续迭代*/return(beta); /*满足精度要求之后返回beta值*/double fanmifa() /*反幂法计算程序*/double u0N,u1N,u2N;double temp,yita,beta=0,beta0; int i,t; for(i=1;iN;i+) /*选取迭代初始向量*/u0i=1;do beta0=beta; temp=0; for(i=1;iN;i+) temp=temp+u0i*u0i; yita=sqrt(temp); for(i=1;iN;i+)u1i=u0i/yita;u2i=u1i; /* 对应书本上y=u1*/ fuzhi(); Doolittle(); /*Au=y中,为求解u先对A进行L-U分解*/ for(i=2;iN;i+) temp=0;for(t=max(1,i-r);t=1;i-)temp=0;for(t=i+1;t=min(i+s,N-1);t+)temp=temp+ci-t+s+1t*u0t;u0i=(u2i-temp)/cs+1i;temp=0; for(i=1;i=accuracy); /*迭代运行条件判断*/return(beta);main() /*主函数*/ double u40; /*定义数组，用于计算第2问*/double lambda_1,lambda_501,lambda_s,a,b,d,cond,det;int i,j,k;fuzhi();a=mifa(); /*幂法计算按模最大值*/fuzhi();for(j=1;jb) /*比较两个按模最大值大小，并相应输出最大特征值和最小特征值*/ lambda_1=b; lambda_501=a; printf(Lmd_1=%13.11en,lambda_1); printf(Lmd_501=%13.11en,lambda_501); else lambda_1=a; lambda_501=b; printf(Lmd_1=%13.11en,lambda_1); printf(Lmd_501=%13.11en,lambda_501); fuzhi();d=fanmifa(); /*反幂法计算按模最小值*/printf(Lmd_s=%13.11en,d); /*输出按模最小特征值s*/printf(n); for(k=1;k40;k+) /*对每一个进行移项反幂法运算，求出最接近k的特征值并输出*/ uk=(lambda_501-lambda_1)*k/40+lambda_1; fuzhi(); for(j=1;jN;j+) c3j=c3j-uk; lambda_s=fanmifa()+uk; i=k; printf(u%d=%13.11e Lmd_i%d=%13.11en,i,uk,i,lambda_s); printf(n); cond=fabs(a/d); printf(cond(A)2=%13.11en,cond); /*计算A条件数并输出*/ fuzhi(); /*计算A的行列式值并输出*/ Doolittle(); det=1; for(i=1;iN;i+) det=det*c3i; printf(detA=%13.11en,det); printf(n);
展开阅读全文