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分解。程序实现过程 以下均为程序的关键部分：1、 输入矩阵Afor(i=1;iN+1;i+) %控制行数printf(请输入矩阵第%d行元素：n,i);for(j=1;jN+1;j+) %控制列数scanf(%f,&Aij); %输入每行的矩阵元素2、 输入矩阵A后先判断A的各阶顺序主子式（即）是否全大于0，若不成立则不能进行LU分解for(t=1,k=d11;tN+1;t+) %控制阶数 if(d11=0) &判断矩阵第一行是否为0，若为0则不分解 break; eLse for(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+1t+1; %取主元n k*=n; %计算各阶的顺序主子式 3、 若可分解，进行LU分解。L为单位下三角矩阵，U为上三角矩阵。for(r=2;rN+1;r+) %控制U的行数（L的列数）for(j=r;jN+1;j+) %控制U的第r行的j列变化 x=0; %初始化中间变量xfor(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+=Lik*Ukr; %待求元素之前对应的行和列的元素和Lir=(cir-y)/Urr; %计算L的第r列元素实验结果及分析1、 运行程序后，输入待分解矩阵A（4*4） A=， 得到如下结果：2、 重新运行程序，输入待分解矩阵B(4*4)B=， 得到如下结果:3、 从上述两个不同矩阵的分解，可以看出：1） 只有矩阵为方阵且其各阶顺序主子式不为0（即为非奇异矩阵），才可以进行LU分解，并且分解式是唯一的；否则，不能进行LU分解。2） 分解的L矩阵为单位下三角矩阵，U矩阵为上三角矩阵。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) 通过程序的设计增强了我编程的意识，采用模块化，若出现错误则只需修改某个函数即可，及时修改一些隐蔽错误，会使程序的修改比较容易。参考文献1.计算机数值方