北航數(shù)學學院數(shù)值分析作業(yè)1

1、數(shù)值分析作業(yè)1一、算法設(shè)計方案1、 計算1和501時，先用冪法求出A矩陣模最大的特征值，記為max1。然后對矩陣A 平移，求矩陣A-maxI 的模最大的特征值max2。若max2>0，則1=max1，501=max2+max1；否則，1=max1+max2；501=max1。2、 求s時，用反冪法。由于反冪法需要解方程，故先對A三角分解（本例采取的是不列主元的Doolittle 分解），然后每解一次方程，就回代一次。3、 對于ik （k=1，2，···，39），利用帶原點平移的反冪法，求矩陣A-kI的模最小的特征值即可。4、 由于已經(jīng)對A進行了三角分解，則求

2、A的行列式det(A)時，只需求A所有對角元的乘積。其次因為A是非奇異的實對稱矩陣，故A條件數(shù)Cond（A）2=max1 / s。二、源程序（Visual Studio 環(huán)境下的C語言）#include<stdio.h>#include<stdlib.h>#include"math.h"#include"malloc.h"#define N 501 #define s 2 /*對矩陣大小、 /#define s 2 帶寬進行預定義，默認上下帶寬一致 */double Tmax(double (*p)N+1) ; /* 求矩陣的模最

3、大的特征值，采用冪法 */void Fjie( double (*p)N+1 ); /* 對矩陣進行三角分解，采用不列主元 Doolittle分解法 */double *Hdai( double (*p)N+1,double *q); /* 對已三角化的方程的系數(shù)矩陣求解，即為回代過程 */double Tmin(double (*p)N+1) ; /* 求矩陣的模最小的特征值 ，采用反冪法 */ void main() double a6502,b6502,x40; double t1,t2,max,min,det_a,cond_a,y; /* t1-矩陣最小特征值； t2-矩陣最大特征值；

4、max-模最大特征值； min-模最小特征值；det_a-a的行列式；cond_a-a的條件數(shù)；a存儲的是題設(shè)矩陣A的非零元，ai-j+s+1j=Aij; */ int i,j; for(i=1;i<=502;i+) /* 對矩陣進行填充 */ a3i=(1.64-0.024*i)*sin(0.2*i)-0.64*exp(0.1/i); if(i>=3) a1i=-0.064;b1i=-0.064; if(i>=2) a2i=0.16;b2i=0.16; if(i<=500) a4i=0.16;b4i=0.16; if(i<=499) a5i=-0.064;b5i

5、=-0.064; max=Tmax(a); /* 求A模最大的特征值 */ for(i=1;i<=501;i+) b3i=a3i-max; t1=Tmax(b); if(max>0) t1+=max; t2=max; else t2=t1+max; t1=max; /* 此時求得的t1，t2就是1，501 */ Fjie(a); for(i=1;i<=39;i+) y=t1+i*(t2-t1)/40;for(j=1;j<=501;j+) b3j=a3j-y; xi=Tmin(b)+y; /* xi的值就是 ii的值 */ for(i=1,det_a=1;i<=50

6、1;i+) det_a=det_a*a3i; min=Tmin(a); /* min的值就是s的值 */ cond_a=fabs(max/min); printf("第一題答案：n"); printf("1=%.12e,501=%.12e,s=%.12en",t1,t2,min); printf("第二題答案：n"); for(i=1,j=1;i<=39;j+,i+) printf("i%d=%.12e ",i,xi); if(j%2=0) printf("n"); /*每兩個換行 */

7、printf("n第三題答案：n"); printf("det_A=%.12e,cond_A=%.12en",det_a,cond_a); system("pause"); /* 顯示結(jié)果窗口 */double Tmax(double (*p)N+1) /* 求矩陣的模最大的特征值，采用冪法 */* p-指向矩陣 */ double uN+1,yN+1,a,b1=1,b2,sum,t=1; int i,j,k1,k2; for(i=1;i<=N;i+) /*對U向量初始化 */ ui=1; while(fabs(t)>1e

8、-12) /* "1e-12 "為精度 */ for(i=1,sum=0;i<=N;i+)sum+=ui*ui; a=sqrt(sum); for(i=1;i<=N;i+) yi=ui/a; for(i=1;i<=N;i+) k1=i+s; k2=i-s; if(k1>N) k1=501;if(k2<1) k2=1;for(j=k2,sum=0;j<=k1;j+) sum+=*(*(p+i-j+3)+j)*yj; ui=sum; for(i=1,b2=0;i<=N;i+) b2+=yi*ui; t=(b2-b1)/b2; b1=b2

