选主元的三角分解法

2012-2013（1）专业课程实践论文选主元（列）的三角分解法范俊，0818180124，R数学08-1班毛勇，0818180117，R数学08-1班一、算法理论从直接三角分解公式可看出，当时计算将中断或者当绝对值很小时，按分解公式计算可能要求舍入误差的累积，但如果非奇异，就可通过交换的行实现矩阵的分解，因此可采用与列主元素消去法类似的方法（可以证明下述方法与列主元素消去法等价），将直接三角分解法修改为（部分）选主元的三角分解法。设第步分解已完成，这时有第步分解需要用到式及式为了避免用小的数作除数，引进量于是有,,用作为，交换的行与行元素(将位置的新元素仍记作及)，于是有。由此在进行第步分解计算。该程序即利用选主元的三角分解法解方程求方程的根。先选择列主元，再构造LU矩阵，通过求解LY=B和UX=Y得出需要的解。注：方程维数在程序中需按题意自定义。二、算法框图输入方程组选列主元是否是该列最大绝对值交换行构造LU矩阵求解LY=B求解UX=Y否是三、算法程序源代码源代码：LU_Decomposition.cpp#includeiostream.h#includemath.h#defineN4//矩阵维数，可自定义staticdoubleA[N][N];//系数矩阵staticdoubleB[N];//右端项staticdoubleY[N];//中间项staticdoubleX[N];//输出staticdoubleS[N];//选取列主元的比较器inti,j,k;//计数器voidmain(){cout请输入线性方程组（ai1，ai2，ai3......ain,yi）：endl;for(i=0;iN;i++){for(intj=0;jN;j++)cinA[i][j];cinB[i];}for(k=0;kN;k++){//选列主元intindex=k;for(i=k;iN;i++){doubletemp=0;for(intm=0;mk;m++){temp=temp+A[i][m]*A[m][k];}S[i]=A[i][k]-temp;if(S[index]S[i]){index=i;}}//交换行doubletemp;for(i=k;iN;i++){temp=A[index][i];A[index][i]=A[k][i];A[k][i]=temp;}temp=B[index];B[index]=B[k];B[k]=temp;//构造L、U矩阵for(j=k;jN;j++){doubletemp=0;for(intm=0;mk;m++){temp=temp+A[k][m]*A[m][j];}A[k][j]=A[k][j]-temp;//先构造U一行的向量}for(i=k+1;iN;i++){doubletemp=0;for(intm=0;mk;m++){temp=temp+A[i][m]*A[m][k];}A[i][k]=(A[i][k]-temp)/A[k][k];//再构造L一列的向量}}//求解LY=BY[0]=B[0];for(i=1;iN;i++){doubletemp=0;for(intj=0;ji;j++){temp=temp+A[i][j]*Y[j];}Y[i]=B[i]-temp;}//求解UX=YX[N-1]=Y[N-1]/A[N-1][N-1];for(i=N-2;i=0;i--){doubletemp=0;for(intj=i+1;jN;j++){temp=temp+A[i][j]*X[j];}X[i]=(Y[i]-temp)/A[i][i];}//打印Xcout线性方程组的解（X1，X2，X3......Xn）为：endl;for(i=0;iN;i++){coutX[i];}}四、算法实现例1.用选主元（列）的三角分解法解解：利用软件的解答方法为：（1）先将程序中的维数N改为3。（2）输入12314回车，输入25218回车，输入31520回车。（3）显示结果123例2.用选主元（列）的三角分解法解解：利用软件的解答方法为：（1）先将程序中的维数N改为3。（2）输入0341回车，输入1-112回车，输入2123回车。（3）显示结果1.27273-1.181820.818182