河北工業(yè)大學(xué)數(shù)值分析實(shí)驗(yàn)五實(shí)驗(yàn)六實(shí)驗(yàn)報(bào)告

1、 數(shù)值分析實(shí)驗(yàn)報(bào)告 指導(dǎo)老師： 宛艷萍 姓 名： 班 級(jí)： 學(xué) 號(hào)：實(shí)驗(yàn)五矩陣的LU分解法，雅可比迭代1. 實(shí)驗(yàn)名稱(chēng) 矩陣的LU分解法，雅可比迭代2. 實(shí)驗(yàn)?zāi)康?1）熟悉求解線性方程組的有關(guān)理論和方法； 2）會(huì)編制列主元消去法、LU 分解法、雅可比及高斯塞德?tīng)柕ǖ鲁绦颍?3）通過(guò)實(shí)際計(jì)算，進(jìn)一步了解各種方法的優(yōu)缺點(diǎn)，選擇合適的數(shù)值方法。3. 實(shí)驗(yàn)內(nèi)容會(huì)編制列主元消去法、LU 分解法、雅可比及高斯塞德?tīng)柕ǖ鲁绦?，進(jìn)一步了解各種方法的優(yōu)缺點(diǎn)。4. 算法描述 1）列主元高斯消去法算法：將方程用增廣矩陣Ab=(表示1) 消元過(guò)程對(duì)k=1,2,n-1選主元，找使得=如果，則矩陣A奇異，程序結(jié)束；

2、否則執(zhí)行。如果，則交換第k行與第行對(duì)應(yīng)元素位置， j=k,n+1消元，對(duì)i=k+1, ,n計(jì)算對(duì)j=l+1, ,n+1計(jì)算2) 回代過(guò)程若，則矩陣A奇異，程序結(jié)束；否則執(zhí)行。；對(duì)i=n-1, ,2,1,計(jì)算 2）矩陣直接三角分解法算法：將方程組Ax=b 中的A分解為A=LU，其中L為單位下三角矩陣，U為上三角矩陣，則方程組Ax=b化為解2個(gè)方程組Ly=b，Ux=y，具體算法如下：對(duì)j=1,2,3,n計(jì)算對(duì)i=2,3,n計(jì)算對(duì)k=1,2,3,n:a. 對(duì)j=k,k+1,n計(jì)算b. 對(duì)i=k+1,k+2,n計(jì)算，對(duì)k=2,3,n計(jì)算,對(duì)k=n-1,n-2,2,1計(jì)算注：由于計(jì)算u的公式于計(jì)算y的公

3、式形式上一樣，故可直接對(duì)增廣矩陣Ab=施行算法，此時(shí)U的第n+1列元素即為y。 3）雅可比迭代法算法：設(shè)方程組Ax=b系數(shù)矩陣的對(duì)角線元素,M為迭代次數(shù)容許的最大值,為容許誤差。取初始向量x=，令k=0。對(duì)i=1,2,n 計(jì)算如果，則輸出，結(jié)束；否則執(zhí)行。如果kM,則不收斂，終止程序；否則,轉(zhuǎn)。 4）高斯-塞爾德迭代法算法:設(shè)方程組Ax=b的系數(shù)矩陣的對(duì)角線元素,M為迭代次數(shù)容許的最大值,為容許誤差取初始向量,令k=0。對(duì)i=1,2,n計(jì)算如果，則輸出,結(jié)束；否則執(zhí)行。如果則不收斂，終止程序；否則,轉(zhuǎn)。5. 程序與實(shí)例 1）列主元高斯消去法源程序： 1）例1#include<stdio.

4、h>#include<math.h>void main(void)float x10,A10,B,amax,S;int n=3,m,i,j,max; float a33=0.101,2.304,3.555,-1.347, 3.712 ,4.623,-2.835 ,1.072, 5.643;float b3=1.183,2.137,3.035;for(i=0;i<=n-1;i+) /*求上三角矩陣*/amax=aii;max=i;for(m=i+1;m<=n-1;m+)if(fabs(ami)>fabs(amax)amax=ami;max=m;for(m=0;

5、m<=n-1;m+)Am=aim;aim=amaxm;amaxm=Am;B=bi;bi=bmax;bmax=B;for(m=i+1;m<=n-1;m+)Am=ami/aii;for(m=i+1;m<=n-1;m+)for(j=i;j<=n-1;j+)amj=amj-Am*aij; bm=bm-Am*bi;xn-1=bn-1/an-1n-1; /*求x*/for(i=n-2;i>=0;i-) S=0.0;for(j=n-1;j>=i+1;j-)S=S+aij*xj;xi=(bi-S)/aii;for(i=0;i<=n-1;i+) /*輸出結(jié)果*/prin

