本节书摘来自华章出版社《数值分析（原书第2版）》一 书中的第2章，第2.4节，作者：（美）Timothy Sauer，更多章节内容可以访问云栖社区“华章计算机”公众号查看。

2.4　PA=LU分解

由于两个严重的缺陷：0主元以及淹没问题，之前我们考虑的高斯消去法称为“朴素”问题.对于一个非奇异矩阵，二者通过改进算法都可以避免.改进策略的重点是交换系数矩阵的行，该方法被称为部分主元.

2.4.1　部分主元

为了使用经典高斯消去法处理n个方程n个未知变量问题，第一步使用对角线元素a11作为主元消去第一列.部分主元包含在每一步消去步骤之前的比较.找到第一列中最大的一个元素，其对应行和主元行进行交换，在当前情况下，主元行是第一行.
换句话说，在高斯消去开始后，部分主元方法要求我们选择第p行，其中ap1≥ai1，　1≤i≤n(2.21)交换第1行和第p行.然后和平常一样，消去第1列，这时候使用的是新版本的a11作为主元.用于消去ai1的乘子如下：mi1=ai1a11mi1≤1.
在算法中，每一次选择主元时都需要进行相同的检查.当决定第2个主元时，我们使用当前的a22并检查在该元素下方的所有元素.我们选择第p行满足ap2≥ai2其中2≤i≤n，如果p≠2交换第2行和第p行.在本步计算中不再包含第1行.如果a22本身就是最大值，则不需要进行行交换.
在消去过程中，对每一列都实施相同的策略.在消去第k列时，找到第p行，k≤p≤n，定位最大的apk，必要时交换第p行和第k行，然后继续进行消去.注意到使用部分主元方法保证所有乘子，或者L的元素的绝对值不大于1.通过这种对于高斯消去方法所进行的小的改变，例2.13中出现的淹没问题可以完全避免.
例2.14　使用高斯消去的部分主元方法求解（2.1）方程组.
方程可以写成如下的表格形式：11|3
3-4|2根据部分主元方法，我们比较a11=1和它下面的所有元素，在这种情况下，唯一的元素a21=3.由于a21＞a11，我们必须交换第1行和第2行.新的表格如下3-4|2
11|3→第2行减去第1行的13倍→3-4|2
073|73
95在回代之后，对应的解是x2=1，x1=2，和我们前面计算的一致.当我们第一次求解这个问题时，乘子是3，但是使用部分主元策略时，则不会出现这样的情况.
例2.15　使用高斯消去的部分主元方法求解方程组:x1-x2+3x3=-3
-x1-2x3=1
2x1+2x2+4x3=0这个例子可以写为表格形式1-13|-3
-10-2|1
224|0在部分主元方法中，我们比较a11=1与a21=1以及a31=2，选择a31作为新的主元.这通过交换第1行和第3行实现:1-13|-3
-10-2|1
224|0→交换第1行和第3行→224|0
-10-2|1
1-13|-3
→第2行减去第1行的-12倍→224|0
010|1
1-13|-3
→第3行减去第1行的12倍→224|0
010|1
0-21|-3在消去第2列时，我们必须比较当前的a22和当前的a32.由于后者更大，我们再一次进行行交换:224|0
010|1
0-21|-3→交换第2行和第3行→224|0
0-21|-3
010|1
→第3行减去第2行的-12倍→224|0
0-21|-3
0012|-12注意到所有的乘子绝对值都比1小.
当前的方程非常容易求解，从如下方程12x3=-12
-2x2+x3=-3
2x1+2x2+4x3=0我们得到其对应解为x=［1，1，-1］.
注意到部分主元方法也可以解决0主元问题.当遇到一个潜在的0主元，比如说，如果a11=0，它立刻被当前列中的一个非零元素替换.如果在当前主元位置以及下方的任何位置没有非零元素，则矩阵是奇异矩阵，在这种情况下以任何方式，高斯消去方法都不能得到解.96

2.4.2　置换矩阵

