基于C语言实现大规模节点方程的LU分解算法.doc

基于C语言实现大规模节点方程的LU分解算法班级 : 0403205班学号 : 040320510姓名 : 邵汉钦完成日期 : 2006-03-25基于C语言实现大规模节点方程的LU分解算法实验目的及主要内容通过编程实现基于C语言的大规模节点方程的LU分解算法，熟悉在电路仿真中运用的节点分析法计算机实现过程，对算法的收敛性、效率、计算过程中节省内存的手段和重要性有清晰的认识，从而能更好地运用商用大规模仿真软件或自己动手编写仿真软件。实验原理对复杂网络所建立的节点方程组，常常是高阶线性方程组，所以说，对线性电路的稳态分析，最终可归结为高阶线性方程组得求解问题。线性方程组的解法很多，一般可以分为两大类：直接法和迭代法。直接法是通过对方程组直接进行一系列的运算而得到方程的解；迭代法是一开始就给出一个近似的解，然后通过逐次逼近来求LAE较好的近似解，迭代过程要一直进行到近似解达到预定的精度为止（即收敛到精确解）。直接法的两个基本方法是高斯消元法和LU分解法，它们也是稀疏矩阵技术的基础。LU分解是矩阵三角化分解法之一，也称三角分解法。它是把方程中的系数矩阵分解成下三角矩阵L和上三角矩阵U的乘积A=LU式中，L为主对角线以上的元素均为零的下三角矩阵；U为主对角线以下的元素均为零，且主对角线元素均为1的上三角矩阵。上式展开就是：=则方程变为： AX=LUX=B一旦A分解成功，我们就可以设UX=Y，方程就变成一个下三角形系统LY=B，解此方程可以用所谓“前代”的方法，即从第1个方程求出，进而从第2个方程求出，依次类推，直到求出，把，代入第n个方程求出。然后求解上三角系统UX=Y，这样就可以用回代的方法求出X来。LU分解法的优点主要表现在当系数矩阵A不变，即只是网络元件参数 不变，仅仅是右手边向量B变化时，即外加激励信号变化时，可以减少解方程的计算工作量。因为一旦系数矩阵分解成L和U的乘积，就可以把L，U存储起来，每次右手边向量变化时，只需做前代和回代，即可完成方程的求解。它在灵敏度计算中需要求转置的方程组时也方便得多。下面推导对A矩阵进行LU分解的方法。在以上推导过程中得到下式：=将其按矩阵乘法展开，就可以推导出L矩阵的第1列元素 U矩阵的第1行元素 L矩阵的第2列元素 U矩阵的第2行元素 若按此方式对列和行依次交换计算下去，最终即将求出 L和U矩阵的全部元素，其通用公式为：以上公式给出了对A矩阵进行LU分解的基本方法，在此基础上就可以编程实现LU分解法解线性方程组。程序设计按以上公式可以归纳出如下的求解线性代数方程组得LU分解算法：输入：A，非奇异矩阵nn。1. 令m=1。2. 计算L的第m列：3. 若m=n，则分解结束。4. 计算U的第m行：5. m=m+1，返回第2步。6. 前代： 7. 回代：LU分解法的计算工作量约为。由上述步骤可以得出如下的程序流程图：输入A矩阵开始计算L的第k列令k=1计算U的第k行k=k+1k=n？回代求Xk前代求Yk结束NY程序源代码/*/-基于C语言实现大规模节点方程的LU分解算法-/- 0403205班 040320510 邵汉钦 -/- 完成日期:2006年3月25日 -/*/程序功能描述：/设n元线性方程组的矩阵形式为 AX=B/输入矩阵A和B,根据LU分解法求解方程组，得到矩阵X/*/-#include #pragma hdrstop#include Unit1.h#include #include #include #include #include #include /-#pragma package(smart_init)#pragma resource *.dfm#define N 230TForm1 *Form1;float AN+1N+1,XN+1,YN+1,BN+1,s,s1,s2,det;/det为矩阵的行列式int i,j,k,p,num;/-_fastcall TForm1:TForm1(TComponent* Owner) : TForm(Owner)/-void _fastcall TForm1:Button1Click(TObject *Sender) int button; String edtOpenFileName;fstream Afile; num=StrToInt(InputBox(方程阶数,输入方程阶数n:,); OpenDialog1-Filter=Text Files(*.txt)|*.txt; OpenDialog1-InitialDir=Application-ExeName; OpenDialog1-Title=打开文件; OpenDialog1-Options.Clear(); OpenDialog1-OptionsofFileMustExist ofHideReadOnly ofNoChangeDir Execute() edtOpenFileName=OpenDialog1-FileName; TMsgDlgButtons buttons; buttonsmbOK; /-从文件中读取A矩阵- Afile.open(edtOpenFileName.c_str(),ios:in); if(!Afile)cout文件打开失败!n;exit(0);for(i=1;i=num;i+) for(j=1;jAij; /-从文件中读取B矩阵- char ch; Afilech; if(ch!=#)button=Application-MessageBoxA(输入数据有错误!n请检查存放数据的文本文档是否有错误!,错误提示,64); if(button=1)Close();for(i=1;iBi; button=Application-MessageBoxA(数据已成功输入!,数据输入,64);/-void _fastcall TForm1:Button2Click(TObject *Sender) int button=0; float temp; int re=1; /-选主元- temp=A11; for(i=2;itemp) temp=Ai1; re=i; if(re!=1) for(i=1;i=num;i+) temp=Arei; Arei=A1i; A1i=temp; temp=Bre; Bre=B1; B1=temp; /-计算L矩阵和U矩阵-for(k=1;k=num;k+)for(i=k;i=num;i+)s=0;for(p=1;p=k-1;p+)s=Aip*Apk+s;Aik=Aik-s;for(j=k+1;j=num;j+)s=0;for(p=1;p=k-1;p+)s=Akp*Apj+s;Akj=(Akj-s)/Akk;/-计算A矩阵的行列式值-s1=1;s2=1;for(i=1;iImage1-Canvas-TextOut(0,0,矩阵A的行列式值为:); Form1-Image1-Canvas-TextOut(10,20,det);/-矩阵A的行列式值为0时方程组没有唯一解的情况-if(det=0)button=Application-MessageBoxA(矩阵A的行列式值为0n方程组没有唯一解!,无解情况,48); if(button=1)Close(); /-计算并输出Y矩阵-Y1=B1/A11;for(i=2;i=num;i+)s=0;for(j=1;j=1;i-)s=0;for(j=i+1;jImage1-Canvas-TextOut(0,50,最后结果-X矩阵为:); for(i=1;iImage1-Canvas-TextOut(10,40+30*i,Xi);/-void _fastcall TForm1:Button3Click(TObject *Sender) int button; button=Application-MessageBoxA(确认结束吗？,结束,0); if(button=1)Close();/-void _fastcall TForm1:FormCreate(TObject *Sender)Label1-Caption=LU分解法求解n元线性方程组nn设方程组的矩阵形式为AX=Bn输入矩阵A和B，得到矩阵X;/-程序运行结果程序的调试和运行都在C+ Builder中进行，采用了可视化的程序界面。程序的运行情况演示如下：双击可执行文件“LU.exe”，出现如下的程序界面。主界面共有三个按钮，分别为“读取A和B矩阵”、“求解最后结果”和“结束”按钮，依次作为程序的输入、输出和关闭程序的功能使用。单击“读取A和B矩阵”按钮，按照程序的提示分别输入方程的未知数个数（即A矩阵的阶数），并从打开文件对话框中选择存放A和B矩阵的.txt文件。若打开文件正确，则会出现“数据已成功输入”的提示框，单击“确定”按钮继续。以下是各步的程序界面。在输入框中输入方程的阶数： 从打开对话框中选择要打开的存放A和B矩阵的文件： 程序提示“数据已成功输入”，单击“确定”按钮继续程序的运行。注意：若输入的文本文档中有错误或者输入的阶数与文件中的不符，程序则会提醒出错，如下图所示。单击“确定”按钮关闭程序，重新启动输入。 程序回到主界面，单击“求解最后结果”按钮，程序将给出矩阵A的行列式值和未知数X向量的值。如下图所示：最后，若要关闭程序，单击“结束”按钮，出现“确认结束”对话框，单击“确定”关闭程序，如下图所示：注意：txt文件中A、B矩阵的存放格式必须如下图所示。前面按矩阵的形式存放A，后面以列向量的形式存放B，中间必须以“#”隔开，否则会出错。输入不同的方程组时，只需在txt文件中改变相应的A、B矩阵即可。算法及程序的改进1.方程组无解情况的讨论在线性方程组求解的过程中，首先应该判断方程组在什么情况下无解。若能判断出方程组无解，这样就可以直接得到结论，而避免进行一系列的复杂运算而浪费时间。因此，在求解之前对方程组进行有无解得初步判断非常重要。对于n元一次线性方程组，若系数矩阵A的行列式值为零，则方程组没有唯一解；若A的行列式不为零，则不能肯定方程组是否有解，需要进一步的判断。从以上讨论中可以看出，在求解方程组之前要首先对系数矩阵A的行列式值进行计算和判断，需要在程序中添加一段求解行列式值的源代码。而根据行列式值计算的公式，求行列式涉及到大量计算，这对算法的时间和空间复杂度都不利。能不能有更好的方法呢？从A矩阵的LU分解中我们得到了启发。根据LU分解的原则，得到下式：=上式中三个矩阵分别称为A、L、U矩阵，它们都是nn的矩阵。根据线性代数中的理论，有以下定理：若A与B都是n阶矩阵，则由上式中A=LU的关系，可得。而L为下三角矩阵，U为上三角矩阵，它们的行列式值计算都很简单，即为其对角线元素的乘积。这样就可以利用LU分解的结果来简化行列式值的计算，从而对方程组有无解的情况作出初步判断。在本程序中就使用了这种方法计算A的行列式值。2.选主元在LU分解法的过程中，都是假定主元素（主元）不为零，从而使计算过程得以继续进行下去。但在实际上会出现主元为零或者值很小的情况，对于这种病态方程，出现主元为零会使计算无法进行下去；主元值很小会由于除法运算时的舍入误差而严重影响计算精度或者数值溢出，从而可能导致算法的失败。为避免上述情况的出现，保证消元法过程中的数值稳定性，需要进行选择主元的工作。事实上，只需交换方程组中有关方程的顺序，也就是所谓的选主元，就能使问题得到解决。对于LU分解法，可采取以下的选主元方法：在系数矩阵A中选取本列绝对值最大的元素作为主元素，这可以通过交换两行的元素使所选的元素处于主元的位置，然后在B矩阵中交换同样位置的两行元素。这就相当于在方程组中交换两个方程的顺序，不会使未知数向量中未知数的顺序发生变化，对最后的解的结果没有影响。在本程序中采用了以上选主元的方法。 3.稀疏矩阵技术在分析大规模电路时，代求的方程组系统中的未知数可能会数以千计。这样，方程组得系数矩阵就会十分庞大。另外，电路方程一般都是稀疏的，含有许多零元素，这是因为在通常的电路中，一个节点不可能与其他所有节点都有电学连接，而是只与少数几个连接，因此，关联矩阵中大量的交叉项为零。基于上述原因，有必要在求解线性方程组时，引入稀疏矩阵技术来节省存储量和计算时间。目前，比较常用的存储稀疏矩阵的数据结构是双链表结构（十字链表结构）。这种数据结构及其有关的一系列算法比较复杂，程序实现比较困难，因此在本程序中未采用。本程序中采用的是一种如下的简单的数据压缩存储方法。在方