北航研究生数值分析编程大作业

1、WORD格式整理版数值分析大作业一、算法设计方案1、矩阵初始化矩阵A二bj 501 501的下半带宽r=2,上半带宽s=2,设置矩阵C 'r s 1 fed 1, 在矩阵C中检索矩阵A中的带内元素讦的方法是：aj =cKj.这样所需要的 存储单元数大大减少,从而极大提升了运算效率.2、利用幕法求出,501幕法迭代格式：非零向量 u0 e Rn=JukUky k J - u k J / k _1u k = Ay k3 k = yku k'J当|悅-B|/Bk兰10上时,迭代终止.首先对于矩阵A利用幕法迭代求出一个,然后求出矩阵B,其中 B二A -1 ( |为单位矩阵),对矩阵B进

2、行幕法迭代,求出,之后令"', 比拟与,的大小,大者为'501,小者为 1.3、利用反幕法求出-s, ik反幕法迭代格式：'非零向量u0ERn4 = JUk 4Uk* yk=Uk/nk J-Auk =yz和=yLu当|丘-円/陈兰10孔时,迭代终止,丸s=1Rk.每迭代一次都要求解一次线性方程组 Auyk4,求解过程为：(1)作分解A二LU对于k =1,2,., n执行k 二cks.1,j：=cks 1,j -' Ck_t -s 1,tcts 1,jt=max(i,k_r,j _s)l.j = k, k 1,., min( k s, n) Ikci _

3、k s 1,k ： 二(ci _k来申,k 一 为 Ci_t衣*tct_k韦申,k )/Cs 1,kt maxOi_r ,k_s) = k 1,k2,., min( k - r, n); k : s 1(2)求解Ly二b,Ux =y (数组b先是存放原方程组右端向量,后来存放中间向量y)i 二bi 二 bi ci _k；s 1,tbt (i = 2,3,.,n)t max(1 ,i _r)xn =bn / cs 1.nmin(i -s,n)Xi = (bi -' Ci_t s 1.tXt)/Cs 1,i(i = n -1, n -2,.,1)t丄卅使用反幕法,直接可以求得矩阵按模最小的

4、特征值s O求与数1输二011(k =1,2,.,39)最接近的特征值 h,对矩阵A-SI40实行反幂法,即可求出对应的= 1/ k -k O4、求出A的条件数和行列式根据cond(A)2 =|甲彳,其中分子分母分别对应按模最大和最小的特征值.det(A)的计算：由于A二LU ,其中L为下三角矩阵,且对角线元素为1,故det(L) =1,所以有A = LU = U,又U为上三角矩阵,故det(U )为对其对角线上各元素的乘积,最后可得det(A)二det(U )学习参考好帮手二、程序源代码(1) 定义所需要的函数：#i nclude <stdio.h>#in elude vconi

5、 o.h>#in elude <math.h>#defi ne N 501#defi ne R 2#defi ne S 2int min (i nt a,i nt b); /求最小值int max(i nt a,i nt b,i nt c); /求最大值double Fan_two(double xN);计算二范数void FenjieLU(double (*C)N);解线性方程组的LU分解过程void Solve(double (*C)N, double *b,double *x);解线性方程组的求解过程double PowerMethod(double CN,double

6、 uN,double yN,double bta,double D);幕法double In versePowerMethod(double CN,double uN,doubleyN,double bta,double D);反幕法;(2) 程序的主函数,Main.cpp代码如下：void mai n()double CR+S+1N;double uN;double yN;double miu39;double C1R+S+1N;double bta = 1.0;double Namda1,Namda501,NamdaS;double Namda39;double Co ndA2;doubl

