实验五矩阵的lu分解法，雅可比迭代法班级：学号：：实验五 矩阵的LU 分解法，雅可比迭代一、目的与要求：➢ 熟悉求解线性方程组的有关理论和方法；➢ 会编制列主元消去法、LU 分解法、雅可比及高斯—塞德尔迭代法德程序；➢ 通过实际计算，进一步了解各种方法的优缺点，选择合适的数值方法。

二、实验容：➢ 会编制列主元消去法、LU 分解法、雅可比及高斯—塞德尔迭代法德程序，进一步了解各种方法的优缺点。

三、程序与实例➢ 列主元高斯消去法算法：将方程用增广矩阵[A ∣b ]=(ij a )1n (n )+⨯表示1) 消元过程对k=1,2,…,n-1①选主元，找{}n ,,1k ,k i k +∈使得k ,i k a =ik a ni k max ≤≤ ②如果0a k ,i k =，则矩阵A 奇异，程序结束；否则执行③。

③如果k i k ≠，则交换第k 行与第k i 行对应元素位置，j i k j k a a ↔ j=k,┅,n+1④消元，对i=k+1, ┅,n 计算kk ik ik a a l /=对j=l+1, ┅,n+1计算kj ik ij ij a l a a -=2) 回代过程①若0=nn a ，则矩阵A 奇异，程序结束；否则执行②。

②nn n n n a a x /1,+=；对i=n-1, ┅,2,1,计算ii ni j j ij n i i a x a a x /11,⎪⎪⎭⎫ ⎝⎛-=∑+=+程序与实例程序设计如下：#include <iostream>#include <cmath>using namespace std;void disp(double** p,int row,int col){for(int i=0;i<row;i++){for(int j=0;j<col;j++)cout<<p[i][j]<<' ';cout<<endl;}}void disp(double* q,int n){cout<<"====================================="<<endl; for(int i=0;i<n;i++)cout<<"X["<<i+1<<"]="<<q[i]<<endl;cout<<"====================================="<<endl; }void input(double** p,int row,int col){for(int i=0;i<row;i++){cout<<"输入第"<<i+1<<"行:";for(int j=0;j<col;j++)cin>>p[i][j];}}int findMax(double** p,int start,int end){int max=start;for(int i=start;i<end;i++){if(abs(p[i][start])>abs(p[max][start]))max=i;}return max;}void swapRow(double** p,int one,int other,int col){double temp=0;for(int i=0;i<col;i++){temp=p[one][i];p[one][i]=p[other][i];p[other][i]=temp;}}bool dispel(double** p,int row,int col){for(int i=0;i<row;i++){int flag=findMax(p,i,row); //找列主元行号if(p[flag][i]==0) return false;swapRow(p,i,flag,col); //交换行for(int j=i+1;j<row;j++){double elem=p[j][i]/p[i][i]; //消元因子 p[j][i]=0;for(int k=i+1;k<col;k++){p[j][k]-=(elem*p[i][k]);}}}return true;}double sumRow(double** p,double* q,int row,int col){ double sum=0;for(int i=0;i<col-1;i++){if(i==row)continue;sum+=(q[i]*p[row][i]);}return sum;}void back(double** p,int row,int col,double* q){for(int i=row-1;i>=0;i--){q[i]=(p[i][col-1]-sumRow(p,q,i,col))/p[i][i]; }}int main(){cout<<"Input n:";int n; //方阵的大小cin>>n;double **p=new double* [n];for(int i=0;i<n;i++){p[i]=new double [n+1];}input(p,n,n+1);if(!dispel(p,n,n+1)){cout<<"奇异"<<endl;return 0;}double* q=new double[n];for(int i=0;i<n;i++)q[i]=0;back(p,n,n+1,q);disp(q,n);delete[] q;for(int i=0;i<n;i++)delete[] p[i];delete[] p;}1. 用列主元消去法解方程⎪⎩⎪⎨⎧=++-=++-=++035.3643x .5072x .1835x .2137.2623x .43712x 347x .1 1.1833.555x 2.304x 0.101x 321321321例2 解方程组⎪⎪⎪⎩⎪⎪⎪⎨⎧=++++=++++=++++=++++=++++-12.041.0F 1.02E 3.47D 1.04C 3.54B -6.301.0F2.01E 2.51D 4.04C 5.05B -8.531.0F 1.21E2.92D 1.46C3.53B -20.071.0F 1.10E4.48D 1.21C 4.93B -32.041.0F 1.55E5.66D 2.40C 8.77B 计算结果如下B=-1.461954C= 1.458125D=-6.004824E=-2.209018F= 14.719421➢矩阵直接三角分解法算法：将方程组A x=b 中的A分解为A=LU，其中L为单位下三角矩阵，U为上三角矩阵，则方程组A x=b化为解2个方程组Ly=b，Ux=y，具体算法如下：①对j=1,2,3,…,n计算jjau11=对i=2,3,…,n计算1111/aalii=②对k=1,2,3,…,n:a. 对j=k,k+1,…,n计算∑-=-=11kqqjkqkjkjulaub. 对i=k+1,k+2,…,n计算kkkqqkiqikikuulal/)(11∑-=-=③11by=，对k=2,3,…,n计算∑-=-=11kqqkqkkylby④nnnnuyx/=,对k=n-1,n-2,…,2,1计算kknkqqkqkkuxuyx/)(1∑+=-=注：由于计算u 的公式于计算y 的公式形式上一样，故可直接对增广矩阵[A ∣b ]=⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡+++1,211,2222211,111211n n nn n n n n n n a a a a a a a a a a a a 施行算法②，③，此时U 的第n+1列元素即为y 。

