1、實(shí)驗(yàn)五矩陣的lu分解法，雅可比迭代法班級(jí)： 學(xué)號(hào)：姓名：實(shí)驗(yàn)五矩陣的LU分解法，雅可比迭代一、目的與要求： 熟悉求解線性方程組的有關(guān)理論和方法； 會(huì)編制列主元消去法、LU 分解法、雅可比及高斯塞德爾迭代法德程序； 通過實(shí)際計(jì)算，進(jìn)一步了解各種方法的優(yōu)缺點(diǎn)，選擇合適的數(shù)值方法。二、實(shí)驗(yàn)內(nèi)容： 會(huì)編制列主元消去法、LU 分解法、雅可比及高斯塞德爾迭代法德程序，進(jìn)一步了解各種方法的優(yōu)缺點(diǎn)。三、程序與實(shí)例 列主元高斯消去法算法：將方程用增廣矩陣Ab=(表示1) 消元過程對(duì)k=1,2,n-1選主元，找使得=如果，則矩陣A奇異，程序結(jié)束；否則執(zhí)行。如果，則交換第k行與第行對(duì)應(yīng)元素位置， j=k,n+1消元

2、，對(duì)i=k+1, ,n計(jì)算對(duì)j=l+1, ,n+1計(jì)算2) 回代過程若，則矩陣A奇異，程序結(jié)束；否則執(zhí)行。；對(duì)i=n-1, ,2,1,計(jì)算程序與實(shí)例程序設(shè)計(jì)如下：#include #include using namespace std;void disp(double* p,int row,int col) for(int i=0;irow;i+) for(int j=0;jcol;j+) coutpij ; coutendl; void disp(double* q,int n) cout=endl; for(int i=0;in;i+) coutXi+1=qiendl; cout=end

