利用直接发求解线性方程组及其应用

1、 利用直接法解线性方程组及其应用摘要: 在科学研究和工程技术中有许多问题可归结为求解线性代数方程组的问题和常微分方程的定解问题。本文主要讨论了解线性方程组的各种直接解法，常微分方程的边值问题及其应用，最后介绍了病态方程和及其求解方法。关键字: 线性方程组 直接法 应用举例一、引言随着科学技术的发展，提出了大量复杂的数值计算问题，在建立电子计算机成为数值计算的主要工具以后，它以数字计算机求解数学问题的理论和方法为研究对象。有可靠的理论分析，要有数值实验，并对计算的结果进行误差分析。数值分析的主要内容包括插值法，函数逼近，曲线拟和，数值积分，数值微分，解线性方程组的直接方法，解线性方程组的迭代法，

2、非线性方程求根，常微分方程的数值解法。运用数值分析解决问题的过程包括：实际问题数学建模数值计算方法程序设计上机计算求出结果。在自然科学研究和工程技术中有许多问题可归结为求解线性代数方程组的问题,线性方程组求解是科学计算中最常遇到的问题。如在应力分析、电路分析、分子结构、测量学中都会遇到解线性方程组问题。在很多广泛应用的数学问题的数值方法中，如三次样条、最小二乘法、微分方程边值问题的差分法与有限元法也都涉及到求解线性方程组。直接法是在没有舍入误差的情况下，通过有限步四则运算来求的方程组精确解的方法。直接法包括高斯列主元消去法，三角分解法，追赶法等的基本思想和原理，学习它们各自的优缺点及适用范围。

3、能够利用直接法解线性方程组解决一些如道路的交通流量，电网的电流流量，解析几何等的实际问题。二、解线性方程组的方法关于线性代数方程组的数值解法一般分为两大类：直接法和迭代法。1、 直接法 直接法就是经过有限步算术运算，可求的线性方程组精确解的方法（若计算过程没有舍入误差），但实际犹如舍入误差的存在和影响，这种方法也只能求得近似解，这类方法是解低阶稠密矩阵方程组级某些大型稀疏矩阵方程组的有效方法。直接法包括高斯消元法，矩阵三角分解法、追赶法、平方根法。2、 迭代法迭代法就是利用某种极限过程去逐步逼近线性方程组精确解的方法。将方程组的解看作是某极限过程的极限值，且计算这一极限值的每一步是利用前一步所

4、得结果施行相同的演算步骤而进行。迭代法具有需要计算机的存储单元少，程序设计简单，原始系数矩阵在计算过程始终不变等优点，但存在收敛性级收敛速度问题。迭代法是解大型稀疏矩阵方程组（尤其是微分方程离散后得到的大型方程组）的重要方法。迭代法包括Jacobi法SOR法、SSOR法等多种方法。本文主要介绍用高斯消元法，矩阵三角分解法求解线性方程组及其应用。三、具体解法及其实例1、Gauss消元法通过一系列的加减消元运算，也就是代数中的加减消去法，以使A对角线以下的元素化为零，将方程组化为上三角矩阵；然后，再逐一回代求解出x向量。现举例说明如下：1.1消元过程第一步：将(1)/3使x1的系数化为1,再将(2

5、)、(3)式中x1的系数都化为零，即由(2)-2×(1)(1)得 由(3)-4×(1)(1)得 第二步：将(2)(1)除以2/3，使x2系数化为1，得再将(3)(1)式中x2系数化为零，由(3)(1)-(-14/3)*(2)(2) ，得第三步：将(3)(2)除以18/3，使x3系数化为1，得经消元后，得到如下三角代数方程组：1.2 回代过程由(3)(3)得 x3=1,将x3代入(2)(2)得x2=-2,将x2 、x3代入(1)(1)得x2=1,所以，本题解为x=1,2,-1T 1.3 用矩阵演示进行消元过程第一步： 先将方程写成增广矩阵的形式第二步：然后对矩阵进行初等行变换

