电路方程的矩阵形式 c语言,基于C语言实现大规模节点方程的LU分解算法.doc

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

第 1 页 共 18 页基于 C 语言实现 大规模节点方程的 LU 分解算法班级 0403205 班学号 040320510姓名 邵汉钦完成日期 2006-03-25第 2 页 共 18 页基于 C 语言实现大规模节点方程的 LU 分解算法实验目的及主要内容通过编程实现基于 C 语言的大规模节点方程的 LU 分解算法，熟悉在电路仿真中运用的节点分析法计算机实现过程，对算法的收敛性、效率、计算过程中节省内存的手段和重要性有清晰的认识，从而能更好地运用商用大规模仿真软件或自己动手编写仿真软件。实验原理对复杂网络所建立的节点方程组，常常是高阶线性方程组，所以说，对线性电路的稳态分析，最终可归结为高阶线性方程组得求解问题。线性方程组的解法很多，一般可以分为两大类直接法和迭代法。直接法是通过对方程组直接进行一系列的运算而得到方程的解；迭代法是一开始就给出一个近似的解，然后通过逐次逼近来求 LAE 较好的近似解，迭代过程要一直进行到近似解达到预定的精度为止(即收敛到精确解) 。直接法的两个基本方法是高斯消元法和 LU 分解法，它们也是稀疏矩阵技术的基础。LU 分解是矩阵三角化分解法之一，也称三角分解法。它是把方程中的系数矩阵分解成下三角矩阵 L 和上三角矩阵 U 的乘积ALU式中，L 为主对角线以上的元素均为零的下三角矩阵；U 为主对角线以下的元素均为零，且主对角线元素均为 1 的上三角矩阵。上式展开就是第 3 页 共 18 页nnnaa 212112 nnlll 2110123112 nu则方程变为 AXLUXB一旦 A 分解成功，我们就可以设 UXY，方程就变成一个下三角形系统LYB，解此方程可以用所谓“前代”的方法，即从第 1 个方程求出 ，进1y而从第 2 个方程求出 ，依次类推，直到求出 ，把 ， 代2yny2n入第 n 个方程求出 。然后求解上三角系统 UXY，这样就可以用回代的方n法求出 X 来。LU 分解法的优点主要表现在当系数矩阵 A 不变，即只是网络元件参数 不变，仅仅是右手边向量 B 变化时，即外加激励信号变化时，可以减少解方程的计算工作量。因为一旦系数矩阵分解成 L 和 U 的乘积，就可以把L，U 存储起来，每次右手边向量变化时，只需做前代和回代，即可完成方程的求解。它在灵敏度计算中需要求转置的方程组时也方便得多。下面推导对 A 矩阵进行 LU 分解的方法。在以上推导过程中得到下式nnnaa 212112 nnlll 2110123112 nu将其按矩阵乘法展开，就可以推导出L 矩阵的第 1 列元素 ilii ,1U 矩阵的第 1 行元素 njlaujj 32L 矩阵的第 2 列元素 iliii ,12U 矩阵的第 2 行元素 jlujjj ,42第 4 页 共 18 页若按此方式对列和行依次交换计算下去，最终即将求出 L 和 U 矩阵的全部元素，其通用公式为 nmiulalmkiii ,11illukimimi ,2,1 以上公式给出了对 A 矩阵进行 LU 分解的基本方法，在此基础上就可以编程实现 LU 分解法解线性方程组。程序设计按以上公式可以归纳出如下的求解线性代数方程组得 LU 分解算法输入A，非奇异矩阵 nn。1. 令 m1。2. 计算 L 的第 m 列 nmiulalkmiii ,113. 若 mn，则分解结束。4. 计算 U 的第 m 行 nmilulaumkiii ,2,11 5. mm1，返回第 2 步。6. 前代 nmilulbymmki ,2,11 7. 回代第 5 页 共 18 页1,1niuxyxnikikmLU 分解法的计算工作量约为 。3由上述步骤可以得出如下的程序流程图第 6 页 共 18 页输入A 矩阵开始计算 L 的第 k 列令 k1计算 U 的第 k 行kk1kn回代求 Xk前代求 Yk结束NY程序源代码第 7 页 共 18 页基于 C 语言实现大规模节点方程的 LU 分解算法- 0403205 班 040320510 邵汉钦 - 完成日期2006 年 3 月 25 日 -程序功能描述设 n 元线性方程组的矩阵形式为 AXB输入矩阵 A 和 B,根据 LU 分解法求解方程组，得到矩阵 X-include pragma hdrstopinclude “Unit1.h“include include include include include include -pragma packagesmart_initpragma resource “*.dfm“define N 230T1 *1;float AN1N1,XN1,YN1,BN1,s,s1,s2,det;det 为矩阵的行列式int i,j,k,p,num;-fastcall T1T1TComponent* Owner TOwner-void fastcall T1Button1ClickTObject *Senderint button;String edtOpenFileName;fstream Afile;numStrToIntBox“方程阶数“,“输入方程阶数 n“,;OpenDialog1-Filter“Text Files*.txt|*.txt“;OpenDialog1-InitialDirApplication-Name;OpenDialog1-Title“打开文件“;第 8 页 共 18 页OpenDialog1-Options.Clear;OpenDialog1-OptionscuteedtOpenFileNameOpenDialog1-FileName;TMsgDlgButtons buttons;buttonsmbOK;-从文件中读取 A 矩阵-Afile.openedtOpenFileName.c_str,iosin;ifAfilecoutAij;-从文件中读取 B 矩阵-char ch;Afilech;ifchbuttonApplication-MessageBoxA“输入数据有错误n 请检查存放数据的文本文档是否有错误“,“错误提示“,64;ifbutton1Close;fori1;iBi;buttonApplication-MessageBoxA“数据已成功输入“,“数据输入“,64;-void fastcall T1Button2ClickTObject *Senderint button0;float temp;int re1;-选主元-tempA11;fori2;itemptempAi1;rei;ifre1fori1;iImage1-Canvas-TextOut0,0,“矩阵 A 的行列式值为“;1-Image1-Canvas-TextOut10,20,det;-矩阵 A 的行列式值为 0 时方程组没有唯一解的情况-ifdet0第 10 页 共 18 页buttonApplication-MessageBoxA“矩阵 A 的行列式值为 0n 方程组没有唯一解“,“无解情况“,48;ifbutton1Close;-计算并输出 Y 矩阵-Y1B1/A11;fori2;i1;is0;forji1;jImage1-Canvas-TextOut0,50,“最后结果X 矩阵为“;fori1;iImage1-Canvas-TextOut10,4030*i,Xi;-void fastcall T1Button3ClickTObject *Senderint button;buttonApplication-MessageBoxA“确认结束吗“,“结束“,0;ifbutton1Close;-void fastcall T1CreateTObject *SenderLabel1-Caption“LU 分解法求解 n 元线性方程组nn 设方程组的矩阵形式为 AXBn 输入矩阵A 和 B，得到矩阵 X“;-程序运行结果程序的调试和运行都在 C Builder 中进行，采用了可视化的程序界面。第 11 页 共 18 页程序的运行情况演示如下双击可执行文件“LU.” ，出现如下的程序界面。主界面共有三个按钮，分别为“读取 A 和 B 矩阵” 、 “求解最后结果”和“结束”按钮，依次作为程序的输入、输出和关闭程序的功能使用。单击“读取 A 和 B 矩阵”按钮，按照程序的提示分别输入方程的未知数个数(即 A 矩阵的阶数) ，并从打开文件对话框中选择存放 A 和 B 矩阵的.txt 文件。若打开文件正确，则会出现 “数据已成功输入”的提示框，单击“确定”按钮继续。以下是各步的程序界面。在输入框中输入方程的阶数第 12 页 共 18 页从打开对话框中选择要打开的存放 A 和 B 矩阵的文件程序提示“数据已成功输入” ，单击“确定”按钮继续程序的运行。第 13 页 共 18 页注意若输入的文本文档中有错误或者输入的阶数与文件中的不符，程序则会提醒出错，如下图所示。单击“确定”按钮关闭程序，重新启动输入。 程序回到主界面，单击“求解最后结果”按钮，程序将给出矩阵 A 的第 14 页 共 18 页行列式值和未知数 X 向量的值。如下图所示最后，若要关闭程序，单击“结束”按钮，出现“确认结束”对话框，单击“确定”关闭程序，如下图所示注意txt 文件中 A、B 矩阵的存放格式必须如下图所示。前面按矩阵的形式第 15 页 共 18 页存放 A，后面以列向量的形式存放 B，中间必须以“”隔开，否则会出错。输入不同的方程组时，只需在 txt 文件中改变相应的 A、B 矩阵即可。算法及程序的改进1.方程组无解情况的讨论在线性方程组求解的过程中，首先应该判断方程组在什么情况下无解。若能判断出方程组无解，这样就可以直接得到结论，而避免进行一系列的复杂运算而浪费时间。因此，在求解之前对方程组进行有无解得初步判断非常重要。对于 n 元一次线性方程组，若系数矩阵 A 的行列式值为零，则方程组没有唯一解；若 A 的行列式不为零，则不能肯定方程组是否有解，需要进一步的判断。从以上讨论中可以看出，在求解方程组之前要首先对系数矩阵 A 的行第 16 页 共 18 页列式值进行计算和判断，需要在程序中添加一段求解行列式值的源代码。而根据行列式值计算的公式，求行列式涉及到大量计算，这对算法的时间和空间复杂度都不利。能不能有更好的方法呢从 A 矩阵的 LU 分解中我们得到了启发。根据 LU 分解的原则，得到下式nnnaa 212112 nnlll 2110123112 nu上式中三个矩阵分别称为 A、L、U 矩阵，它们都是 nn 的矩阵。根据线性代数中的理论，有以下定理若 A 与 B 都是 n 阶矩阵，则 BA由上式中 ALU 的关系，可得 。而 L 为下三角矩阵，U 为上三L角矩阵，它们的行列式值计算都很简单，即为其对角线元素的乘积。这样就可以利用 LU 分解的结果来简化行列式值的计算，从而对方程组有无解的情况作出初步判断。在本程序中就使用了这种方法计算 A 的行列式值。2.选主元在 LU 分解法的过程中，都是假定主元素(主元)不为零，从而使计算过程得以继续进行下去。但在实际上会出现主元为零或者值很小的情况，对于这种病态方程，出现主元为零会使计算无法进行下去；主元值很小会由于除法运算时的舍入误差而严重影响计算精度或者数值溢出，从而可能导致算法的失败。为避免上述情况的出现，保证消元法过程中的数值稳定性，需要进行第 17 页 共 18 页选择主元的工作。事实上，只需交换方程组中有关方程的顺序，也就是所谓的选主元，就能使问题得到解决。对于 LU 分解法，可采取以下的选主元方法在系数矩阵 A 中选取本列绝对值最大的元素作为主元素，这可以通过交换两行的元素使所选的元素处于主元的位置，然后在 B 矩阵中交换同样位置的两行元素。这就相当于在方程组中交换两个方程的顺序，不会使未知数向量中未知数的顺序发生变化，对最后的解的结果没有影响。在本程序中采用了以上选主元的方法。 3.稀疏矩阵技术在分析大规模电路时，代求的方程组系统中的未知数可能会数以千计。这样，方程组得系数矩阵就会十分庞大。另外，电路方程一般都是稀疏的，含有许多零元素，这是因为在通常的电路中，一个节点不可能与其他所有节点都有电学连接，而是只与少数几个连接，因此，关联矩阵中大量的交叉项为零。基于上述原因，有必要在求解线性方程组时，引入稀疏矩阵技术来节省存储量和计算时间。目前，比较常用的存储稀疏矩阵的数据结构是双链表结构(十字链表结构) 。这种数据结构及其有关的一系列算法比较复杂，程序实现比较困难，因此在本程序中未采用。本程序中采用的是一种如下的简单的数据压缩存储方法。在方程组得求解过程中，需要用到三个 nn 的矩阵，即 A、L 和 U 矩阵，因此在程序中就分别要为它们分配 个单元，这样做比较浪费存储空间。由 A 矩阵和2nL、U 矩阵的关系，想到节省内存的方法。由于将 A 矩阵分解成下三角矩阵L 和上三角矩阵 U，而 L、U 矩阵都是只有一半元素非零，一般元素为零，并且是互补的，因此可以考虑将 LU 分解的结果直接存放在 A 矩阵中，而不第 18 页 共 18 页必在内存中为 L、U 矩阵分配空间，这样就能节省三分之二的内存空间。事实证明，这种方法是可行的，对计算结果无影响，只是在需要输出 L 和 U矩阵时将它们的那些零元素补充出来即可。本程序即采用了这种方法。