2022电力系统分析中的计算方法.docx

《2022电力系统分析中的计算方法.docx》由会员分享，可在线阅读，更多相关《2022电力系统分析中的计算方法.docx（249页珍藏版）》请在淘文阁 - 分享文档赚钱的网站上搜索。

1、电力系统分析中的计算方法（ 原书第 2 版）目录第 1 章引论 1第 2 章线性方程组的解法 32. 1 Gauss 消 去法42. 2 LU 分 解法82. 2. 1 采用部分选主元的 LU 分解法152. 2. 2 采用完全选主元的 LU 分解法192. 3 条件数与误差传播202. 4 松 弛法212. 5 共轭梯度法262. 6 广义最小残差法302. 7 问题35第 3 章非线性方程组的解法393. 1 不动点迭代法393. 2 Newton - Raphson 迭 代法453. 2. 1 收 敛性473. 2. 2 用于求解非线性方程组的 Newton - Raphson 法483

2、. 2. 3 Newton - Raphson 法的改进513. 3 连 续 法3. 4 割 线 法 52 553. 5 数值微分法573. 6 在电力系统中的应用603. 6. 1 潮 流 计算613. 6. 2 调节变压器683. 6. 3 解耦潮流算法723. 6. 4 快速解耦潮流算法743. 6. 5 PV 曲线与连续潮流计算77目录3. 6. 6 三相潮流计算823. 7 问题83第 4 章稀疏矩阵求解技术864. 1 存 储 方法864. 2 稀疏矩阵的表示方法904. 3 排序方案4. 3. 1 方 案 0 92 984. 3. 2 方 案994. 3. 3 方 案1044.

3、3. 4 其 他 方案1064. 4 在电力系统中的应用1074. 5 问题110第 5 章数值积分1145. 1 单 步法1155. 1. 1 基于 Taylor 级数的算法1155. 1. 2 向 前 Euler 法1155. 1. 3 Runge - Kutta 法1165. 2 多 步法1165. 2. 1 Adams 算法1225. 2. 2 Gear 法1255. 3 精度与误差分析1265. 4 数值稳定性分析1295. 5 刚 性 系统1335. 6 步 长 选择1375. 7 微分代数方程1405. 8 在电力系统中的应用1415. 8. 1 暂态稳定性分析1425. 8.

4、2 中期稳定性分析1495. 9 问题152第 6 章优化方法1576. 1 最小二乘状态估计1576. 1. 1 加权最小二乘状态估计1616. 1. 2 坏数据的检测1636. 1. 3 非线性最小二乘状态估计1666. 2 线 性 规划167电力系统分析中的计算方法 （ 原书第 2 版）6. 2. 1 单 纯 形法1686. 2. 2 内 点法1716. 3 非线性规划1766. 3. 1 二 次 规划1766. 3. 2 最速下降法1786. 3. 3 序贯二次规划法1826. 4 在电力系统中的应用1846. 4. 1 最 优 潮流1846. 4. 2 状 态 估计1946. 5 问

5、题198第 7 章特征值问题2027. 1 幂法2027. 2 QR 法2057. 2. 1 移 位 QR 法2117. 2. 2 缩 减法2127. 3 Arnoldi 法2127. 4 奇异值分解2197. 5 模 态 辨识2227. 5. 1 Prony 法2237. 5. 2 矩 阵 束法2257. 5. 3 Levenberg - Marquardt 法2267. 5. 4 Hilbert 变换2287. 5. 5 举例2297. 6 在电力系统中的应用2357. 7 问题236第 1 章引论在今天放松管制的大环境下， 本国电网已不能按照当初设计时设想好的方式运行。 因此， 系统分析

6、和计算对于预测并不断地更新电网的运行状态就显得非常重要。 系统分析和计算的内容包括估算当前的潮流和电压 （ 潮流分析和状态估计）， 确定系统的稳定极限 （ 连续潮流、 用于暂态稳定计算的数值积分和特征值分析）， 以及费用最小化计算 （ 最优潮流） 等。 本书讲述各种计算方法的入门知识， 这些计算方法构成了电力系统及其他工程和科学领域很多分析计算的基础， 也是无数商业软件所使用算法的基础。 通过了解支撑这些算法的理论， 读者可以更好地使用那些商业软件并做出更明智的决策 （ 即在仿真软件中选择更好的数值积分方法和确定更合适的积分步长等）。由于电网的规模巨大， 电力系统分析计算依靠手工几乎是不可能的