在讲述LU分解方法可以利用行交换来处理高斯消去问题之前，我们将讨论置换矩阵的主要性质.
定义2.7　置换矩阵是一个n×n的矩阵，其在每一行、每一列仅有一个1，其他全部为0.
等价地，置换矩阵P可以通过对n×n的单位矩阵应用任意的行交换（或者任意的列交换）得到.例如10
01,01
10是仅有的两个2×2的置换矩阵，100
010
001,010
100
001,100
001
010
001
010
100,001
100
010,010
001
100是六个3×3的置换矩阵.
下一个定理描述在一个矩阵的左侧乘上置换矩阵会发生的情况.
定理2.8（置换矩阵的基础定理）　令P是通过对单位矩阵实施一组特定的行交换后得到的一个n×n的置换矩阵.则对于任意的n×n矩阵A，PA对应于对矩阵A实施同样的行交换得到的结果.
例如，置换矩阵100
001
010通过交换单位矩阵的第2行和第3行后得到.在矩阵的左侧乘上这样的置换矩阵P相当于交换矩阵的第2行和第3行:100
001
010abc
def
ghi=abc
ghi
def记住定理2.8的较好的方式是想象把P和单位矩阵I相乘:100
001
010100
010
001=100
001
010有两种方式来看上面的这种等价:首先，作为单位矩阵相乘（我们在右侧得到置换矩阵）;其次，作为置换矩阵对于单位矩阵的行所进行的交换.定理2.8的内容是通过和矩阵P相乘得到的行置换，与构造P所进行的行置换相同.97

2.4.3　PA=LU分解

在本节中，我们把关于高斯消去所知道的一切组合起来，进行PA=LU分解.这是部分主元的高斯消去的矩阵形式.PA=LU分解是求解线性方程组的主要方法.
正如名字中所暗示的，PA=LU是对于矩阵A包含行交换的LU分解.在部分主元中，开始的时候我们并不知道需要进行置换的列，所以我们必须谨慎地把行置换的信息加入分解中，我们需要记录高斯消去法中所有的主元.首先，我们来看一个例子.
例2.16　找出矩阵A的PA=LU分解.A=215
44-4
131首先，根据部分主元，第1行和第2行需要进行交换:P=010
100
001
215
44-4
131→交换第1行和第2行→44-4
215
131我们将使用置换矩阵P记录所有进行的行置换.现在我们进行两个行操作即→第2行减去第1行的12倍→44-4
-17
131
→第3行减去第1行的14倍→44-4
-17
22消去第1列.我们现在做得和以前不一样，没有仅仅在消去位置放上0，而是把0作为存储位置.在位置（i，j）的每个0里，保存用于消去该位置元素的乘子mij.这样做是有原因的.使用该机制，当使用更多的行交换，乘子也可以和其对应的行在一起.
然后必须通过比较选择第2个主元.由于a22=1＜2＝a32，在消去第2列之前必须进行行置换.注意到先前的乘子随着置换出现了移动:P=010
001
100→交换第2行和第3行→44-4
22
-17最后，随着另一次的行操作，消去过程结束了:→第3行减去第2行的-12倍→44-4
22
8

98这是最后一次消去.现在我们可以输出PA=LU分解:010
001
100215
44-4
131=100
1410
12-12144-4
022
008
P　　　　　　A　　　　　　L　　　　　　　U(2.22)L矩阵的元素在下三角矩阵（对角线以下）的0中，U则对应上三角矩阵.最终（累计）的置换矩阵保存为P.
使用PA=LU求解方程组Ax=b，这是A=LU的一个小的变体.在方程Ax=b两侧左乘P，并如前所述进行处理:PAx=Pb
LUx=Pb(2.23)求解1.Lc=Pb得到c
2.Ux=c得到x（2.24）如前所述，计算的重点是计算中代价最大的部分。在当前问题中，该部分计算是要确定的PA=LU，这部分计算不需要知道b.由于是对系数矩阵的行置换PA做LU分解，那么在做回代之前有必要对于右侧的向量b做相同的置换.该置换过程可以通过在回代开始计算Pb得到.看到高斯消去的矩阵形式的值很明显：所有消去和主元的细节都自动保存在矩阵方程中.
例2.17　使用PA=LU分解，求解系统Ax=b，其中A=215
44-4
131,b=5
0
6PA=LU分解可以从（2.22）得到.还需要进行两步回代.
1.Lc=Pb:100
1410
12-121c1
c2
c3=010
001
1005
0
6=0
6
5从顶部开始，我们得到c1=0
14(0)+c2=6c2=6
12(0)-12(6)+c3=5c3=82.Ux=c:44-4
022
008x1
x2
x3=0
6
8

