C语言实现方阵的LU分解.doc

对n阶方阵A的LU分解 实验目的 掌握对n阶方阵A的LU分解方法，并用C语言编程实现；观察L、U矩阵的特征，并在此基础上得出求解方程方程组AX=b的方法。 实验原理 设矩阵的各阶顺序主子式均不为零，即：，则由定理： 设n阶方阵A的顺序主子式均不为零，则A的LU分解存在且唯一。 可得：A有分解式 =LU （1） 对比等式左边和右边乘积矩阵LU的第r行主对角元（含主对角元）的对应元素，得 （2） 再对比等式两边第r列主对角元以下（不含主对角元）的对应元素，得 （3） 当r=1时，由（2）和（3）式分别解出 （4） （5） 假设已求出U的第1至r-1行，L的第1至r-1列，由式（2）和式（3）分别得出计算U的第r行、L的第r列元素的计算公式： （6） （7） 称由（4）~（7）式所表示的矩阵分解为Doolittle分解，也称为LU分解。 程序实现过程 以下均为程序的关键部分： 输入矩阵A for(i=1;iN+1;i++) %控制行数 { printf(请输入矩阵第%d行元素：\n,i); for(j=1;jN+1;j++) %控制列数 scanf(%f,A[i][j]); %输入每行的矩阵元素 } 输入矩阵A后先判断A的各阶顺序主子式（即）是否全大于0，若不成立则不能进行LU分解 for(t=1,k=d[1][1];tN+1;t++) %控制阶数 { if(d[1][1]==0) 判断矩阵第一行是否为0，若为0则不分解 break; eLse { for(i=t+1;iN+1;i++) %进行第t次（行优先）所有元素的消元 { m=d[i][t]/d[t][t]; %计算行乘数 for(j=t+1;jN;j++) d[q][j]=d[q][j]-m*d[t][j]; %计算每次不为0的元素的值 m=0; %将m归零 } n=d[t+1][t+1]; %取主元n k*=n; %计算各阶的顺序主子式 } } 若可分解，进行LU分解。L为单位下三角矩阵，U为上三角矩阵。 for(r=2;rN+1;r++) %控制U的行数（L的列数） { for(j=r;jN+1;j++) %控制U的第r行的j列变化 { x=0; %初始化中间变量x for(k=1;kr;k++) x+=L[r][k]*U[k][j]; %待求元素之前对应的行和列的元素和 U[r][j]=c[r][j]-x; %计算U的第r行元素 } %往下是对L的第r列进行分解 for(i=r+1;iN+1;i++) %控制L的第r列的i行变化 { y=0; %初始化中间变量y for(k=1;kr;k++) y+=L[i][k]*U[k][r]; %待求元素之前对应的行和列的元素和 L[i][r]=(c[i][r]-y)/U[r][r]; %计算L的第r列元素 } } 实验结果及分析 运行程序后，输入待分解矩阵A（4*4） A=， 得到如下结果： 重新运行程序，输入待分解矩阵B(4*4) B=， 得到如下结果: 从上述两个不同矩阵的分解，可以看出： 只有矩阵为方阵且其各阶顺序主子式不为0（即为非奇异矩阵），才可以进行LU分解，并且分解式是唯一的；否则，不能进行LU分