7、， 计算机是唯一真正可行的手段。 电力工业界是计算机技术的最大用户之一， 也是在大型计算机出现后最早接受计算机分析的工业部门之一。 虽然电力系统分析的第一批算法是在 20 世纪 40 年代开发出来的， 但直到 60 年代计算机才在电力工业界得到广泛使用。 今天用于大系统仿真和分析的很多技术和算法起初是为电力系统应用而开发的。随着电力系统更多地运行在重载状态下， 计算机仿真将会在电力系统的控制和安全性评估方面发挥更大的作用。 将商业化软件包用于仿真重载的电力系统时， 往往不能运行或者给出错误的结果。 为了正确理解和解释商业化软件包的仿真结果， 必须了解其内在的数值计算方法。 例如， 是系统确实呈

8、现出了仿真得到的行为特征， 还是仅仅因为数值计算不精确而导致的伪结果？ 如何弥补商业化软件包存在的算法缺陷？ 是通过选择更好的仿真参数， 还是通过将问题描述成更适合于数值计算的形式？ 对诸如此类的问题， 受过数值计算训练的用户可以做出更好的判断。 本书将给出多种广泛使用的数值计算方法的背景知识， 这些计算方法构成了众多电力系统分析和设计商业化软件包的基础。本书是为研究生计算方法课程而写的教材， 课程长度是一个学期。 尽管本书中的大多数例子是基于电力系统应用的， 但数值计算的理论是按照一般性的方式阐述的， 因而可应用于大范围的工程系统。 虽然完全领会本书中某些内容的微妙之处需要具有电力系统工程的

9、一些知识， 但这些知识对理解算法本身并不是前提条件。 相对于内容广泛的计算方法领域， 本书的叙述和例子仅仅是一个入门的导引， 并不能包罗万象。 本书阐述的很多算法一直存在着众多的改进算法， 并且仍然是当前研究的目标； 由于本书的意图是打基础， 因此很多新的进展并没有明确地包含在叙述中， 而是作为参考文献列出， 以方便感兴趣的读者进一步学习。 为了让读者能够容易地将书中的算例复现出来， 本书特意挑选了简单而足够完整的问题作为算例。 很多 “ 真实世界” 的问题其规模和范围要大得多， 然而本书给出的处理问题的方法， 应该可以为读者提供足够的工具， 以应对其所遇到的困难。本书中的大部分算例是通过 M

10、atlab 编程实现的， 虽然这是本书作者所采用的平台， 但实际上， 任何计算机语言都可用来实现本书中的算法， 没有理由偏好某个特定的平台或语言。第 2 章线性方程组的解法在很多工程和科学领域， 希望根据一组物理学上的关系式， 能够用数学的方法来确定系统的状态。 这些物理学上的关系式可以是描述诸如电路拓扑、 质量、重量、 力等物理特性的方程， 例如， 电路系统中针对每个节点的由注入电流、 网络拓扑和支路阻抗构成的约束方程。 在很多情况下， 描述已知量 （ 即输入） 与未知状态 （ 即输出） 之间的关系式是线性的。 因此， 一般性地可以将线性系统描述为Ax = b（2. 1）式中， b 是一个

11、n 1 向量， 为已知量； x 是一个 n 1 的未知状态向量； 而 A 是将 x 与 b 联系起来的 n n 阶矩阵。 暂时我们假定矩阵 A 是可逆的， 即非奇异的， 这样， 每个向量 b 将对应一个唯一的向量 x。 也就是矩阵 A - 1 存在， 且x = A - 1 b（2. 2）是式 （2. 1） 的唯一解。对式 （2. 1） 进行求解的自然方法是直接计算矩阵 A 的逆， 并将它与向量 b相乘。 计算 A - 1 的方法之一是采用 Cramer 法则：（A - 1 （ i，j） = det1 A）（ Aij ） T对 于 i = 1，n， j = 1，n （2. 3） 式中， A -

12、1 （ i，j） 是 A - 1 的第 i 行、 第 j 列元素， 而 Aij 是矩阵 A 的元素 aij 的余子式。这种方法需要计算 （ n + 1） 个行列式， 从而需要 2 （ n + 1 ）！ 次乘法运算。 当 n 值较大时， 计算量爆炸式增长， 不可控制。 因此， 必须开发替代的算法。一般地说， 求解式 （2. 1） 有两种基本的方法：1） 直接法， 也称为消去法。 该方法通过有限次算术运算找到上述方程的精确解 （ 指在计算机的精度内）。 如果不计计算机的舍入误差， 采用直接法得到的解 x 将是完全精确的。2） 迭代法。 这种方法的每次计算都基于相同的计算过程， 试图产生一系列不断逼