7、e detA = 1.0;double D = 1.0e-12;int i, j, k;FILE * fp;fp = fope n( "Namda.txt","w");/对数组进行初始化/int i, j;for (i = 0; i < N; i+)ui = 1;for (i = 0;i< R + S + 1;i+)for (j = 0;j< N;j+)if (i=0|i=4)Cij=-0.064;else if (i=1|i=3)Cij=0.16;else if (i=2)Cij=(1.64-0.024*(j+1)*si n(0.2*

8、(j+1) -0.64*exp(0.1/(j+1);/幕法求Namda1Namda1 = PowerMethod(C, u, y, bta, D);prin tf("n=n"); prin tf("Namda1 = %12.11e", Namdal);prin tf("n=n");/ 幕法求 Namda501bta = 1.0;for (i = 0; i < R + S + 1; i+)for (j = 0; j < N; j+)if (i = 2)C1ij =Cij -Namdal;elseC1ij = Cij;Namd

9、a501 = algorism.PowerMethod(C1, u, y, bta, D) +Namda1;prin tf("n=n");prin tf("Namda501 = %12.11e", Namda501);prin tf("n=n"); /反幕法求NamdaS/bta = 1.0;NamdaS = In versePowerMethod(C, u, y, bta, D);prin tf("n=n"); prin tf("NamdaS = %12.11e", NamdaS);prin

10、tf("n=n"); /反幕法求Namdakprin tf("n=n"); for (k = 0; k < 39; k+)miuk = Namda1 + (k + 1) * (Namda501 - Namda1) / 40.0;bta = 1.0;for (i = 0; i < R + S + 1; i+)for (j = 0; j < N; j+)if (i = 2)C1ij =Cij - miuk;elseC1ij = Cij;Namdak = In versePowerMethod(C1, u, y, bta, D) + miuk

11、; fprin tf(fp,"与%12.11e最接近的特征值为:12.11en",miuk,Namdak);printf(" 求与miuk最接近的Namdak的计算结果已经输出到文件Namda.txt 中");prin tf("n=n"); /求A的谱范数/prin tf("n=n");printf("A 的谱范数为:%12.11e", sqrt(Namda501);prin tf("n=n"); /求A的条件数/Con dA2 = fabs( Namda1 / NamdaS

12、);prin tf("n=n");printf("A 的谱范数的条件数 Cond(A)2为:%12.11e",CondA2);prin tf("n=n"); / 求det(A)2 的值/for (j = 0; j < N; j+)detA *= C2j;prin tf("n=n");printf("行列式 A的值为:%12.11e",detA);prin tf("n=n"); fclose(fp);_getch();return;(3)成员函数的实现int min (i

13、 nt a,i nt b)return a < b ? a : b;int max(i nt a,i nt b,i nt c)int temp;temp = a > b ? a : b;return temp > c ? temp : c;double Fan_two(double xN) _double sum = 0.0;int i;for (i = 0; i < N; i+)sum += pow(xi,2);return sqrt(sum);void FenjieLU(double (*C)N)double sum = 0;int i, j, k,t;for (k

14、 = 0; k < N; k+)j = k;i = k + 1;while (1)if (j = mi n(k + S + 1, N)break;for (t = max(0, k - R, j - S); t <= k - 1; t+) sum += Ck-t+St * Ct-j+Sj;Ck-j+Sj = Ck-j+Sj - sum;sum = 0.0;j+；if (k = N-1)break;if (i = mi n(k + R + 1, N) break;for (t = max(0, i - R,k - S); t <= k - 1; t+)sum += Ci-t+S

15、t * Ct-k+Sk;Ci-k+Sk = (Ci-k+Sk - sum) / CSk;sum = 0;i+;void Solve(double (*C)N, double *b,double *x)double sum = 0;int i, t;sum = 0;for (i = 1; i < N; i+)for (t = max(0, i - R); t <= i - 1; t+)sum += Ci-t+St * bt;bi = bi - sum;sum = 0;xN-1 = bN-1 / CSN-1;for (i = N - 2; i >= 0; i-)for (t =

16、i+1; t <= mi n(i + S, N - 1); t+)sum += Ci-t+St * xt;xi = (bi - sum) / CSi;sum = 0;double PowerMethod(double CN,double uN,double yN,double bta,double D)double ita;double sum = 0;double temp = 0.0;int i,j,k = 0;while (fabs(bta - temp) / fabs(bta) > D)temp = bta;ita = Fan _two(u);for (i = 0; i &

17、lt; N; i+)yi = ui / ita;for (i = 0; i < N; i+)for (j = max(0,i - R); j < mi n(i + S + 1,N); j+)sum += Ci - j + Sj * yj;ui = sum;sum = 0;for (i = 0; i < N; i+)sum += yi * ui;bta = sum;sum = 0;k+;return bta;double In versePowerMethod(double CN,double uN,double yN,double bta,double D)double TC

18、R+S+1N;double tyN;double ita;double sum = 0;double temp = 0.0;int i,j,k = 0;FenjieLU(C);while (abs(1/bta - 1/temp) / abs(1/bta) > D)temp = bta;ita = Fan _two(u);for (i = 0; i < N; i+)yi = ui / ita;/用到临时存储数组TC和ty是由于函数Solve执行过程中会改变 A和 yfor (i = 0; i < R + S + 1; i+)for (j = 0; j < N; j+)TC

19、ij = Cij;for (i = 0; i < N; i+)tyi = yi;Solve(C, y, u);for (i = 0; i < R+S+1; i+)for (j = 0; j < N; j+)Cij = TCij;for (i = 0; i < N; i+)yi = tyi;for (i = 0; i < N; i+)sum += yi * ui;bta = sum;sum = 0;k+;bta = 1.0 / bta;return bta;三、程序运行结果下列图为主程序运行结果其中'ik的结果输出在Namda.txt文件中,结果如下:与Y

20、018949*i9222e + 001> 与-9.67887622933e*000S： -9.16825753648e*00eft _-8 .65763S8436Ue + caoJr -ft_lH702 QISOSQe + OOsW -7-63640145795etO8Ow 4-7-1257fi276511et00fl® 4-&-&15164ft7226e + 0O0£ 与-6.1045537942e + 0O_ -5.59392668657e+O00jg； _-5.O833B799373e + 00Bft -4-572fi893 0O89e + O0

21、04. 06207 06 08 0Ue + 0Bfl< -3.551H519152Oe + 000<： 4-3 O083322235e + 00Oft -2.53 021452951e + 0O01 近的牛 近的 近的 近的 近的-2.01959583666e + O0Cb|； 写亠T. 50S977143«2e + O0GS； -9.98358U5 0975e-O0lW -4-87739758131e-001 2.28789347135e- 002> 与 $ ,33M7627558e-O0lfc 1 .C4411632 0t|0e+00R 与X .554735O1

22、325e+0O0M -2.06535370609? +000 _2.57597239894e + 000ft 与8-O86S91S9178e+BO&fe 3-59720978ii62e+O0&S 4>107«2fiV77e+O0B 4-618717 0319+001 与 5.125>06586316e+0OBft 5 _63968H-55600 +BOB% 6.1S0363248859 + 与血-660921 ?i*l69e+00&fe 7.1715UB63ii53e+00B 7.6«215932738e+00B 4e-192778 02

23、 022e+ 006 与8,70339671307e+OOft& -9_21£»8154O591e+000近的征566891I近的 近的 近的 近的 近的卡征等征近的恃征 近的轉 近的彳 近的 近的 近的匕 近的占 近的 近的 近舸 近的 遊的 近的 近的: 近的 近的 近的： 近的 近的占 近飾-1.01829340331+091 -9-5857O742S07e+008 -9+0OO ：-8-652284007909+000 :-fi _fl93483B0868e+0Q0 :-7.659U05iiG769e+ftO0 ：-7.11968a

24、»6Xi869e+00O ：-6.61176i»3394fle+00O ：-6.06610322660P+000 ：-5.585101 05263e*090 ：-5,11408352981e+000 ：->»-57887217687e+000 为:-.0978293G79Oe+flO0 ：-3.55421121575e+000：-3 . B4109 001813e+ 080 ：-2-526430591?+DOB ：-2,00323076956?+000 ：-1.0355761123?+000 ：-9 -4386060089-061 :-.73UU062771

25、1e"ftO1 2征 征 征 征 征 征 征 征 征 征 征 征 征 征 征 征22332317362U957e-002 B11&3921fi46e-0n OS289896269&+0O0 615786874956+ODO B6Q330U6027?+B00 6 034577ii527e*B0a 08024050931?+000 61362 0867«9«*O0O 09137851045e+000 C03 03537028+000+SCO 59490634&Ofle+OOO 08093385703P+000 6S035ii 09211?+00a 2938774iU813e+0O0 71711171ii2Jie+000 2