数值分析第一次计算实习大作业.docx

数值分析第一次计算实习大作业数值分析B第一次计算实习大作业学院名称：宇航学院专业名称：航天工程学生姓名：学号：2017年 11月一、题目：设有的实对称矩阵A，其中，。矩阵A的特征值为，并且有1.求，和的值。2.求A的与数最接近的特征值。3.求A的(谱范数)条件数和行列式detA。说明：1.在所用的算法中，凡事要给出精度水平的，都取=10-12。2.选择算法时，应使矩阵A的所有零元素都不存储。3.打印以下内容：（1）全部源程序；（2）特征值，以及，detA的值。4.采用e型输出实型数，并且至少显示12位有效数字。二、算法设计方案：第一问：求，和的值。，则是A的最大特征值，是A的最小特征值，是A的按模最小特征值。通过幂法可以计算矩阵的按模最大特征值，然后通过判断特征值的正负，若为正，则按模最大特征值为，若为负，则按模最大特征值为。由此计算出，其中一个后，再通过平移的方法求另外一个特征值。则直接通过反幂法计算得到。第二问：求A的与数最接近的特征值。先通过公式计算，将作为偏移量对A进行平移获得新矩阵。然后对新矩阵通过反幂法计算得到按模最小特征值。将所得结果加上，就是所求。第三问：求A的(谱范数)条件数和行列式detA。因为矩阵A为非奇异的实对称矩阵，则有以下公式：其中为按模最大特征值，为按模最小特征值，在第一问中都已得到结果，可直接计算公式得到。矩阵A为五对角矩阵，采用追赶法对矩阵A进行LU分解。A=LU，因为LU均为三角矩阵，则有以下公式：detA=i=1nu(i,i)其中u(i,i)为主对角线上元素，通过此公式计算detA。三、算法函数计算：（1）幂法计算使用第二种幂法迭代格式：终止迭代的控制理论使用。由于A的所有零元素都不存储，因此我采取一个数组ai保存主对角线元素，定义数值常数b和c。因此需要将上式中的进行简化运算，最后得到以下结果：（2）反幂法计算使用如下迭代格式：其中，通过追赶法对矩阵A进行LU分解，求解方程组得到。LU分解如下：推出分解公式如下：回代求解公式如下：（3）追赶法LU分解计算detA=i=1501d(i,i)分解公式与反幂法计算中所采用的分解公式相同。四、源程序代码：#include #include #include using namespace std;const double epsilon=1e-12,b=0.16,c=-0.064;int i,j,k;double power_method(double f501); /定义幂法函数double inv_power_method(double f501); /定义反幂法函数double det(double f501);/定义追赶法LU分解后求行列式的函数int main() /主程序cout数值分析第一次大作业：endl学院：宇航学院 姓名：何海涛 学号：ZY1715114endl计算结果如下：endl;double a501,d501,lambda_1,lambda_501,lambda_s,lambda_k,lambda_x;double move;for(i=0;i501;i+)ai=(1.64-0.024*(i+1)*sin(0.2*(i+1)-0.64*exp(0.1/(i+1);/第一问lambda_x=power_method(a); /通过幂法求出矩阵的按模最大特征值coutendl第一问结果为:endlendl;cout 按模最大特征值 lambda_x=setw(20)setprecision(12)scientificlambda_xendl;if(lambda_x0)/按模最大特征值为负，则为lambda_1lambda_1=lambda_x;cout 最小特征值 lambda_1=setw(20)setprecision(12)scientificlambda_1setw(15)迭代次数：jendl;for(i=0;i501;i+)di=ai-lambda_1;/通过平移求另一个特征值lambda_501=power_method(d)+lambda_1;cout 最大特征值 lambda_501=setw(20)setprecision(12)scientificlambda_501setw(15)迭代次数：jendl;else/按模最大特征值为正，则为lambda_501lambda_501=lambda_x;for(i=0;i501;i+)di=ai-lambda_501;lambda_1=power_method(d)+lambda_501; /通过平移求另一个特征值cout 最小特征值 lambda_1=setw(20)setprecision(12)scientificlambda_1setw(15)迭代次数：jendl;cout 最大特征值 lambda_501=setw(20)setprecision(12)scientificlambda_501setw(15)迭代次数：jendl;lambda_s=inv_power_method(a); /通过反幂法求矩阵的按模最小特征值cout 按模最小特征值 lambda_s=setw(20)setprecision(12)scientificlambda_ssetw(15)迭代次数：jendl;/第二问coutendlendl第二问结果为:endlendl;for(k=1;k=39;k+)for(i=0;i501;i+)move=lambda_1+k*(lambda_501-lambda_1)/40;di=ai-move;lambda_k=inv_power_method(d)+move;if(k10)cout 特征值 lambda_i0k=setw(20)setprecision(12)scientificlambda_ksetw(15)迭代次数：jendl;elsecout 特征值 lambda_ik=setw(20)setprecision(12)scientificlambda_ksetw(15)迭代次数：jendl;/第三问coutendlendl第三问结果为:endlendl;coutsetw(19) A的条件数cond(A)2=setw(20)setprecision(12)scientificlambda_x/lambda_sendl; /A的条件数等于按模最大特征值除以按模最小特征值coutsetw(19) 行列式detA=setw(20)setprecision(12)scientificdet(a)endlendl;double power_method(double f501) /幂法函数int n=5000;double u501,y501;double max=1,lambda;for(i=0;i501;i+)/任取非零向量uiui=1;for(j=0;jn;j+)for(i=0;i501;i+)/计算yiyi=ui/fabs(max);u0=f0*y0+b*y1+c*y2;u1=b*y0+f1*y1+b*y2+c*y3;u499=c*y497+b*y498+f499*y499+b*y500;u500=c*y498+b*y499+f500*y500;for(i=2;i499;i+)/计算u(i)ui=c*yi-2+b*yi-1+fi*yi+b*yi+1+c*yi+2;lambda=0;for(i=0;i=fabs(lambda)lambda=ui;if(lambda0)if(fabs(fabs(lambda)-fabs(max)/fabs(lambda)epsilon)break;if(fabs(lambda-max)/fabs(lambda)epsilon)break;max=lambda;return lambda;double inv_power_method(double f501) /反幂法函数double d501,e501,g501,h501,u501,y501;double max=1,lambda;int n=5000;/LU分解d0=f0;g0=b/d0;e1=b;d1=f1-e1*g0;h0=c/d0;h1=c/d1;g1=(b-e1*h0)/d1;for(i=2;i501;i+)ei=b-c*gi-2;di=fi-c*hi-2-ei*gi-1;hi=c/di;gi=(b-ei*hi-1)/di;/反幂法迭代计算for(i=0;i501;i+)ui=1;for(j=0;jn;j+)/计算yifor(i=0;i501;i+)yi=ui/fabs(max);/计算uiu0=y0/d0;u1=(y1-e1*u0)/d1;for(i=2;i=0;i-)ui=ui-gi*ui+1-hi*ui+2;lambda=0;for(i=0;i=fabs(lambda)lambda=ui;if(lambda0)if(fabs(fabs(lambda)-fabs(max)/fabs(lambda)epsilon)break;if(fabs(lambda-max)/fabs(lambda)epsilon)break;max=lambda;return 1/lambda;double det(double f501) /追赶法LU分解后求行列式的函数double det_A=1;double d501,e501,g501,h501;d0=f0;g0=b/d0;e1=b;d1=f1-e1*g0;h0=c/d0;h1=c/d1;g1=(b-e1*h0)/d1;for(i=2;i501;i+)ei=b-c*gi-2;di=fi-c*hi-2-ei*gi-1;hi=c/di;gi=(b-ei*hi-1)/di;for(i=0;i501;i+)det_A=det_A*di;return det_A;五、程序运行结果：1-1.070011361488e+0015019.724634101479 e+000s-5.557910794214 e-003i1-1.018293403315 e+001i2-9.585707425068 e+000i3-9.172672423928 e+000i4-8.652284007898 e+000i5-8.093483808675 e+000i6-7.659405407692 e+000i7-7.119684648691 e+000i8-6.611764339397 e+000i9-6.066103226595 e+000i10-5.585101052628 e+000i11-5.114083529812 e+000i12-4.578872176865 e+000i13-4.097829307903 e+000i14-3.554211215751 e+000i15-3.041090018133 e+000i16-2.526430459138 e+000i17-2.003230769564 e+000i18-1.503557611227 e+000i19-9.935586060075 e-001i20-4.870426738850 e-001i212.231736249575 e-002i225.324174742067 e-001i231.052898962693 e+000i241.589445881881 e+000i252.060330460274 e+000i262.558075597073 e+000i273.080240509307 e+000i283.580798705450 e+000i294.091378510451 e+000i304.603035378279 e+000i315.132924283898 e+000i325.594906348083 e+000i336.080933857027 e+000i346.680354092112 e+000i357.293877448126 e+000i367.717111714236 e+000i378.225220014050 e+000i3883648666065194 e+000i399.254200344575 e+0001.925204273883 e+003detA2.772786141752 e+118六、计算分析：在计算过程中，不同的初始向量收敛速度不同，计算结果也不会完全相同。最开始选取初始向量为u0=(1,0,0,0),结果存在较大误差。改变初始向量为u0=(1,1,1,1)后，计算结果较准确。因此通过改变幂法函数中选取的初始向量，分析对结果的影响。ui(1,1,1,1)-1.070011361488e+0019.724634101479e+000(1,2,3,501)-1.070011361488e+0019.724634101485e+000(1,1,1,1,0,0)-1.07