13、近真解的近似解。 当所获得的近似解达到预先设定的精度指标或者已确定迭代过程不再对精度有改进时， 终止迭代。求解方法的选择通常取决于所研究系统的结构， 某些系统更适合于采用特定的方法进行求解。 一般地， 直接法更适合于满矩阵， 而迭代法更适合于稀疏矩阵。 但是， 正像对于大多数一般性结论一样， 上述经验法则存在很多的例外。第 2 章线性方程组的解法92. 1Gauss 消去法Gauss 消去法是一种不需要显式计算 A - 1 而对式 （2. 1） 进行求解的方法。 因为 x 是直接求出的， 因此 Gauss 消去法是一种求解线性方程组的直接法， 并且是一种常用的直接法。 Gauss 消去法的基本

14、思路是利用第 1 个方程消去其余方程中的第 1 个未知量， 依次重复这个过程就能消去第 2 个未知量、 第 3 个未知量 直到整个消去过程完成， 此时， 第 n 个未知量可以直接从输入向量 b 中求出。 然后， 对上述方程进行递归回代就能求出所有的未知量。Gauss 消去法是这样一个过程， 它通过一系列初等行变换运算， 将 n （ n + 1）阶的增广矩阵 A | b变换为 n （ n + 1） 阶的矩阵其中，I | bAx = bA - 1 Ax = A - 1 bIx = A - 1 b = b x = b这样， 如果存在一系列初等行变换运算可以将矩阵 A 变换成一个单位矩阵 I， 那么，

15、 对向量b 进行相同系列的初等行变换运算就能将向量b 变换为方程组的解x。初等行变换运算包括如下 3 种矩阵运算： 1） 交换矩阵中的任意两行。2） 用一个常数乘任意一行。3） 将任意行的线性组合加到另一行上。选择一系列初等行变换运算将矩阵 A 变换成一个上三角矩阵， 该上三角矩阵的对角元为 1， 而所有下三角元素为零。 这个过程被称为 “ 向前消去过程”。向前消去过程的每一步可以通过对矩阵 A 乘以一个初等矩阵 来获得， 这里的初等矩阵 指的是通过对单位矩阵进行初等行变换运算就能得到的矩阵。例 2. 1 寻找一系列初等矩阵将如下矩阵变换成上三角矩阵。1348 A =2123 4358 927

16、4 解 2. 1为了使上述矩阵上三角化， 所采用的行变换运算应能使此矩阵对角线以下每一列的元素为零。 这可以通过如下的行变换运算来实现： 将对角线以下的每一行用其本身减去一个常数乘对角线行来代替， 该常数的选择应使此列中对角线以下的元素为零。 因此， A 的第 2 行应该用 （ 行 2 - 2 行 1） 来代替， 而对应此初等行变换运算的初等矩阵是1000 1 = - 2100 0010 0001 和1348 A =10- 5- 6- 13 4358 9274 注意除了第 2 行其他的行没有变化， 而第 2 行在第 1 个对角元下的元素变为了零。 类似地， 完成消去第 1 列任务的另外 2 个

17、初等矩阵是2 = 1000 0100 - 4010 3 = 0001 1000 0100 0010 - 9001 并且有1348 A =3 2 10- 5- 6- 13 0- 9- 11- 24 0- 25- 29- 68 （2. 4）现在向前消去过程将针对第 2 列展开， 首先是将第 2 个对角元下面的元素全部消去， 然后将第 2 个对角元变换为 1。 具体过程是1000 0100 4 = 0- 91050001 1000 0100 5 = 0010 250- 501 100-100 00 6 = 50010 0001 并且有1348 01613 6 5 4 3 2 1 A = 55 13

18、（2. 5）001- 3 00- 5- 5 继续向前消去过程， 有1000 01007 = 从而有0010 0051 1000 =80100 00- 50 0001 1348 01613 8 7 6 5 4 3 2 1 A = 55 （2. 6）0013 000- 6 最后有1000 0100 9 = 0010 1000- 6 和1348 01613 9 8 7 6 5 4 3 2 1 A = 55 （2. 7）0013 0001 从而完成了向前消去过程， 将矩阵变换成了上三角矩阵。一旦将矩阵变换成了上三角矩阵， 式 （2. 1） 的解向量x的逐次代换 （ 也称为 “ 回代”） 而求得。例 2

