北航數(shù)值分析大作業(yè)一

1、數(shù)值分析第一次作業(yè)設(shè)有的實(shí)對(duì)稱矩陣A，其中，。矩陣A的特征值為，并且有1.求，和的值。2.求A的與數(shù)最接近的特征值。3.求A的(譜范數(shù))條件數(shù)和行列式detA。一 方案設(shè)計(jì)1 求，和的值。 是矩陣A最小的特征值，為矩陣最大的特征值。對(duì)矩陣A使用冪法可以求出、其中的一個(gè)，但不確定是哪一個(gè)，設(shè)求出的值為。然后再構(gòu)造矩陣B=A-I，對(duì)矩陣B使用冪法求出其按模最大的特征值，可得、中另外一個(gè)值，為+。比較和+，其中較小者為為，較大者為。為按模最小特征值，,直接利用反冪法可以求得。2 求A的與數(shù)最接近的特征值。 題目可看成求以為偏移量后，按模最小的特征值。即以為偏移量做位移，使用反冪法求出按模最小特征值后

2、，加上，即為所求。3 求A的(譜范數(shù))條件數(shù)和行列式detA。 矩陣A為非奇異對(duì)稱矩陣，可知，(1-1) 其中為按模最大特征值，為按模最小特征值。 detA可由LU分解得到。因LU均為三角陣，則其主對(duì)角線乘積即為A的行列式。二 算法實(shí)現(xiàn)1 冪法 使用如下迭代格式：(2-1) 終止迭代的控制理論使用， 實(shí)際使用(2-2) 由于不保存A矩陣中的零元素，只保存主對(duì)角元素a501及b,c值。則上式中簡(jiǎn)化為：(2-3)2 反冪法 使用如下迭代格式：(2-4) 其中，解方程求出。對(duì)于的求解算法我采用了LU分解法和JACOBI迭代法兩種。 （1）LU分解由于A為5對(duì)角矩陣，選擇追趕法求取LU分解。求解過(guò)程如

3、下：(2-5) 由上式推出分解公式如下：(2-6)推導(dǎo)出回代求解公式如下：(2-7)(2-8)（2）JACOBI迭代法迭代公式如下所示：u1k=-bu2k-1+cu3k-1+y1k-1/a11u2k=-(bu1k-1+bu3k-1cu4k-1)+ y2k-1/a22u500k=-bu498k-1+bu499k-1+cu501k-1+y500k-1/a500500 u501k=-cu499k-1+bu500k-1+y501k-1/a501501uik=-cui-2k-1+bui-1k-1+bui+1k-1+cui+2k-1+yik-1/aii(2-9)3 及A行列式求解 (2-10) 由式(2-

4、5)可得： (2-11)三 源程序#include #include double mifa(double a501); /冪法double fanmifa1(double a501); /反冪法（三角分解法求解U(k)） double det(double a501) ; /求det double fanmifa(double a501);/反冪法（jacobi求解）double ep=1e-12,b=0.16,c=-0.064; int j=0;int main()int i,k;double A501,B501,beta_1,beta_501,beta_s,beta_k,beta_m;d

