北航数值分析大作业一(纯原创-高分版)

数值分析上机实习作业一题目：设有501×501的实对称矩阵A，𝐴=[𝑎1𝑏𝑐𝑏⋱⋱⋱𝑐⋱⋱⋱⋱⋱⋱𝑐𝑐𝑏𝑎501]其中𝑎𝑖=(1.64−0.024𝑖)sin(0.2𝑖)−0.64𝑒0.1𝑖,(𝑖=1,2,⋯,501),𝑏=0.16,𝑐=−0.64。矩阵A的特征值为，𝜆𝑖(𝑖1,2,⋯,501,并且有𝜆1≤𝜆2≤⋯≤𝜆501,|𝜆𝑠|=𝑚𝑖𝑛1≤𝑖≤501|𝜆𝑖|1.求λ1,和λ501的值。2.求A的与数𝜇𝑘=𝜆1+𝑘𝜆501−𝜆140最接近的特征值𝜆𝑖𝑘(𝑘=1,2,⋯,39).3.求A的（普范数）条件数cond(A)和行列式detA。说明1.在所用的算法中，凡是要给出精度水平的ε，都取𝜀=10−12。2.选择算法的时候应使矩阵A的所有零元素都不存储。3.打印以下内容：(1)算法设计方案和思路。(2)全部源程序。(3)特征值𝜆1,𝜆501,𝜆𝑠,𝜆𝑖𝑘(𝑘=1,2,⋯,39)以及𝑐𝑜𝑛𝑑(𝐴)2,det𝐴的值(采用e型输出实型数，并且至少显示12位有效数字)。(4)发现的现象与问题。算法设计思路和方案因为题目只要求求部分特征值，及最大的特征值和最小的特征值以及距离某些给定实数最近的特征值，而矩阵A是实对称矩阵，所有特征值均为实数，故我们可以用幂法和反幂法结合适当的平移量求解。其中求解最大的特征值和最小的特征值，采用幂法并结合适当的平移量求解，而求解λs，和距离μk最近的特征值采用带平移的反幂法求解。算法分析一、幂法求解𝝀𝟏,𝝀𝟓𝟎𝟏我们知道幂法适用于求解矩阵A的绝对值最大的特征值，设矩阵A的特征值为|λ1|≥|𝜆2|≥⋯≥|𝜆𝑛|，对于情况|λ1||𝜆2|，幂法是适用的，并可以求出特征值𝜆1，然而对于|λ1|=|𝜆2|时，幂法只对λ1=𝜆2适用(也即是当绝对值最大的特征值是重根的情况，幂法是适用的)，然而对于λ1=−𝜆2，绝对值最大的特征值为互为相反数的情况并不适用。所以在用幂法求解𝜆1,𝜆501时，我们分两种情况讨论。1.𝜆1≠−𝜆501此时，由于𝜆1≤𝜆2≤⋯≤𝜆501，我们知道矩阵A绝对值最大的特征值必然为𝜆1,𝜆501中的一个，而且幂法适用，通过幂法迭代一定次数后就可以得到满足精度要求的特征值λ，这时做如下平移：𝐵=𝐴−𝜆𝐼然后对矩阵B用幂法。由线性代数知识知，矩阵B的特征值为𝜆1−𝜆,𝜆2−𝜆,⋯,𝜆501−𝜆，若𝜆=𝜆1，则有0=𝜆1−𝜆≤𝜆2−𝜆≤⋯≤𝜆501−𝜆，对B用幂法即可求出其绝对值最大的特征值𝜆′=𝜆501−𝜆，也即可以得到A的特征值𝜆501=𝜆′+𝜆。对于另外一种情况𝜆=𝜆501，则有𝜆1−𝜆≤𝜆2−𝜆≤⋯≤𝜆501−𝜆=0，对B用幂法即可求出其绝对值最大的特征值λ′=𝜆1−𝜆，同样也即可以得到A的特征值𝜆1=λ′+𝜆。综上，对于𝜆1≠−𝜆501，只需要对A用幂法即可得到A的绝对值最大的特征值𝜆，然后做平移𝐵=𝐴−𝜆𝐼，对B用幂法即可求出其绝对值最大的特征值𝜆′，那么有{𝜆,𝜆′+𝜆}={𝜆1,𝜆501}。2.𝜆1=−𝜆501由于这种情况，不能直接对矩阵A使用幂法。此时有两种解决办法，一种是做适当平移𝐵=𝐴−𝑝𝐼，对用幂法迭代求出其绝对值最大的特征值𝜆后，则{−(𝑃+𝜆)，𝑃+𝜆}={𝜆1,𝜆501}。另一个解决思路是考虑矩阵𝐵=𝐴2，由线性代数知识知，矩阵B的特征值为𝜆12,𝜆22,⋯,𝜆5012，对其用幂法迭代即可求出其绝对值最大的特征值𝜆=𝜆12=𝜆5012。综合1，2，求解矩阵A的特征值𝜆1,𝜆501算法思路为设定最大迭代次数和收敛精度，然后对矩阵A使用幂法，如果收敛，说明𝜆1≠−𝜆501，则可以得到一个特征值𝜆，然后做平移𝐵=𝐴−𝜆𝐼，对矩阵B使用幂法得到B的一个特征值𝜆′，则有{𝜆1,𝜆501}={𝜆,𝜆′+𝜆}；若达到最大迭代次数还没有达到收敛要求的精度，说明𝜆1=−𝜆501，则做适当平移𝐵=𝐴−𝑝𝐼，对用幂法迭代求出其绝对值最大的特征值𝜆后，则{𝜆1,𝜆501}={−(𝑃+𝜆)，𝑃+𝜆}；或者考虑矩阵𝐵=𝐴2，对B用幂法迭代，则可求出B绝对值最大的特征值𝜆=𝜆12=𝜆5012。二、反幂法求解𝝀𝒔,𝝀𝒊𝒌设𝐵𝜇=(𝐴−𝜇𝐼)−1，由于𝐵𝜇的特征值为1λ1−𝜇,1λ2−𝜇,⋯,1λ501−𝜇，可知以次取𝜇=0,𝑢1,𝜇2,⋯,𝜇39时，矩阵𝐵𝜇绝对值最大的特征值依次为𝜆(𝜇)=1𝜆𝑠,1λ𝑖1−𝜇1,1λ𝑖2−𝜇2,⋯,1λ𝑖39−𝜇39，而𝜆(𝜇)的值可以对矩阵𝐴−𝜇𝐼使用反幂法求出。可知，若取𝜇0=0，从𝑘=0,1,2,⋯,39，对矩阵𝐵𝜇=(𝐴−𝜇𝐼)−1使用幂法，也即对𝐴𝜇𝑘=𝐴−𝜇𝑘𝐼使用反幂法，求得𝜆(𝜇𝑘)，则可得矩阵A的特征值𝜆𝑖𝑘=1𝜆(𝜇𝑘)+𝜇𝑘,(𝑘=0,1,2,⋯,39),𝜆𝑠=𝜆𝑖0，迭代过程中求解线性方程组𝐵𝜇𝑖𝑢𝑘=𝑦𝑘−1，因为𝐵𝜇𝑖是带状矩阵，上下带宽r=s=2，可以采用DoolittleLU分解法求解。三、求𝒄𝒐𝒏𝒅(𝑨)𝟐、𝐝𝐞𝐭𝑨由定义知，𝑐𝑜𝑛𝑑(𝐴)2=||𝐴||2||𝐴−1||2=|𝜆𝑚||𝜆𝑠|,|𝜆𝑚|=𝑚𝑎𝑥1≤𝑖≤501|𝜆𝑖|，𝜆𝑚,𝜆𝑠在前面都已经求出，可以直接用来算出𝑐𝑜𝑛𝑑(𝐴)2。而对于det𝐴，在用法幂法求解𝜆𝑠时，需用求解线性方程组𝐴𝜇𝑘=𝑦𝑘−1，因为A是带状矩阵，上下带宽r=s=2，可以采用DoolittleLU分解法求解。有A=LU，所以det𝐴=det𝐿⋅det𝑈，而L为对角线上元素全是1的下三角矩阵，U为上三角矩阵，可知det𝐿=1,det𝑈=𝛱𝑖=1501𝑈𝑖𝑖，即有det𝐴=𝛱𝑖=1501𝑈𝑖𝑖。最后对于矩阵A的存储，考虑到矩阵A的特殊性，只需定义一个一维数组用来存储矩阵A的对角线上的元素，和两个浮点数存储b，c即可。算法实现流程图幂法求解𝝀𝟏,𝝀𝟓𝟎𝟏流程图反幂法求𝜆𝑠,𝜆𝑖𝑘流程图令k=1选定初始迭代向量uk-1βk-1=0计算yk-1=uk-1/||uk-1||计算uk=Ayk-1计算βk=yk-1Tuk|βk-βk-1|/|βk|10-12?是λ=βkA=A-λIk10000?否是k=k+1否令k=1选定初始迭代向量uk-1βk-1=0计算yk-1=uk-1/||uk-1||计算uk=AAyk-1计算βk=yk-1Tuk|βk-βk-1|/|βk|10-12?k=k+1k=1选定初始迭代向量uk-1βk-1=0计算yk-1=uk-1/||uk-1||计算uk=Ayk-1计算βk=yk-1Tuk|βk-βk-1|/|βk|10-12?k=k+1否开始是λ0=βk+λλ=βk是输出λ1=min(λ0,λ)λ501=max(λ0,λ)输出λ1=-sqrt(λ)λ501=sqrt(λ)结束否幂法求A的最大和最小特征值流程图k10000?是否输出请调整初始迭代向量k10000?是否输出请调整初始迭代向量令k=1选定初始迭代向量uk-1βk-1=0计算yk-1=uk-1/||uk-1||利用前面算好的LU分解反解Auk=yk-1计算βk=yk-1Tuk|βk-βk-1|/|βk|10-12?k10000?否是k=k+1开始令μ0=0，t=0μi=λ1+i*(λ501-λ1)/40计算A=A-μtI对A进行LU分解得矩阵L,U是输出λt=1/βk+μtt=t+1t=39?是否结束否请调整初始迭代向量尝试重新计算λtt=t+1反幂法求λsλik流程图计算结果𝜆1-1.070011361502E+01𝜆𝑖19-9.935586060075E-01𝜆5019.724634099619E+00𝜆𝑖20-4.870426738850E-01𝜆𝑠-5.557910794229E-03𝜆𝑖212.231736249575E-02𝜆𝑖1-1.018293403315E+01𝜆𝑖225.324174742069E-01𝜆𝑖2-9.585707425068E+00𝜆𝑖231.052898962693E+00𝜆𝑖3-9.172672423928E+00𝜆𝑖241.589445881881E+00𝜆𝑖4-8.652284007898E+00𝜆𝑖252.060330460274E+00𝜆𝑖5-8.093483808675E+00𝜆𝑖262.558075597073E+00𝜆𝑖6-7.659405407692E+00𝜆𝑖273.080240509307E+00𝜆𝑖7-7.119684648691E+00𝜆𝑖283.613620867692E+00𝜆𝑖8-6.611764339397E+00𝜆𝑖294.091378510451E+00𝜆𝑖9-6.066103226595E+00𝜆𝑖304.603035378279E+00𝜆𝑖10-5.585101052628E+00𝜆𝑖315.132924283898E+00𝜆𝑖11-5.114083529812E+00𝜆𝑖325.594906348083E+00𝜆𝑖12-4.578872176865E+00𝜆𝑖336.080933857027E+00𝜆𝑖13-4.096470926260E+00𝜆𝑖346.680354092112E+00𝜆𝑖14-3.554211215751E+00𝜆𝑖357.293877448127E+00𝜆𝑖15-3.041090018133E+00𝜆𝑖367.717111714236E+00𝜆𝑖16-2.533970311130E+00𝜆𝑖378.225220014050E+00𝜆𝑖17-2.003230769563E+00𝜆𝑖388.648666065193E+00𝜆𝑖18-1.503557611227E+00𝜆𝑖399.254200344575E+00𝑐𝑜𝑛𝑑(𝐴)21.925204273902e+003det𝐴2.772786141752e+118初始迭代向量为：𝒖𝟎⃗⃗⃗⃗=(1,1,1,⋯,1)。(计算𝜆501所使用的𝒖𝟎⃗⃗⃗⃗=(1,2,3,⋯,501)分析初始向量对结果的影响分析初始向量对计算结果的影响，以计算𝜆1,𝜆501为例(迭代步数上限为K=10000)。初始迭代向量𝝀𝟏计算结果迭代次数正确与否是否收敛𝒖𝟎⃗⃗⃗⃗=(1,1,1,⋯,1)-1.070011361502e+001341是是𝒖𝟎⃗⃗⃗⃗=(0,1,0,⋯,0)=𝑒2-2.080981085336e+000158否是𝝁𝟎⃗⃗⃗⃗=𝑒63-2.080981085338e+000184否是𝒖𝟎⃗⃗⃗⃗=𝑒64-1.070011361509e+001412是是初始迭代向量𝝀𝟓𝟎𝟏计算结果迭代次数正确与否是否收敛𝝁𝟎⃗⃗⃗⃗=(1,1,1,⋯,1)9.723340141234e+00010000是否𝒖𝟎⃗⃗⃗⃗=(1,2,3,⋯,501)9.724634098771e+000641是是𝒖𝟎⃗⃗⃗⃗=𝑒29.724634099653e+00081否是𝒖𝟎⃗⃗⃗⃗=𝑒649.724634099652e+00073是是𝒖𝟎⃗⃗⃗⃗=𝑒5019.724634098776e+0001