19、. 2利用例 2. 1 的上三角矩阵， 求如下方程的解。就可以通过状态量1348 x1 1 2123x2 1 = 4358 x3 1 9274 x4 1 解 2. 2注意一系列下三角矩阵的乘积仍然是下三角矩阵。 因此， 矩阵乘积W = 9 8 7 6 5 4 3 2 1（2. 8）是一个下三角矩阵。 将 W 与矩阵 A 相乘将得到一个上三角矩阵，WA = U（2. 9）即 U 是一个通过向前消去过程而得到的上三角矩阵。 在式 （2. 1） 的两边同时左乘 W， 可以得到WAx = Wb（2. 10）Ux = Wb（2. 11）= b（2. 12）式中， Wb = b。根据例 2. 11000

20、21W = 5- 500 29- 50 1- 1 6且1461 1 56 1 1 6 b = W = 5 6 1 3 1 2 这样，1348 x1 1 52 01613 x 5 1 5 x = 6（2. 13）0013 3 x4 3 0001 通过观察， x4 = 3 / 2。 由第 3 行可得 2 x3 = 6 - 3x4（2. 14）将 x4 的值代入到式 （2. 14） 中， 得到 x3 = 3 / 2。 类似地，x2 = 1 - 6 x3 - 13x4（2. 15）555将 x3 和 x4 的值代入到式 （2. 15） 中， 得到 x2 = - 11 / 2。 以类似的方法求解 x1

21、，2x1 = 1 - 3x2 - 4x3 - 8x4（2. 16）从而得到= - 1x1 - 1 （2. 17）x2 = 1 - 11 x3 2 3 x4 3 逐次将已求得的 x 值代回到方程中进行求解的步骤， 在 Gauss 消去法的求解过程中被称为 “ 回代”。 因此， Gauss 消去法包括两个主要步骤： “ 向前消去” 过程和 “ 回代” 过程。 向前消去过程是将矩阵 A 变换成其三角因子， 而回代过程则根据输入向量 b 和 A 的矩阵因子求解未知向量 x。 Gauss 消去法同时也为 LU 因子分解过程提供了框架。2. 2LU 分解法Gauss 消去法的向前消去过程产生了一系列与矩阵

22、 A 相关的上三角矩阵和下三角矩阵， 如式 （2. 9） 所示。 其中， 矩阵 W 是一个下三角矩阵， 而矩阵 U 是一个对角元为 1 的上三角矩阵。 考虑到下三角矩阵的逆仍然是一个下三角矩阵， 因此如果定义L = W - 1那么A = LU正是基于矩阵 L 和矩阵 U， 将矩阵的因子分解或消去算法命名为矩阵的 “ LU 分解”。 事实上， 对于任何非奇异的矩阵 A， 存在一些置换矩阵 P （ 有可能 P = I）， 使得LU = PA（2. 18）式中， U 是对角元为 1 的上三角矩阵， L 是对角元为非零的下三角矩阵， 而 P 是由单位矩阵 I 通过行交换或列交换而形成的元素仅为 0 或

23、 1 的矩阵。 一旦选定一个合适的矩阵 P， 那么上述的 LU 分解就是唯一的6 。 如果 P、 L 和 U 已经确定， 那么方程组Ax = b（2. 19）就可以迅速得到求解。 在式 （2. 19） 两边同时左乘矩阵 P 得PAx = Pb = b（2. 20）LUx = b（2. 21）式中， b仅仅是对向量 b 的一个重新排列。 如果引入一个 “ 哑” 向量 y， 使得Ux = y（2. 22）那么考虑到式 （2. 23） 的结构为Ly = b（2. 23）l11000 y1 l21l2200 y2b1 b2 l31l32l330 y3 = b3 ln1ln2ln3lnn yn bn 向

24、量 y 的元素可以通过直接代换得到：y1 =b1 l11y2 =y3 =l1 （ b2 - l21 y1 ）2233l1 （ b3 - l31 y1 - l32 y2 ）lyn = 1nn（bnn-1-lnj yjj = 1求出 y 以后， x 就可以很容易根据如下方程求得：第 2 章线性方程组的解法39 1u12u13u1n x1 01u23u2n x2 y1 y2 n 001u3 x3 = y3 0001类似地， 解向量 x 可以通过回代得到：xn = ynxn-1 = yn-1 - un-1，n xn xn yn xn-2 = yn-2 - un-2，n xn - un-2，n-1 xn