6、第三步：将增广矩阵变换成上三角矩阵，即主对角线全为1，左下三角矩阵全为0，形式如下：即原方程组被等价转化成为上三角方程组，然后，逐步回代得原方程组的解即可。1.4高斯消元的公式综合以上讨论，不难看出，高斯消元法解方程组的公式为第一步，消元（1） 令aij(1) = aij , （i,j=1,2,3,n）bi(1) =bi ,（i=1,2,3,n）（2） 对k=1到n-1，若akk(k)0，进行lik = aik (k) / akk(k) ,（i=k+1,k+2,n）aij(k+1) = aij(k) - lik * akj(k), （i,j= k+1,k+2,n）bi(k+1) = bi(k)

7、 - lik * bk(k), （i= k+1,k+2,n）第二步，回代若ann(n) 0xn = bn(n) / ann(n)xi = (bi(i) sgm(aij(i) * xj )/- aii(i) ,（i = n-1,n-2,1）,( j = i+1,i+2,n )1.5 解的判断 设方程组的增广矩阵记为，则经过初等行变换可化为如下的阶梯形矩阵(必要是可重新排列未知量的顺序)：其中cii¹0(i=1，2，r)于是可知：(1). 当dr+1=0，且r=n时，原方程组有唯一解(2). 当dr+1=0，且r<n时，原方程组有无穷多解(3). 当dr+1¹0，原方程组

8、无解2、LU分解法求解线性代数方程组除了高斯消元法外，还常用LU分解法（三角形分解法）。LU分解法的优点是当方程组左端系数矩阵不变，仅仅是方程组右端列向量改变，即外加激励信号变化时，能够方便地求解方程组。矩阵的三角分解法可分为直接三角分解法,列主元三角分解法,平方根法,三对角方程组的追赶法。下面讨论直接三角分解法。设n阶线性方程组Ax=b 。假设能将方程组左端系数矩阵A，分解成两个三角阵的乘积，即A=LU ，式中，L为主对角线以上的元素均为零的下三角矩阵, 且主对角线元素均为1的上三角矩阵；U为主对角线以下的元素均为零所以有，LUx=b 令 Ux=y， 则 Ly=b由A=LU,由矩阵的乘法公式

9、： a1j = u1j , j=1,2,nai1 = li1u11 , i=1,2,n 推出 u1j = a1j, j=1,2,nli1 = ai1/u11, i=1,2,n 这样就定出了U的第一行元素和L的第一列元素。设已定出了U的前k-1行和L的前k-1列，现在确定U的第k行和L的第k列。由矩阵乘法：当r>k时，lkr=0, 且lkk=1，因为所以，同理可推出计算L的第k列的公式：因此得到如下算法杜利特（Doolittle）算法:(1) 将矩阵分解为A=LU,对k=1,2,n；j=k,k+1,n； i=k,k+1,n；公式1 (2) 解Ly=b (3) 解Ux=y对大规模稀疏问题，如

10、果能够通过调整方程及未知量的顺序使得方程组的系数矩阵成带状结构，则对系数矩阵使用通常的LU分解，可以保障单位下三角矩阵L及上三角矩阵U仍为带状结构.3. 直接法解线性方程组的应用例1：确定下图电网中的回路电流。解 在回路1中,电流I1流过三个电阻,且电压降为；在回路2中的电流也流经回路1的一部分,即从D到A的分支,对应的电压降为6I2伏特然而,回路1中电流在DA段的方向与回路2中选定的方向相反,因此,回路1中所有电压降的代数和为由于回路1中的电压为+60伏特,由基尔霍夫电压定律,可得回路1的方程为,同理,可得回路2的方程为 , 其中, -6I1是回路1中流经DA分支的电流(因为电流与回路2中的

11、电流方向相反,所以电压为负)；12I2是回路2中所有的电阻乘上回路电流的和； -2I3是回路3中流经CB分支上2欧姆电阻的电流,方向与回路2中该段的电流方向相反。回路3的方程为 注意,在CB分支上10伏特的电池被当作是回路2和回路3中的一部分,但是由于回路3中电流方向,电池在回路3中为-10伏特出于同样的道理,40伏特的电池也应取负值.综合上述讨论,上述电网的回路电流满足下列线性方程组 写成矩阵形式为 (*)对增广矩阵进行行变换,得从而解得I1=3安培, I2=1安培, I3=-8安培I3取负值说明回路3中的实际电流与图中显示的电流方向相反.在方程组(*)中,如果将其系数矩阵记为R,右端列向量

12、记为u,i =(I1,I2,I3)T,则可得到以矩阵形式表示的欧姆定律:.4 由于计算机的发展，人们开始考虑在计算机上求解线性方程组Ax=b的近似解，并用某种极限过程去逐渐逼近精确解。如初等函数及其他特殊函数，当函数只在有限的点集上给定函数值，要在包含该点集的区间上用公式写出函数简单表达式，这些都涉及在去见a,b上用简单的函数逼近一直复杂函数的问题，这也就是函数逼近问题。下面通过实例说明问题。在提拉单晶硅的过程中，假定晶体导电，熔化物和晶体的平衡方程式是：该模型还包括坩埚的传导,熔化物、晶体、坩埚、加热元件表面和环境的热辐射,以及融化物和周围的环境的对流。利用格林第二定理，将面积分变成线积分X

13、i是边界节点i出的方向矢量。是边界线，A是横截面积，对轴对称的几何解决方法，i扩散方程式：Km是椭圆部分。（ri,zi）是节点i的坐标。在数值实现方法、边界被分为温度变化和接近温度梯度的每个值，通过插值逼近中间节点值:是k的函数，是坐标，Ne是每个元素的节点数。方程(10) 、(12)代入(9)，形成N每个节点的微分/代数方程：T是节点处的温度矢量，Q是节点的变化矢量。为了获得一个边界积分相当于瞬时的域积分,一般利用插值(23)。在这种方法中,对时间的导数形式是:X是位置矢量，x= r z T,N是节点数，fkx; xk是几何的已知函数，kt是时间的未知系数，方程（13）变为：矩阵M的整数数域

14、，注意到这矩阵的计算在一个固定的领域只计算一次,方程(16)应用到N结点,结果是:矩阵F是通过fkx; xk在节点处形成的，方程（18）代入方程（17），整数域里方程（9）变成：另外,为了避免整数域M在较复杂的几何域,利用双重定理，在这种方法中,节点fk x;Xk 是取代拉普拉斯节点的一个新的节点, ：再将方程（21）代入（17），得到一个积分方程：如果节点，估测有相同的插值多项式k，利用温度，矩阵M可写成：, （23）,矩阵 ,是通过节点形成的。该模型的非线性介绍是通过自由表面传导结合热对流及热辐射效应:例2 计算模型本文运用有限元的方法对生长的硅单晶和晶体生长的热场作为整体进行研究,得到了

15、不同条件下单晶炉内的热场分布。对晶体实际生产有一定的指导作用。 根据热动力学第一定律即热能守恒定律和热扩散定律,可以得到(1)式:将(1)式进行展开得到其常见形式:(2)式是有限元计算中用到的偏微分方程。式中为物质密度,C为热容量,T为温度(T=T(x, y,z, t), t为时间, L是矢量运算符, V是热量物质输运矢量, q为热流矢量,q为单位体积内热量产生率, D=是热量传导矩阵,Kxx, Kyy, Kzz是热量在各坐标轴方向的导热系数。在计算中采取了近似方法5,假定: (1)固-液界面处的温度为硅的熔点,且界面形状为平面;(2)物性参数如热传导系数、密度、热容量等不随温度变化,都为常数

16、; (3)单晶炉内热场主要以热传导6和辐射7,8为主,不考虑熔体的热对流问题。计算过程中用到的材料属性和处理参数如表1所示。4、方程的性态与误差分析实际上由于误差的存在，直接法只能求得近似解。误差分析是计算方法中一个既重要复杂的问题。 因为几乎每一步运算都有误差，而实际问题往往需要进行千百万次计算，所以每步运算即分析误差几乎是不可能的，也是不必要的。设Axb（矩阵A非奇异）中，系数矩阵有小扰动和b，扰动后的方程组为A(A+A)(x +x) = b +b若A为n阶非奇异矩阵，称Cond(A) = A1 A为矩阵A的条件数。当条件数Cond(A)比较大时，系数矩阵和右端向量的小扰动会引起解的相对误差较大。如果矩阵A的条件数比较大，就说方程组是“病态”的；如果矩阵A的条件数比较小，就说方程组是“良态”的；方程组的性态是方程组的固有性质，它与求解方法无关。在计算机上求解的方程组都是所给方程组的扰动方程（字长限制将产生舍入误差）。对于良态方程组，只要求解方法稳定，即可得到比较满意的计算结果。但对于病态方程组，即使使用稳定性好的算法求解也未必理想。求解病态方