99从底部开始，8x3=8x3=1
2x2+2(1)=6x2=2
4x1+4(2)-4(1)=0x1=-1(2.25)因而，解是x=［-1，2，1］.
例2.18　使用部分主元的PA=LU分解，求解系统2x1+3x2=4,3x1+2x2=1.
在矩阵形式中，对应方程如下23
32x1
x2=4
1首先我们忽略右侧的b.根据部分主元,必须交换第1行和第2行（由于a21>a11）.消去过程是P＝01
10　　　　　　　　　　　　　　　　　　　
A=23
32→交换第1行和第2行→32
23
→第2行减去第1行的23倍→32
53因而，PA=LU分解如下01
1023
32=10
23132
053
P　　　A　　　L　　　U第一次回代Lc=Pb如下10
231c1
c2=01
104
1=1
4从顶部开始，我们得到c1=1
23(1)+c2=4c2=103第二次回代Ux=c如下32
053x1
x2=1
103从底部开始，我们得到53x2=103x2=2
3x1+2(2)=1x1=-1(2.26)因而，对应的解是x=［-1，2］.100
每一个n×n矩阵都具有一个PA=LU分解.我们简单遵循部分主元原则，如果得到的主元是0，这意味着所有需要消去的元素都已经是0了，因而这一列不需要任何操作.
到目前为止，所有描述的技术在MATLAB中都已经实现.我们所讨论过的关于高斯消去法最复杂的形式是PA=LU分解.MATLAB的lu命令接受方块系数矩阵A并返回P、L和U.后面的MATLAB脚本定义了例2.16中的矩阵并计算其对应的分解:

2.4节习题
1.找出下面矩阵的PA=LU分解（使用部分主元）:
(a) 13
23　　　　　　　(b) 24
13　　　　　　　(c)15
512　　　　　　　(d) 01
10
2.找出下面矩阵的PA=LU分解（使用部分主元）:
(a) 110
21-1
-11-1(b) 013
211
-1-12(c) 12-3
242
-103(d) 010
102
-210
3.通过PA=LU分解求解系统，并实施两步回代.
(a) 37
61x1
x2=1
-11(b) 312
634
315x1
x2
x3=0
1
3
4.通过PA=LU分解求解系统，并实施两步回代.
(a) 420
442
223x1
x2
x3=2
4
6(b) -101
211
-120x1
x2
x3=-2
17
3