25、-1nx1 = y1 - j = 2 u1j xjLU 分解的价值在于， 一旦矩阵 A 被分解成上三角矩阵和下三角矩阵， 那么寻找解向量x 就是一件直截了当的事了。 注意上述求解过程并没有显式地对矩阵A 求逆。存在多种方法对矩阵进行 LU 因子分解， 每种方法具有各自的优缺点。 一种常用的 LU 分解方法被称为 Croutb 算法6 。 定义矩阵 Q 为l11u12u13u1n l21l22u23u2n QL + U - I = l31l32l33u3n （2. 24） ln1ln2ln3lnn Crout 算法逐列逐行地计算 Q 的元素， 首先计算 Q 的列元素， 然后计算 Q 的行元素，

26、如图 2. 1 所示。 Q 中的每个元素 qij 只依赖于 A 的元素 aij 和前面已经计算得到的 Q 的元素。对矩阵 A 进行 LU 分解的 Crout 算法（1） 对矩阵 Q 进行初始化， 即清零， 并令 j = 1；（2） 根据下式计算 Q 的第 j 列 （ 即矩阵L 的第 j 列）：图 2. 1 矩阵 Q 列和行元素的计算次序qkj= akjj -1-qkii = 1qij对 k = j，n（2. 25）（3） 如果 j = n， 则停止；（4） 假定 qjj 0， 根据下式计算 Q 的第 j 行 （ 即矩阵 U 的第 j 行）：jkqjjikq= 1 （a- j -1 q q）对

27、k = j + 1，n（2. 26）jki = 1ji（5） 令 j = j + 1， 转到步骤 （2）。一旦完成 LU 分解， 通过 “ 前代” 就能得到哑向量 y：q-qyk = 1kk（bkk -1kj yjj = 1）对 k = 1，n（2. 27）n类似地， 通过 “ 回代” 就能得到解向量 x：xk = yk - j= k +1qkj xj对 k = n，n - 1，1（2. 28）衡量 LU 分解过程计算量的一种指标是， 获得解向量所需要的乘法和除法运算次数， 因为两者都是浮点数运算。 计算 Q 的第 j 列 （ L 的第 j 列） 需要的乘法和除法次数为nnj = 1 k =

28、j （ j - 1）类似地， 计算 Q 的第 j 行 （ U 的第 j 行） 需要的乘法和除法次数为 前代所需要的乘法和除法次数为回代所需要的乘法和除法次数为n-1n jj = 1 k = j +1nj = 1 jnj = 1 （ n - j）这样， 可以算出 LU 分解所需要的乘法和除法次数为31 （ n3 - n）而前代和回代所需要的乘法和除法次数为 n2 。 因此， 求解式 （2. 1） 线性方程组所需要的总的乘法和除法次数为31 （ n3 - n） + n2（2. 29）将这个乘除法次数与采用Cramer 法则所需要的2（ n + 1）！ 次乘除法相比， 对于任何较大阶数的矩阵， 显然

29、采用 LU 分解和前代/ 回代方法的计算效率要高得多。例 2. 3采用 LU 分解和前代/ 回代算法， 求解例 2. 2 中方程的解。解 2. 3第 1 步是求矩阵 A 的 LU 因子：1348 A =2123 4358 9274 从 j = 1 开始， 式 （2. 25） 表示 Q 的第 1 列与 A 的第 1 列完全相同。 类似地， 根据式 （2. 26）， Q 的第 1 行变为q12 = a12 = 3 = 3q111q13 = a13 = 4 = 4q111q14 = a14 = 8 = 8q111这样， 对于 j = 1， 矩阵 Q 变为1348 Q =249对于 j = 2， 将分

30、别计算矩阵 Q 主对角元下面的列元素和右面的行元素。 对于 Q 的第 2 列有q22 = q22 - q21 q12 = 1 - （2）（3） = - 5q32 = a32 - q31 q12 = 3 - （4）（3） = - 9q42 = a42 - q41 q12 = 2 - （9）（3） = - 25Q 中的每个元素通过 A 中的对应元素和 Q 中已经求得的元素得到。 注意乘积的第 2 个下标总是相同的， 而第 1 个下标与正要计算的元素下标相同。 这个规则对于列计算和行计算都是成立的。 Q 的第 2 行计算如下：q23 = q1 （ a23 - q21 q13 ） = 1 （2 - （