3、l;void input(double* p,int row,int col) for(int i=0;irow;i+) cout輸入第i+1行:; for(int j=0;jpij; int findMax(double* p,int start,int end) int max=start; for(int i=start;iabs(pmaxstart) max=i; return max;void swapRow(double* p,int one,int other,int col) double temp=0; for(int i=0;icol;i+) temp=ponei; pon

4、ei=potheri; potheri=temp; bool dispel(double* p,int row,int col) for(int i=0;irow;i+) int flag=findMax(p,i,row); /找列主元行號(hào) if(pflagi=0) return false; swapRow(p,i,flag,col); /交換行 for(int j=i+1;jrow;j+) double elem=pji/pii; /消元因子 pji=0; for(int k=i+1;kcol;k+) pjk-=(elem*pik); return true;double sumRow(d

5、ouble* p,double* q,int row,int col) double sum=0; for(int i=0;i=0;i-) qi=(picol-1-sumRow(p,q,i,col)/pii; int main() coutn; double *p=new double* n; for(int i=0;in;i+) pi=new double n+1; input(p,n,n+1); if(!dispel(p,n,n+1) cout奇異endl; return 0;double* q=new doublen; for(int i=0;in;i+) qi=0; back(p,n,

6、n+1,q); disp(q,n); delete q; for(int i=0;in;i+) delete pi; delete p; 1. 用列主元消去法解方程例2 解方程組計(jì)算結(jié)果如下B=-1.461954C= 1.458125D=-6.004824E=-2.209018F= 14.719421 矩陣直接三角分解法算法：將方程組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

7、計(jì)算，對(duì)k=2,3,n計(jì)算,對(duì)k=n-1,n-2,2,1計(jì)算注：由于計(jì)算u的公式于計(jì)算y的公式形式上一樣，故可直接對(duì)增廣矩陣Ab=施行算法，此時(shí)U的第n+1列元素即為y。程序與實(shí)例例3 求解方程組Ax=bA=, b=結(jié)果為X0= 3.000001X1=-2.000001X2= 1.000000X3= 5.000000#include using namespace std;double* newMatrix(int row,int col) double *p=new double* row; /行 for(int i=0;irow;i+) /列 pi=new double col; retu

8、rn p;void delMatrix(double* p,int row) for(int i=0;irow;i+) delete pi; delete p; void inputMatrix(double* p,int row,int col) for(int i=0;irow;i+) cout輸入第i+1行:; for(int j=0;jpij; void dispMatrix(double* p,int row,int col) for(int i=0;irow;i+) for(int j=0;jcol;j+) coutpij ; coutendl; void dispVector(d

9、ouble* q,int n) cout=endl; for(int i=0;in;i+) coutXi+1=qiendl; cout=endl;void initMatrix(double* p,int row,int col) for(int i=0;irow;i+) for(int j=0;j=col?col:row); double sum=0; for(int i=0;imin;i+) sum+=(Uicol*Lrowi); return sum;void resolve(double* A,double* L,double* U,int row,int col) for(int i

10、=0;icol;i+) U0i=A0i; /把A的第一行給U的第一行 L00=A00; for(int i=1;irow;i+) /填充L的第一列 Li0=Ai0/A00; for(int i=1;irow;i+) for(int j=1;jcol;j+) if(i=j) Uij=Aij-sum(L,U,i,j); else Lij=(Aij-sum(L,U,i,j)/Ujj; for(int i=0;irow;i+) Lii=1;double sumRowY(double* L,double* y,int row) double sum=0; for(int i=0;irow;i+) sum

11、+=(Lrowi*yi); return sum;void backY(double* L,double* b,double* y,int row) for(int i=0;irow;i+) yi=0; /初始化y向量全為零 for(int i=0;irow;i+) yi=bi-sumRowY(L,y,i);double sumRowX(double* U,double* x,int row,int col) double sum=0; for(int i=row+1;icol;i+) sum+=(Urowi*xi); return sum;void backX(double* U,doubl

12、e* y,double* x,int row) for(int i=0;i=0;i-) xi=(yi-sumRowX(U,x,i,row)/Uii;int main() coutn; cout=endl; cout開始錄入方陣數(shù)據(jù).endl; double* A=newMatrix(n,n); /開辟矩陣A inputMatrix(A,n,n); /錄入數(shù)據(jù)到A中 cout錄入方陣數(shù)據(jù)完畢.endl; cout=endlendl; cout開始錄入列陣.endl; cout輸入列陣:; double* b=new doublen; /開辟向量b for(int i=0;ibi; /錄入數(shù)據(jù)到b

13、中 cout錄入列陣數(shù)據(jù)完畢.endl; double* L=newMatrix(n,n); /開辟方陣L initMatrix(L,n,n); /初始化L全為零 double* U=newMatrix(n,n); /開辟方陣U initMatrix(U,n,n); /初始化U resolve(A,L,U,n,n); /分解A為L與U double* y=new doublen; /開辟向量y backY(L,b,y,n); /回代求出y double* x=new doublen; /開辟向量X backX(U,y,x,n); /回代求出X dispVector(x,n); /釋放空間 de

14、lMatrix(A,n); delMatrix(L,n); delMatrix(U,n); delete b; delete y; delete x; 迭代法雅可比迭代法算法：設(shè)方程組Ax=b系數(shù)矩陣的對(duì)角線元素,M為迭代次數(shù)容許的最大值,為容許誤差。取初始向量x=，令k=0。對(duì)i=1,2,n 計(jì)算如果，則輸出，結(jié)束；否則執(zhí)行。如果kM,則不收斂，終止程序；否則,轉(zhuǎn)。程序與實(shí)例例4 用雅可比迭代法解方程組結(jié)果為迭代次數(shù)為20X0= 1.000000X1= 2.000000X2=-1.000000#include #include using namespace std;double* newM

15、atrix(int row,int col) double *p=new double* row; /行 for(int i=0;irow;i+) /列 pi=new double col; return p;void delMatrix(double* p,int row) for(int i=0;irow;i+) delete pi; delete p; void inputMatrix(double* p,int row,int col) for(int i=0;irow;i+) cout輸入第i+1行:; for(int j=0;jpij; void dispMatrix(double

16、* p,int row,int col) for(int i=0;irow;i+) for(int j=0;jcol;j+) coutpij ; coutendl; void dispVector(double* q,int n) for(int i=0;in;i+) coutXi+1=qit; coutendl;double sumRow(double* A,double* x1,int row,int col) double sum=0; for(int i=0;icol;i+) if(row=i) continue; sum+=(Arowi*x1i); return sum;void i

17、terat(double* A,double* b,double* x1,double* x2,int row) for(int i=0;irow;i+) x2i=(bi-sumRow(A,x1,i,row)/Aii;bool blow_error(double* x1,double* x2,int row,double e) double sum=0; for(int i=0;irow;i+) sum+=abs(x2i-x1i); if(sume) return true; else return false;int main() coutn; cout=endl; cout開始錄入方陣數(shù)據(jù)

18、.endl; double* A=newMatrix(n,n); /開辟矩陣A inputMatrix(A,n,n); /錄入數(shù)據(jù)到A中 cout錄入方陣數(shù)據(jù)完畢.endl; cout=endlendl; cout開始錄入列陣.endl; cout輸入列陣:; double* b=new doublen; /開辟向量b for(int i=0;ibi; /錄入數(shù)據(jù)到b中 cout錄入列陣數(shù)據(jù)完畢.endl; cout=endl; coutM; coute; cout=endl; double* x1=new doublen; /開辟x1向量 double* x2=new doublen; /開

19、辟x2向量 cout輸入初始向量:; for(int i=0;ix1i; cout=endl; int k=0; /迭代計(jì)數(shù)器 for(;) iterat(A,b,x1,x2,n); /迭代一次 if(blow_error(x1,x2,n,e) dispVector(x2,n); cout迭代次數(shù)為:k=M) cout不收斂endl; return 0; else k+; for(int i=0;in;i+) x1i=x2i; /釋放空間 delMatrix(A,n); delete b; delete x1; delete x2; 高斯-塞爾德迭代法算法:設(shè)方程組Ax=b的系數(shù)矩陣的對(duì)角線元

20、素,M為迭代次數(shù)容許的最大值,為容許誤差取初始向量,令k=0。對(duì)i=1,2,n計(jì)算如果，則輸出,結(jié)束；否則執(zhí)行。如果則不收斂，終止程序；否則,轉(zhuǎn)。程序與實(shí)例例5 用高斯-塞爾德迭代法解方程組結(jié)果為X0=3.000000X1=2.000000X2=1.000000#include #include #include #include #define N 100main()int i;float *x;float c12=8.0,-3.0,2.0,20.0,4.0,11.0,-1.0,33.0,6.0,3.0,12.0,36.0;float *GauseSeidel(float *,int);x=

21、GauseSeidel(c,3);for (i=0;i=2;i+) printf(x%d=%fn,i,xi);getch();float *GauseSeidel(float *a,int n)int i,j,nu=0;float *x,dx,d;x=(float *)malloc(n*sizeof(float);for (i=0;i=n-1;i+) xi=0.0;dofor (i=0;i=n-1;i+)d=0.0;for (j=0;j=N)printf(fold numbern);nu+;while (fabs(dx)1e-6);return x;例6 用雅可比迭代法解方程組迭代4次得解,若

22、用高斯-塞爾德迭代法則發(fā)散。#includevoid main (void)int k,n;double x3=7,2,5;for(k=0;k5;k+)double a,b;a=x0;b=x1;x0=(7-2.0*x1+2*x2)/1;x1=(2-a-x2)/1;x2=(5-2*a-2*b)/1;for(n=0;n3;n+)printf(x%d=%8.6fn,n,xn); 用高斯-塞爾德迭代法解方程組迭代84次得解，若用雅克比迭代法則發(fā)散。#include#includevoid LOOP(float a1010,float b10,float x10,int);void main(void)

23、 float a1010,b10,x10,A10;double S;int M,n,i,j;printf(請(qǐng)輸入方陣階數(shù):);scanf(%d,&n);printf(請(qǐng)輸入最大允許迭代次數(shù):);scanf(%d,&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+)scanf(%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=0&S=0.000001);if(M=0)printf(迭代次數(shù)M=%dn,M);for(i=0;i=n-1;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;