101
5.写出一个5×5的矩阵P，在另一个矩阵左侧乘上P会使其对应的第2行和第5行交换.
6.（a） 写出一个4×4的矩阵P，在另一个矩阵左侧乘上P会使其对应的第2行和第4行交换.
（b） 如果在右边乘上P会得到什么结果？用一个例子来证实.
7.改变最左边矩阵的四个元素使得矩阵乘法等式成立:0000
0000
0000
00001234
3456
5678
7890=5678
3456
7890
12348.找出习题2.3.15中的矩阵A的PA=LU分解.所需的最大乘子lij是多少？
9.（a） 找出矩阵A=1001
-1101
-1-111
-1-1-11的PA=LU分解.
（b）令A是一个与（a）中形式一致的n×n矩阵.描述PA=LU分解矩阵的每一个元素.
10.（a） 假设A是一个n×n矩阵，元素aij≤1，其中1≤i，j≤n.证明:PA=LU分解中的矩阵U对于所有的1≤i，j≤n，满足uij≤2n-1.见习题9（b）.（b）对于任意n×n矩阵A，表达并证实类似的事实.
事实验证2　欧拉伯努利横梁
欧拉伯努利横梁是材质在应力下扭曲的一个基本模型.通过离散化可以把微分系统转化为线性方程组.离散尺度越小，得到的方程组系统就越大.这个例子会给我们提供一个对于系统规模和科学计算中病态问题的非常有趣的实例研究.
在一个长度为L的横梁上，横梁垂直方向的偏移表示为函数y（x），其中0≤x≤L.在计算中我们将使用MKS单位:米、千米、秒.偏移y（x）满足欧拉伯努利方程EIy=f(x)(2.27)其中E是材质的杨氏模量，I是面积惯性模量.E和I对于横梁都是常数.右侧的f（x）是加载负荷，包括横梁重量，该重量以单位长度上受的力来表示.
微分方程的离散化技术将在第5章进行讲述，其中对于一个小的增量h，四阶导数的合理近似是y(x)≈y(x-2h)-4y(x-h)+6y(x)-4y(x+h)+y(x+2h)h4(2.28)离散化的近似误差和h2成比例（见习题5.1.21）.在我们的策略中，横梁被表示为许多长度为h的分段的组合，并对于每个分段应用微分方程的离散形式.102
对于正数n，令h=L/n.考虑均匀分割的格子0=x0跳水板就是一个横梁，它的一端系在支撑物上，而另一端是自由端.这被称为钳制自由横梁或者有时称为悬臂梁.对于钳制（左）端和自由（右）端的边界条件如下：y(0)=y′(0)=y″(L)=y(L)=0特别是，y0=0.但是，注意到计算y1为我们提出一个问题，对于微分方程（2.27）在点x1应用（2.29）近似得到y-1-4y0+6y1-4y2+y3=h4EIf(x1)(2.30)其中y-1没有定义.我们必须在点x1接近钳制端时，使用不同的导数近似.习题5.1.22（a）导出如下近似y(x1)≈16y(x1)-9y(x1+h)+83y(x1+2h)-14y(x1+3h)h4(2.31)当y(x0)=y′(x0)=0时该方程有效.
当前把近似称为“有效”，意味着离散化的近似误差和h2成正比，这与方程（2.28）中的情况一致.理论上，这意味着该方式下近似导数的误差，对于足够小的h，会降低到0.该概念是第5章中讨论数值微分的焦点.对于我们来说的结果是，我们可以使用（2.31）的近似，在i=1处考虑端点条件，得到16y1-9y2+83y3-14y4=h4EIf(x1)横梁自由的右端的计算稍微复杂一些，这是由于我们必须在整个横梁上计算yi.再一次，我们需要在最后的两个点xn-1和xn使用不同的导数近似.习题5.1.22给出如下近似y(xn-1)≈-28yn+72yn-1-60yn-2+16yn-317h4(2.32)
y(xn)≈72yn-156yn-1+96yn-2-12yn-317h4(2.33)

103　　当y″（xn）=y(xn)=0时该方程有效.
现在我们可以写出关于跳水板的n个等式和n个未知变量的方程组.该矩阵方程包含了所有对于微分方程（2.27）在点x1，…，xn上的近似，具有h2的精度:16-983-14
-46-41
1-46-41
1-46-41

1-46-41
1-46-41
1617-60177217-2817
-12179617-156177217y1
y2



yn-1
yn=h4EIf(x1)
f(x2)