31、2）（4） = 622- 55q24 = q1 （ a24 - q21 q14 ） = 1 （3 - （2）（8） = 1322- 55对 j = 2 处理完以后， 矩阵 Q 变为1348 52- 5613 Q = 54- 99- 25继续计算 j = 3， Q 的第 3 列计算如下：6 1q33 = a33 - （ q31 q13 + q32 q23 ） = 5 - （4）（4） + （ - 9） 5 = - 5q43 = a43 - （ q41 q13 + q42 q23 ） = 7 -（9） （4） + （ - 25） 6 = 1而 Q 的第 3 行变为 15 q34 = q33 （ a

32、34 - （ q31 q14 + q32 q24 ）= （ - 5） 8 - （4）（8） + （ - 9） 13 = 3 5 得到1348 52- 5613 Q = 514- 9- 53 9- 251最后， 计算 j = 4， 最后一个对角元为q44 = a44 - （ q41 q14 + q42 q24 + q43 q34 ）= 4 - （9）（8） + （ - 25） 13 + （3）（1） = - 6这样， 5 1348Q = 613 2- 5515 54- 9-3 9- 251- 6 1000 2- 500 L = 4- 9- 109- 251- 6 51348 01613 U =

33、55 0013 0001 检查上述计算是否正确的一种方法是验证 LU 是否等于 A， 对于本例， 结果是成立的。1一旦求得了矩阵 LU， 下一步就是基于矩阵 L 和向量 b 通过前代来计算哑向量 y， 采用前代解方程 Ly = b 得到 y 为y =1=b1 L111 = 1y2 = （ b2 -LL21 y1 ） = （1 - （2）（1） = 122- 55y3 = （ b3 - （ L31 y1 + L32 y2 ） = （ - 5）1 - （4）（1） + （ - 9） 1 = 6L3344y4 = （ b4 - （ L41 y1 +LL42 y2 + L43 y3 ）5 =1 - （

34、9）（1） + （ - 25） 1 + （1）（6） = 这样- 6 5 1 1 y = 5 6 322 3 类似地， 采用回代解方程 Ux = y， 可以得到解向量 x 为x4 = y4 = 32 3 3x3 = y3 - U34 x4 = 6 - （3） 2 = 2x2 = y2 - （ U24 x4 + U23 x3 ） = 1- 13 3 + 6 3 = - 115x1 = y1 - （ U14 x4 + U13 x3 + U12 x2 ） 5 2 5 2 2= 1 - （8） 3 + （4） 3 + （3） - 11 = - 1最终得到解向量 2 2 - 1 2 2x = 1 - 1

35、1 2 3 3 此解与例 2. 2 中采用 Gauss 消去法和回代方法所得到的解是一样的。 检查所得到的解是否正确的一种快速方法是将解向量 x 代回到线性方程组 Ax = b 中， 看是否成立。2. 2. 1采用部分选主元的 LU 分解法10 - 101 x1 = 1 上述的 LU 分解过程假定了对角元是非零的。 实际上不但要求对角元是非零的， 还要求对角元与其他非零元具有相同的数量级。 考察如下线性方程组的解：2 21 x 5 （2. 30）通过观察， 很容易得出此线性方程组的解为x1 2x2 1而 A 的 LU 分解因子为L = 10 - 100 2（1 - 2 1010 ） U = 1

36、1010 采用前代方法求解哑向量 y 得到：01 10y1 = 10y =12（5 - 2 1010 ）（1 - 2 1010 ）x采用回代方法求解 Uy = x 得到2= y2 1x1 = 1010 - 1010 x2 0这里， x2 的解是正确的， 而 x1 的解是完全不对的。 为什么会发生这种情况呢？ 问题是式 （2. 30） 中的 10 - 10 对大多数计算机来说是太接近于零了。 但是， 如果将式 （2. 30） 进行整理， 变为21x52 10 - 101 x1 = 1 （2. 31）然后， 再进行 LU 分解， 则得到 20L = 110 - 10 1 - 2 10 - 10 11 U = 2 01 哑向量 y 变为y1 = 525- 10 = 1 - 2 101y21 1 - 2 10 - 10 通过回代， 得到 x 为x2 12x1 52- 1 （1） = 2这与通过观察得到的解相同。 因此， 即使对角元不是精确等于零， 将方程组的次序进行重新排列以使最大数值的元素位于对角元上， 仍然是一种好的做法。 这种做法被称为 “ 选主元”， 并导致了式 （2. 18） 中的置换矩阵 P。由于 Crout 算法从第 1 列和第