9、; return b1; void Fjie( double (*p)N+1 ) /* 對矩陣進行三角分解，采用不列主元 Doolittle分解法 */* p-指向矩陣 */ double sum; int i,j,k,t,h1,h2; for(k=1;k<=N;k+) h1=k+s; if(h1>N) h1=N;for(j=k;j<=h1;j+) h2=j-s; if(h2<1) h2=1; for(t=h2,sum=0;t<=k-1;t+) sum+=*(*(p+k-t+3)+t)*(*(*(p+t-j+3)+j); *(*(p+k-j+3)+j)-=sum;

10、 if(k<N) for(i=k+1;i<=h1;i+) h2=i-s; if(h2<1) h2=1; for(t=h2,sum=0;t<=k-1;t+) sum+=*(*(p+i-t+3)+t)*(*(*(p+t-k+3)+k); *(*(p+i-k+3)+k)=(*(*(p+i-k+3)+k)-sum)/(*(*(p+3)+k); double *Hdai( double (*p)N+1,double *q) /* 對已三角化的方程的系數(shù)矩陣求解，即為回代過程 */* p-指向方程系數(shù)矩陣； q-指向方程右端向量；*/double *x; int i,t,h1,h2

11、; double sum; x=(double *)malloc(sizeof(double)*(N+1); *(x+1)=q1; for(i=2;i<=N;i+) h1=i-s; if(h1<1) h1=1; for(t=h1,sum=0;t<=i-1;t+) sum+=*(*(p+i-t+3)+t)*xt; xi=qi-sum; *(x+N)=*(x+N)/(*(*(p+3)+N); for(i=N-1;i>=1;i-) h2=i+s; if(i+s>N) h2=N; for(t=i+1,sum=0;t<=h2;t+) sum+=*(*(p+i-t+3)

12、+t)*(*(x+t); *(x+i)=(*(x+i)-sum)/(*(*(p+3)+i); return x;double Tmin(double (*p)N+1) /* 求矩陣的模最小的特征值 ，采用反冪法 */ /* p指向矩陣 */ double yN+1,a,b1=1,b2,sum,t=1; double *u; int i; u=(double *)malloc(sizeof(double)*(N+1); for(i=1;i<=N;i+) ui=1; while(fabs(t)>1e-12) for(i=1,sum=0;i<=N;i+)sum+=ui*ui; a=

13、sqrt(sum); for(i=1;i<=N;i+) yi=ui/a; u=Hdai(p,y); for(i=1,b2=0;i<=N;i+) b2+=yi*ui;t=(1/b2-1/b1)*b2; b1=b2; return 1/b1;三、輸出結(jié)果四、對初始向量的討論在本例利用冪法求模最大特征值中，不同的初值對結(jié)果的影響比較大，例如若取U0=（1，1，···，1）時，求得max= -1.070011e+001；而當U0=（1，0，0，···，0）或U0=（0，1，0，···，0）時，max=

14、-2.080981，冪法好像“失效”了。下面分別是兩個不同初值的運行情況（次數(shù)是指迭代次數(shù)k，t表示每次迭代(k+1-k)/ k+1的值）：1、取U0=（1，1，···，1）次數(shù)=3,特征值=-1.366284 ,t=0.3725088058072次數(shù)=4,特征值=-1.862396 ,t=0.2663837322315······次數(shù)=340,特征值=-10.700114 ,t=0.0000000000012次數(shù)=341,特征值=-10.700114 ,t=0.0000000000011次數(shù)=342,特

15、征值=-10.700114 ,t=0.0000000000010次數(shù)=343,特征值=-10.700114 ,t=0.00000000000102、U0=（1，0，0，···，0）次數(shù)=2,特征值=-0.255820 ,t=0.7933720926240次數(shù)=3,特征值=-0.358333 ,t=0.2860831422079次數(shù)=4,特征值=-0.341909 ,t=-0.0480354199885······次數(shù)=157,特征值=-2.080981 ,t=0.0000000000015次數(shù)=158,特

16、征值=-2.080981 ,t=0.0000000000012次數(shù)=159,特征值=-2.080981 ,t=0.0000000000010次數(shù)=160,特征值=-2.080981 ,t=0.0000000000009對比可發(fā)現(xiàn),第二個初值由于運行的次數(shù)較少，而t很快就減少到10-12以下，從而較快跳出循環(huán)，使得輸出的值僅為-2.080981。若強制運行次數(shù)達到1400次以上時，輸出的答案就是-10.700114。結(jié)果如下：······次數(shù)=1396,特征值=-10.700113 ,t=0.0000000021784次數(shù)=1397,特

17、征值=-10.700113 ,t=0.0000000018852次數(shù)=1398,特征值=-10.700114 ,t=0.0000000016314次數(shù)=1399,特征值=-10.700114 ,t=0.0000000014119次數(shù)=1400,特征值=-10.700114 ,t=0.0000000012218次數(shù)=1401,特征值=-10.700114 ,t=0.0000000010574事實上，若x1,x2,···x501分別對應(yīng)于A特征值1，2，···，501的特征向量，對于初始向量U0，存在唯一的不全為零的1，2，·&

18、#183;·，501，st：U0=1x1+2x2+···+501x501；k= Uk 2;則k=yk-1T Uk =1x1T+(2/1)k-12x2T+···+（501/1）k-1501x501T* 11x1+2(2/1)k-12x2+···+501（501/1）k-1501x501/ k2;記Uk（2，···，501）= (2/1)k2x2+···+（501/1）k501x501；則k=1（1x1T+UTk-1）（1x1+Uk）/k2 =1 (21x1Tx1+1UTk-1x1+1x1T Uk +UTk-1Uk) /k2不妨假設(shè)取到這樣的初值U0，使得1絕對值很?。ㄉ踔恋扔诹悖?，而2，···，501絕對值都很大，此時可理解為U0幾乎與x1“垂直”了再考慮t（2，···，501）=(k+1-k)/ k+1 當k不是很大時，則t（2，···，501）(UTkUk+1 -UTk-1Uk)/ UTk-1Uk易知其關(guān)于i（i=2,3&