6、tf("x%d=%10fn",i,xi); 2）例2#include<stdio.h>#include<math.h>void main(void)float x10,A10,B,amax,S;int n=5,m,i,j,max; float a55=8.77,2.40,5.66,1.55,1.0,4.93,1.21,4.48,1.10,1.0,3.53,1.46,2.92,1.21,1.0,5.05,4.04,2.51,2.01,1.0,3.54,1.04,3.47,1.02,1.0;float b5=-32.04,-20.07,-8.53,-6.

7、30,-12.04;for(i=0;i<=n-1;i+) /*求上三角*/amax=aii;max=i;for(m=i+1;m<=n-1;m+)if(fabs(ami)>fabs(amax)amax=ami;max=m;for(m=0;m<=n-1;m+)Am=aim;aim=amaxm;amaxm=Am;B=bi;bi=bmax;bmax=B;for(m=i+1;m<=n-1;m+)Am=ami/aii;for(m=i+1;m<=n-1;m+)for(j=i;j<=n-1;j+)amj=amj-Am*aij; bm=bm-Am*bi;xn-1=bn-

8、1/an-1n-1; /*求x*/for(i=n-2;i>=0;i-) S=0.0;for(j=n-1;j>=i+1;j-)S=S+aij*xj;xi=(bi-S)/aii;for(i=0;i<=n-1;i+) /*輸出*/printf("x%d=%10fn",i,xi);例1 解方程組輸出結(jié)果如下：X0=-0.398234X1= 0.013795X2= 0.335144運(yùn)行結(jié)果：例2 解方程組計(jì)算結(jié)果如下B=-1.161954C= 1.458125D=-6.004824E=-2.209018F= 14.719421 運(yùn)行結(jié)果：2） 矩陣直接三角分解法源程

9、序#include<stdio.h>int min(int,int);void main(void)float A1010,B1010,L1010,U1010,y10,b10,x10,S,EX;int i,j,c,k,n;printf("請(qǐng)輸入方陣的階數(shù):");scanf("%d",&n);printf("請(qǐng)按行輸入要分解的方程組各系數(shù)");for(i=0;i<=n-1;i+)for(j=0;j<=n-1;j+)scanf("%f",&Aij);printf("請(qǐng)依

10、次輸入方程組各值");for(i=0;i<=n-1;i+)scanf("%f",&bi);for(j=0;j<=n-1;j+)B0j=A0j;for(i=1;i<=n-1;i+)Bi0=Ai0/A00;for(i=1;i<=n-1;i+)for(j=1;j<=n-1;j+)c=min(i,j);EX=0;for(k=0;k<=c-1;k+)EX=EX+Bik*Bkj;Bij=Aij-EX;if(i>j)Bij=Bij/Bjj;printf("LU聯(lián)合矩陣為:n");for(i=0;i<=n

11、-1;i+)for(j=0;j<=n-1;j+)printf("%10f",Bij);if(j=n-1)printf("n"); printf("L矩陣為:n");for(i=0;i<=n-1;i+)for(j=0;j<=n-1;j+)if(i<j)Lij=0;else if(i=j)Lij=1;elseLij=Bij;for(i=0;i<=n-1;i+)for(j=0;j<=n-1;j+)printf("%10f",Lij);if(j=n-1)printf("n&qu

12、ot;);printf("U矩陣為:n");for(i=0;i<=n-1;i+)for(j=0;j<=n-1;j+)if(i>j)Uij=0;elseUij=Bij;for(i=0;i<=n-1;i+)for(j=0;j<=n-1;j+)printf("%10f",Uij);if(j=n-1)printf("n");y0=b0/L00; /*求y*/for(i=1;i<=n-1;i+)S=0.0;for(j=0;j<=i-1;j+)S=S+Lij*yj;yi=(bi-S)/Lii;xn-1=y

13、n-1/Un-1n-1; /*求x*/for(i=n-2;i>=0;i-)S=0.0;for(j=n-1;j>=i+1;j-)S=S+Uij*xj;xi=(yi-S)/Uii;for(i=0;i<=n-1;i+)printf("x%d=%4fn",i,xi);int min(int i,int j)int c;c=i;if(j<i)c=j;return(c);例3 求解方程組Ax=bA=, b= 結(jié)果為X0= 3.000001X1=-2.000001X2= 1.000000X3= 5.000000 運(yùn)行結(jié)果：3） 雅可比迭代法源程序#include&

14、lt;stdio.h>#include<math.h>void LOOP(float x10,float a10,float b10,int);void main(void)float a1010,b10,x10,A10;double S;int M,n,i,j,N;printf("請(qǐng)輸入方陣階數(shù):");scanf("%d",&n);printf("請(qǐng)輸入最大允許迭代次數(shù):");scanf("%d",&M);N=M;printf("請(qǐng)按行輸入各方程系數(shù):");f