f(xn-1)
f(xn)(2.34)（2.34）中的结构矩阵A是一个带状矩阵，意味着所有和对角线足够远的元素的值都是0.特别是矩阵元素aij除了i-j≤3之外的值都为0.该带状矩阵的带宽是7，由于i-j具有7个非零的元素aij.
最后，我们可以对钳制自由矩阵进行建模.我们考虑一个由道格拉斯杉制成的实木跳水板.假设跳水板的长度是L=2米，30厘米宽，3厘米厚.道格拉斯杉的密度大约是480kg/m3.1牛顿力是1kg·m/s2，该木材的杨氏模量大约是E=1.3×1010帕斯卡，或者牛顿/米2.横梁中心附近面积的惯性模量是wd3/12，其中w是横梁宽度，d是横梁的厚度.
开始时将计算没有负荷时横梁的偏移，因而f（x）仅仅表示横梁自身的重量，单位是力每米.因而f（x）每米的质量480wd乘上向下的重力加速度-g=-9.81m/s2，或者常量f（x）=f=-480wdg.读者应该检查一下（2.27）两侧的单位的匹配.f是常数时（2.27）有解析解，因而可以检查计算的结果.
在检查未加负荷的横梁代码之后，将进一步对两种情况进行建模.第一种情况下，一个正弦负载（或者“堆”）会被加在横梁上.在这种情况下，仍然有一个已知的解析解.但是没有精确的导数近似，因而，你可以对误差进行建模，它是关于格子大小h的函数.对于大的n，看一下条件问题的效果.然后，你将把跳水运动员放在横梁上.
建议活动
1.写出定义（2.34）中结构矩阵A的MATLAB程序.然后，使用MATLAB的＼命令或者自己设计的代码，求解系统，得到偏移yi，使用n=10作为网格大小.
2.同时画出步骤1中的解和正确解，y(x)=(f/24EI)x2(x2-4Lx+6L2)，其中f=f（x）是上面定义的常数.在横梁末端检查误差，此时x=Lm.在这种简单情况下，导数近似是正确的，因而你的误差接近机器舍入的结果.
3.返回步骤1中的计算，对于n=10·2k，其中k=1，…，11.对于每个n在x=L时的误差制表，对于哪一个n这个误差最小？为什么这个误差在某一点后随着n增加？你可以再104做一张关于A的条件数的表，该条件数对应一个关于n的函数，并有助于对最后一个问题的回答.对于一个大的k实施这个步骤，你可能需要让MATLAB把矩阵A保存为一个稀疏矩阵以避免完全消耗内存.使用命令A=sparse（n，n）初始化A，就可以得到稀疏矩阵，其他操作和前面一样.我们将在下一节中更加详细地讨论稀疏矩阵.
4.对横梁加上正弦压力.这意味着给受力项加上一个形如s(x)=-pgsinπLx的函数.证明解y(x)=f24EIx2(x2-4Lx+6L2)-pgLEIπL3π3sinπLx-x36+L2x2-L2π2x满足欧拉伯努利横梁方程以及钳制自由边界条件.
5.重新运行第3步中对于正弦负载的计算.（确保包含了横臂自身的重量.）设p=100kg/m，同时画出你的计算结果和正确的解.回答步骤3中的问题，并回答下面的问题:在x=L点的误差和前面生成的h2是否成比例？你可以在log-log图上画出误差与h，来探求该问题.条件数起作用了吗？
6.现在移除正弦负载并在横梁上加上70kg跳水员，在横梁的最后20cm保持平衡.你必须对f（xi）在单位长度加上力-g乘上70/0.2kg/m，1.8≤xi≤2，使用步骤5中找到的最优的n值，并再次求解问题.画出图并找出跳水板在自由端的偏斜度.
7.如果我们把跳水板的自由端也固定，就得到钳制钳制横梁，两端的边界条件相同y(0)=y′(0)=y′(L)=0.它用于对下陷结构，例如桥梁，进行建模.从有微小差异的均匀步长格点开始0=x0＜x1<…＜xn＜xn+1=L，其中h=xi-xi-1，i=1,…,n，找出具有n个方程n个未知变量的系统，用于确定y1,…,yn.（这和钳制自由版本的模型类似，除了系数矩阵A的最后两行应该是前两行的反转.）对于正弦负载进行求解，并对于横梁中心x=L/2回答步骤5中的问题.在正弦负载下的钳制钳制横梁的精确解是y(x)=f24EIx2(L-x)2-pgL2π4EIL2sinπLx+πx(x-L)8.更多探索的想法:如果跳水板的宽度加倍，跳水板的偏移该如何变化？如果跳水板的厚度加倍这种变化会更多还是更少？（两个板的质量相同.）如果横截面是圆形或者环形，但是面积和矩形截面相同，最大偏移量会如何变化？（面积惯性模量对于半径为r的圆形截面是I=πr4/4，对于内半径为r1外半径为r2的环形截面的面积惯性模量是I=π(r42-r41)/4.）找出对于I横梁的面积惯性模量.不同材质的杨氏模量已知并已经被制表.例如，钢材的密度是7850kg/m3，它对应的杨氏模量大约是2×1011帕斯卡.
欧拉伯努利是一个相对简单、经典的模型.更多最近的模型考虑更加奇异的扭曲，例如Timoshenko横梁，其中横梁的截面与横梁的主轴可能不垂直.105