程序与实例例3 求解方程组A x=bA=⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡-----381265973274581221, b=⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡4911427 结果为X[0]= 3.000001X[1]=-2.000001X[2]= 1.000000X[3]= 5.000000#include <iostream>using namespace std;double** newMatrix(int row,int col){double **p=new double* [row]; //行for(int i=0;i<row;i++) //列p[i]=new double [col];return p;}void delMatrix(double** p,int row){for(int i=0;i<row;i++)delete[] p[i];delete[] p;}void inputMatrix(double** p,int row,int col){for(int i=0;i<row;i++){cout<<"输入第"<<i+1<<"行:";for(int j=0;j<col;j++)cin>>p[i][j];}}void dispMatrix(double** p,int row,int col){for(int i=0;i<row;i++){for(int j=0;j<col;j++)cout<<p[i][j]<<' ';cout<<endl;}}void dispVector(double* q,int n){cout<<"====================================="<<endl;for(int i=0;i<n;i++)cout<<"X["<<i+1<<"]="<<q[i]<<endl;cout<<"====================================="<<endl;}void initMatrix(double** p,int row,int col){for(int i=0;i<row;i++)for(int j=0;j<col;j++)p[i][j]=0;}double sum(double** L,double** U,int row,int col){int min=(row>=col?col:row);double sum=0;for(int i=0;i<min;i++)sum+=(U[i][col]*L[row][i]);return sum;}void resolve(double** A,double** L,double** U,int row,int col){ for(int i=0;i<col;i++)U[0][i]=A[0][i]; //把A的第一行给U的第一行L[0][0]=A[0][0];for(int i=1;i<row;i++) //填充L的第一列L[i][0]=A[i][0]/A[0][0];for(int i=1;i<row;i++)for(int j=1;j<col;j++){if(i<=j)U[i][j]=A[i][j]-sum(L,U,i,j);elseL[i][j]=(A[i][j]-sum(L,U,i,j))/U[j][j];}for(int i=0;i<row;i++)L[i][i]=1;}double sumRowY(double** L,double* y,int row){double sum=0;for(int i=0;i<row;i++)sum+=(L[row][i]*y[i]);return sum;}void backY(double** L,double* b,double* y,int row){for(int i=0;i<row;i++)y[i]=0; //初始化y向量全为零for(int i=0;i<row;i++)y[i]=b[i]-sumRowY(L,y,i);}double sumRowX(double** U,double* x,int row,int col){double sum=0;for(int i=row+1;i<col;i++)sum+=(U[row][i]*x[i]);return sum;}void backX(double** U,double* y,double* x,int row){for(int i=0;i<row;i++)x[i]=0; //初始化x向量全为零for(int i=row-1;i>=0;i--)x[i]=(y[i]-sumRowX(U,x,i,row))/U[i][i];}int main(){cout<<"输入方阵大小 n :";int n=0;cin>>n;cout<<"====================================="<<endl;cout<<"开始录入方阵数据....."<<endl;double** A=newMatrix(n,n); //开辟矩阵AinputMatrix(A,n,n); //录入数据到A中cout<<"录入方阵数据完毕......"<<endl;cout<<"====================================="<<endl<<endl;cout<<"开始录入列阵....."<<endl;cout<<"输入列阵:";double* b=new double[n]; //开辟向量bfor(int i=0;i<n;i++)cin>>b[i]; //录入数据到b中 cout<<"录入列阵数据完毕....."<<endl;double** L=newMatrix(n,n); //开辟方阵LinitMatrix(L,n,n); //初始化L全为零double** U=newMatrix(n,n); //开辟方阵UinitMatrix(U,n,n); //初始化Uresolve(A,L,U,n,n); //分解A为L与Udouble* y=new double[n]; //开辟向量ybackY(L,b,y,n); //回代求出ydouble* x=new double[n]; //开辟向量XbackX(U,y,x,n); //回代求出XdispVector(x,n);//释放空间delMatrix(A,n);delMatrix(L,n);delMatrix(U,n);delete[] b;delete[] y;delete[] x;}➢ 迭代法雅可比迭代法算法：设方程组Ax=b 系数矩阵的对角线元素)n ,,2,1i (0a ii =≠,M 为迭代次数容许的最大值,ε为容许误差。