C语言实现方阵的LU分解

1、对 n 阶方阵 A 的 LU 分解实验目的掌握对 n 阶方阵 A 的 LU 分解方法，并用 C 语言编程实现；观察 L、U 矩阵的特征，并 在此基础上得出求解方程方程组 AX=b 的方法。实验原理ij n*nij设矩阵 A = (a )的各阶顺序主子式均不为零，即：a(i) 0(i = 1, 2, 3L, k) ，则由定理：设 n 阶方阵 A 的顺序主子式 D1 , D2 ,L, Dn-1 均不为零，则 A 的 LU 分解存在且唯一。可得：A 有分解式 a11a12La1n a21a22La2n MMM an1 1an 2Lann u11u12Lu1 n l211 u22L u2 n = MO

2、 LU（1）OM ln1ln 2L1 unn 对比等式左边和右边乘积矩阵 LU 的第 r 行主对角元（含主对角元）的对应元素，得rari = lrkukik =1(i = r,L, n; r = 1, 2,L, n).（2）再对比等式两边第 r 列主对角元以下（不含主对角元）的对应元素，得rair = likukrk =1(i = r +1,L, n; r = 1, 2,L, n -1).（3）当 r=1 时，由（2）和（3）式分别解出u1i = a1i(i = 1, 2,L, n),（4）li1 =ai1u11(i = 2, 3,L, n).（5）假设已求出 U 的第 1 至 r-1 行，L

3、 的第 1 至 r-1 列，由式（2）和式（3）分别得出计算 U的第 r 行、L 的第 r 列元素的计算公式：r -1uri = ari - lrkukik =1(i = r,L, n; r = 2, 3,Ln),（6）lirr -1air - lrk ukr= k =1urr(i = r +1,L, n; r = 2, 3,L n -1).（7）称由（4）（7）式所表示的矩阵分解为 Doolittle 分解，也称为 LU 分解。程序实现过程以下均为程序的关键部分：1、 输入矩阵 Afor(i=1;iN+1;i+)%控制行数printf(请输入矩阵第%d 行元素：n,i); for(j=1;j

4、N+1;j+)%控制列数scanf(%f,&Aij);%输入每行的矩阵元素nkk2、 输入矩阵 A 后先判断 A 的各阶顺序主子式（即 k0，若不成立则不能进行 LU 分解a ( k ) (k = 1, 2,L n ) ）是否全大于for(t=1,k=d11;tN+1;t+)%控制阶数if(d11=0)&判断矩阵第一行是否为 0，若为 0 则不分解break;eLsefor(i=t+1;iN+1;i+)%进行第 t 次（行优先）所有元素的消元m=dit/dtt;%计算行乘数 for(j=t+1;jN;j+)dqj=dqj-m*dtj;%计算每次不为 0 的元素的值m=0;%将 m 归零n=dt

5、+1t+1;%取主元 nk*=n;%计算各阶的顺序主子式3、 若可分解，进行 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+=Lrk*Ukj;%待求元素之前对应的行和列的元素和Urj=crj-x;%计算 U 的第 r 行元素%往下是对 L 的第 r 列进行分解for(i=r+1;iN+1;i+)%控制 L 的第 r 列的 i 行变化y=0;%初始化中间变量 y for(k=1;kr;k+)y+=

6、Lik*Ukr;%待求元素之前对应的行和列的元素和Lir=(cir-y)/Urr;%计算 L 的第 r 列元素实验结果及分析1、 运行程序后，输入待分解矩阵 A（4*4） 2100-3 -3-4-1213 A= ，得到如下结果： 123 4149-4 -132、 重新运行程序，输入待分解矩阵 B(4*4) 2141 3-121 B= ， 得到如下结果: 1232 5062 3、 从上述两个不同矩阵的分解，可以看出：1） 只有矩阵为方阵且其各阶顺序主子式不为 0（即为非奇异矩阵），才可以进行 LU分解，并且分解式是唯一的；否则，不能进行 LU 分解。2） 分解的 L 矩阵为单位下三角矩阵，U 矩

7、阵为上三角矩阵。4、 实验简单心得：1) C 语言中矩阵的存储同二维数组，均采用行优先原则；本程序调用了两个函数 IsA（）和 DisA（）分别判断 A 矩阵是否为奇异矩阵和分解 A 矩阵。2) 判断矩阵 A 是否为奇异矩阵时，若第一个元素为 0 则其必为奇异矩阵，进行每次 消元循环时，必须计算出每一行的行列式，而且必须设一个中间变量 n 记录每次 消元后的主对角元素，最后通过每个主对角元素连乘 k 判断是否为 0。3) 进行矩阵 LU 分解时，必须同时计算 U 的第 r 行和 L 的第 r 列，在计算具体行或列 式则需先计算出 U 的第 r 行，才能计算出 L 的第 r 列。4) 本程序没有考虑部分选主元，若考虑则必须在每次分解前判断主元进行行交换， 增加了程序的设计难度。5) 程序不能实现矩阵大小的实时改变，若要改变则必须在源程序代码中修改宏定义 的 N 值，是一个缺陷。6) 通过程序的设计增强了我编程的意识，采用模块化，若出现错误则只需修改某个 函