5、ouble mu;for(i=0;i501;i+)Ai=(1.64-0.024*(i+1)*sin(0.2*(i+1)-0.64*exp(0.1/(i+1); beta_1=mifa(A); for(i=0;ibeta_1)beta_1=beta_1;beta_501=beta_m;elsebeta_1=beta_m;beta_501=beta_1;printf(1t= %.12etn,beta_1);for(i=0;i501;i+) Bi=Ai-beta_1; beta_501=mifa(B)+beta_1; printf(501=t= %.12etn,beta_501); beta_s=f

6、anmifa(A); printf(st= %.12et迭代次數(shù)：%dn,beta_s,j); for(k=1;k=39;k+) mu=beta_1+k*(beta_501-beta_1)/40; for(i=0;i501;i+) Bi=Ai-mu; beta_k=fanmifa(B)+mu; printf(i%dt= %.12et迭代次數(shù)：%dn,k,beta_k,j); printf(cond(A)2= %.12en,beta_1/beta_s); printf(detAt= %.12en,det(A);return 0;double mifa(double a501) /冪法 int i

7、=0,r;double b=0.16,c=-0.064;double u501,y501;double s,lastbeta,beta,sq;for(i=0;i501;i+) ui=1;j=0;while(1)j+;s=fabs(u0);for(i=1;i501;i+)sq=fabs(ui);if(ssq) s=sq;for(i=0;i501;i+)if(fabs(ui)=s) r=i;for(i=0;i501;i+)yi=ui/s;u0=a0*y0+b*y1+c*y2; /迭代求出下一個(gè)u u1=b*y0+a1*y1+b*y2+c*y3; u499=c*y497+b*y498+a499*y4

8、99+b*y500; u500=c*y498+b*y499+a500*y500; for(i=2;i499;i+)ui=c*yi-2+b*yi-1+ai*yi+b*yi+1+c*yi+2;lastbeta=beta;for(i=0;i501;i+)beta=ur/yr; if(fabs(beta-lastbeta)/fabs(beta)ep) break; return beta;double fanmifa1(double a501) /反冪法(用三角分解法求解AU(k)=Y(k-1)中的U(k) double p501,r501,t501,q501,u501,y501;double bet

9、a,m=1;int i,N=1000;p0=a0;t0=b/p0;r1=b;p1=a1-r1*t0;q0=c/p0;q1=c/p1;t1=(b-r1*q0)/p1;for(i=2;i501;i+)ri=b-c*ti-2;pi=ai-c*qi-2-ri*ti-1;qi=c/pi;ti=(b-ri*qi-1)/pi;for(i=0;i501;i+)ui=1;j=0;while(jN)for(i=0;i501;i+)yi=ui/fabs(m); u0=y0/p0;u1=(y1-r1*u0)/p1;for(i=2;i=0;i-)ui=ui-ti*ui+1-qi*ui+2;beta=0; for(i=0

10、;i=fabs(beta) beta=ui; if(beta0) if(fabs(fabs(beta)-fabs(m)/fabs(beta)ep) break; if(fabs(beta-m)/fabs(beta)ep) break; m=beta;j+; return 1/beta;double fanmifa(double a501) /用jacobi迭代法求解的反冪法double u501,y501,lastu501,beta=0,lastbeta;double s=0,max;int i,r;for(i=0;i501;i+) /將初始向量的分量全賦值為1ui=1;j=0;while(1

11、)j+;s=0;for(i=0;i501;i+)s=ui*ui+s;s=sqrt(s);for(i=0;i501;i+) /將y單位化yi=ui/s;while(1)for(i=0;i501;i+)lastui=ui;u0=(-b*lastu1-c*lastu2+y0)/a0;u1=(-b*lastu0-b*lastu2-c*lastu3+y1)/a1;u499=(-c*lastu497-b*lastu498-b*lastu500+y499)/a499;u500=(-c*lastu498-b*lastu499+y500)/a500;for(i=2;i499;i+)ui=(-c*lastui-2

12、-b*lastui-1-b*lastui+1-c*lastui+2+yi)/ai;max=fabs(u0-lastu0);for(i=1;i501;i+)double tran=fabs(ui-lastui);if(maxtran) max=tran;if(max=ep) break;lastbeta=beta;beta=0;for(i=0;i501;i+)beta+=yi*ui;if(fabs(beta-lastbeta)/fabs(beta)=ep) break;return 1/beta;double det(double a501) /求det double det_A=1;doubl

13、e p501,r501,t501,q501;int i;p0=a0;t0=b/p0;r1=b;p1=a1-r1*t0;q0=c/p0;q1=c/p1;t1=(b-r1*q0)/p1;for(i=2;i501;i+)ri=b-c*ti-2;pi=ai-c*qi-2-ri*ti-1;qi=c/pi;ti=(b-ri*qi-1)/pi;for(i=0;i501;i+)det_A=det_A*pi;return det_A;四 程序結(jié)果（1）在反冪法中使用三角分解法，且使用的初始迭代向量為u=(1,1,1)時(shí)程序運(yùn)行結(jié)果如下所示：(2)在反冪法中使用jacobi迭代法的程序運(yùn)行結(jié)果如下，部分特征值無(wú)法

14、計(jì)算，部分特征值計(jì)算結(jié)果正確。五 計(jì)算過(guò)程中的現(xiàn)象（1）在反冪法程序中，對(duì)方程組Auk=yk-1的求解用了三角分解法和JACOBI迭代法兩種方法。使用三角分解法時(shí)結(jié)果正常，但使用JACOBI迭代法時(shí)，題干中求解與k（k=1，2，39）最接近的特征值這一問(wèn)中，部分特征值計(jì)算不出來(lái)，部分特征值可計(jì)算出結(jié)果，且結(jié)果與使用三角分解法的結(jié)果一致，分析原因?yàn)椴糠謽?gòu)造的矩陣不滿足JACOBI迭代方法的收斂條件。（2）迭代初始向量的選擇對(duì)計(jì)算結(jié)果會(huì)產(chǎn)生一定影響，主要表現(xiàn)在收斂速度上。選用初始迭代向量u=(1,0,0,，0)時(shí)，程序的運(yùn)行結(jié)果如下所示，可以看出，程序運(yùn)行結(jié)果與使用初始向量u=(1,1,1)的結(jié)果差別較大，而且迭代次數(shù)明顯增加，前者迭代次數(shù)最多可到356次，后者最多為129次。使用MATLAB進(jìn)行計(jì)算發(fā)現(xiàn)使用u=(1,1,1)向量更加接近準(zhǔn)確值。1、初始迭代向量中元素等于0的個(gè)數(shù)越少，收斂結(jié)果越穩(wěn)定；初始迭代向量中元素等于0的個(gè)數(shù)越多，收斂結(jié)果越不穩(wěn)定； 2、迭代初始值ui=s(i=1,2,501)且s的絕對(duì)值值極大，收斂結(jié)果可以穩(wěn)定但收斂速度減慢，其原因?yàn)閟的數(shù)量級(jí)