15、or(i=0;i<=n-1;i+)for(j=0;j<=n-1;j+)scanf("%f",&aij);printf("請(qǐng)輸入各方程值:");for(i=0;i<=n-1;i+)scanf("%f",&bi);printf("請(qǐng)依次輸入首次迭代x值:");for(i=0;i<=n-1;i+)scanf("%f",&xi);dofor(i=0;i<=n-1;i+)Ai=xi;LOOP(x,a,b,n);M-;S=0.0;for(i=0;i<

16、;=n-1;i+)S=S+fabs(xi-Ai);while(M>=0&&S>=0.000001); if(M>=0) printf("迭代次數(shù)M=%dn",N); for(i=0;i<=n-1;i+) printf("x%d=%fn",i,xi);elseprintf("該迭代發(fā)散n");void LOOP(float x10,float a1010,float b10,int n)float S,A10;int i,j;for(i=0;i<=n-1;i+)Ai=xi;for(i=0;i

17、<=n-1;i+)S=0.0;for(j=0;j<=n-1;j+)if(j=i)continue;elseS=S+aij*Aj;xi=(bi-S)/aii;例4 用雅可比迭代法解方程組結(jié)果為迭代次數(shù)為20X0= 1.000000X1= 2.000000X2=-1.000000運(yùn)行結(jié)果：4） 高斯-塞爾德迭代法源程序#include<stdio.h>#include<math.h>void LOOP(float a1010,float b10,float x10,int);void main(void) float a1010,b10,x10,A10;doub

18、le S;int M,n,i,j,N;printf("請(qǐng)輸入方陣階數(shù):");scanf("%d",&n);printf("請(qǐng)輸入最大允許迭代次數(shù):");scanf("%d",&M);N=M;printf("請(qǐng)按行輸入各方程系數(shù):");for(i=0;i<=n-1;i+)for(j=0;j<=n-1;j+)scanf("%f",&aij);printf("請(qǐng)輸入各方程值:");for(i=0;i<=n-1;i+)sc

19、anf("%f",&bi);printf("請(qǐng)依次輸入首次迭代x值:");for(i=0;i<=n-1;i+)scanf("%f",&xi);doS=0.0;for(i=0;i<=n-1;i+)Ai=xi;LOOP(a,b,x,n);M-;for(i=0;i<=n-1;i+)S=S+fabs(xi-Ai);while(M>=0&&S>=0.000001);if(M>=0)printf("迭代次數(shù)M=%dn",N);for(i=0;i<=n-1

20、;i+)printf("x%d=%fn",i,xi);elseprintf("該迭代發(fā)散n");void LOOP(float a1010,float b10,float x10,int n)float S1,S2,A10;int i,j;for(i=0;i<=n-1;i+)Ai=xi;for(i=0;i<=n-1;i+)S1=0.0;S2=0.0;for(j=0;j<=i-1;j+)S1=S1+aij*xj;for(j=i+1;j<=n-1;j+)S2=S2+aij*Aj;xi=(bi-S1-S2)/aii;例5 用高斯-塞爾德

21、迭代法解方程組 結(jié)果為X0=3.000000X1=2.000000X2=1.000000 運(yùn)行結(jié)果： 例6 用雅可比迭代法解方程組迭代4次得解,若用高斯-塞爾德迭代法則發(fā)散。運(yùn)行結(jié)果： 1）雅可比迭代法： 2）高斯-塞爾德迭代法：例7 用高斯-塞爾德迭代法解方程組迭代84次得解，若用雅克比迭代法則發(fā)散。運(yùn)行結(jié)果： 1）高斯-塞爾德迭代法： 2）雅克比迭代法：實(shí)驗(yàn)六矩陣的特征值與特征向量的計(jì)算一、 目的與要求：Ø 領(lǐng)會(huì)求矩陣特征值及特征向量的冪法的理論及其方法；Ø 會(huì)編制冪法的計(jì)算程序，并用來(lái)計(jì)算有關(guān)問(wèn)題。二、 實(shí)驗(yàn)內(nèi)容：Ø 編制冪法的計(jì)算程序，并用來(lái)計(jì)算有關(guān)問(wèn)題三、 算法概要Ø 冪法是矩陣主特征值的一種迭代方法。設(shè)AR有n個(gè)線性無(wú)關(guān)的特征向量x，x,x,而相應(yīng)的特征值滿足>,則對(duì)任意非零初始向量Ø 用冪法求矩陣的按模最大的特征值和相應(yīng)的特征向量。精確至6位有效數(shù)字。取。Ø 結(jié)果為A按模取最大的特征值為7.00000；相應(yīng)的特征向量為。四、 源程序#include<stdio.h>#include<math.h>void LOOP(float a55,float u5,int);float MAX(float u5,int);void main(void)float a55,u5,x5,