线性代数方程组的数值解法hyf

1、2斜面上物体的滑动斜面上物体的滑动m1m2m3已知：已知：m1=1.0kgm1=1.0kg，m2=2.0kgm2=2.0kg，m3=3.5kgm3=3.5kg，斜面倾角，斜面倾角=30=30；m1m1与斜面与斜面摩擦系数摩擦系数1=0.11=0.1， m1 m1与与m2m2摩擦系数摩擦系数2=0.82=0.8；不计绳和滑轮的质量与摩擦；不计绳和滑轮的质量与摩擦；求：作用在求：作用在m1m1与与m2m2上的正压力上的正压力F1F1和和F2F2，绳子的张力，绳子的张力F F，三个物体的加速度，三个物体的加速度a1, a2, a3a1, a2, a3；画出画出a1a1和和a2a2与斜面倾角与斜面倾角

2、 的关系曲线，分析在什么倾角范围内的关系曲线，分析在什么倾角范围内m2m2相对相对m1m1是是静止的？静止的？3斜面上物体的滑动斜面上物体的滑动FgmamgmFamFFgmFamgmFFgmF 31322222221111122211sinsincoscosm1m2m3oxyF2Fr2m2gm1gF2Fr2F1Fr1Fm3gF4线性方程组的矩阵解线性方程组的矩阵解FgmamgmFamFFgmFamgmFFgmF 31322222221111122211sinsincoscosbAx 0100000000100100011322121mmmA 2121aaFFFx gmgmgmgmgmb3221

3、1sincossincos5vclearvm1=1;vm2=2; vm3=3.5;vmu1=0.1;vmu2=0.8;vg=9.8;valpha=30*pi/180;vA=1 -1 0 0 0; -mu1, -mu2, 1, -m1, 0; 0 1 0 0 0; 0 mu2 0 0 -m2; 0, 0, 1, m3 0;vb= m1*g*cos(alpha); m1*g*sin(alpha); m2*g*cos(alpha); m2*g*sin(alpha); m3*g;vx=Ab6vclearvsyms F1 F2 F a1 a2;veq1 = F1 = F2 + m1*g*cos(alph

4、a);veq2 = F2 = m2*g*cos(alpha);veq3 = m1*a1 = F - m1*g*sin(alpha) - mu1*F1 - mu2*F2;veq4 = m2*a2 = mu2*F2 - m2*g*sin(alpha);veq5 = m3*a1 = m3*g - F;vF1, F2, F, a1, a2 = solve(eq1, eq2, eq3, eq4, eq5, F1, F2, F, a1, a2)vm1=1;vm2=2; vm3=3.5;vmu1=0.1;vmu2=0.8;vg=9.8;valpha=30*pi/180;vsubs(F)7 在科学计算中，经常

5、需要求解含有在科学计算中，经常需要求解含有n n个未知量个未知量 的的n n个方程构成的线性方程组个方程构成的线性方程组(4.1)(4.1) nnnnnnnnnnbxaxaxabxaxaxabxaxaxaLLLLLLLLLLLLLLL22112222212111212111 nnnnnnnnbbbbxxxxaaaaaaaaaAMMLLLLLLL2121212222111211,方程组还可以用矩阵形式表示为方程组还可以用矩阵形式表示为 bAx 8 当当 n n 较大时，较大时， 这个计算量是惊人的。这个计算量是惊人的。 例例 如，如， 当当 n= 20 n= 20 时，约需乘法次数为时，约需乘法

6、次数为若系数矩阵若系数矩阵A A非奇异，即非奇异，即则方程组有惟一解则方程组有惟一解 如果用每秒百万次的计算机来计算，需要三千万如果用每秒百万次的计算机来计算，需要三千万年时间。可见年时间。可见GramerGramer法则不是一种实用的方法。法则不是一种实用的方法。 因此，必须构造出适合于计算机使用的线性方程因此，必须构造出适合于计算机使用的线性方程组的求解方法。组的求解方法。0)det( ATnxxxx),(21L 根据根据 Gramer Gramer（克莱姆）法则，求解方程组（克莱姆）法则，求解方程组（4.14.1）时，要计算大量的行列式，所需乘法次数大约为时，要计算大量的行列式，所需乘法

7、次数大约为!)1(2nnN 20107.9 N detiiAxA9直接法直接法经有限步算术运算可求得方程组的精确解经有限步算术运算可求得方程组的精确解的方法（若计算过程无舍入误差）。的方法（若计算过程无舍入误差）。 已知的直接解法是克莱姆法则已知的直接解法是克莱姆法则(Gramer)。迭代法迭代法构造某种迭代格式，用其极限过程去逐步构造某种迭代格式，用其极限过程去逐步逼近线性方程组精确解的方法。逼近线性方程组精确解的方法。 求解线性方程组的数值方法可分为两大类：直接求解线性方程组的数值方法可分为两大类：直接方法和迭代方法。方法和迭代方法。10只有一列元素的矩阵只有一列元素的矩阵,21 naaa

8、BM称为列矩阵称为列矩阵( (或列向量或列向量).).全为零的方阵称为上三角矩阵。全为零的方阵称为上三角矩阵。即即主主对对角角线线以以下下元元素素形形如如）（3 nnnnaaaaaaLMMMLL00022211211O矩阵知识矩阵知识11 称为对角矩阵称为对角矩阵(或对角阵）或对角阵）. n LLLLLLL00000021形如形如 的方阵的方阵, ,全为零的方阵称为下三角矩阵。全为零的方阵称为下三角矩阵。即即主主对对角角线线以以上上元元素素形形如如 nnnnaaaaaaLMMMLL21222111000矩矩阵阵的的方方阵阵，即即既既是是上上三三角角又又是是下下三三角角记作记作).,(21ndi

9、agAL12222263422142 C22 16 32 816? )3(42)2(矩阵乘法矩阵乘法13定义线性方程组的初等变换是指下列三种变换定义线性方程组的初等变换是指下列三种变换 用一个非零的数乘某一个方程；用一个非零的数乘某一个方程； 将一个方程的倍数加到另一个方程上；将一个方程的倍数加到另一个方程上； 交换两个方程的位置交换两个方程的位置性质线性方程组经初等变换后，得到的线性性质线性方程组经初等变换后，得到的线性方程组与原线性方程组同解方程组与原线性方程组同解线性方程组的初等变换线性方程组的初等变换 14三角形方程组的解三角形方程组的解 对角形方程组对角形方程组若若aii0,则则 x

10、i=bi/aii, i =1,2,n。运算量运算量 O(n)。1111222222 nna xba xba xbL15 下三角方程组下三角方程组111121122221122112233 iiiiiinnnnnnnl xbl xl xbl xl xl xbl xl xl xl xbLLLL11122,11111From th equation, () /, 1,2, . iiiiiiii iiiijjjiiiijjiijil xbl xl xlxbl xxbl xlinLL16 上三角方程组,1111From th equation ()/ , /, 1,2,2,1.niiiii iiinni

11、ijjj innnnniiijjiij iiu xbuxu xbu xxb uxbu xuinn LL11112213311,111,111,1 nniiii iiinninnnnnnnnnnnu xu xu xu xbu xuxu xbuxuxbu xbLLLL17消元法的基本思想就是通过对方程组作初等变换，把一般形式消元法的基本思想就是通过对方程组作初等变换，把一般形式的线性方程组化为等价的易于求解的三角方程组。的线性方程组化为等价的易于求解的三角方程组。 高斯消元法（高斯消元法（Gaussian EliminationGaussian Elimination） 高斯消元法是一种古老的求解

12、线性方程组的方法，高斯消元法是一种古老的求解线性方程组的方法，它就是通过一系列的初等变换（消元），把线性方它就是通过一系列的初等变换（消元），把线性方程组化为等价的上三角方程组，然后通过回代方法程组化为等价的上三角方程组，然后通过回代方法求出原方程组的解。求出原方程组的解。1 Gauss消元法消元法 Gauss Gauss消元法由消元和回代两个过程组成。消元法由消元和回代两个过程组成。 nnnnnnnnnnbxaxaxabxaxaxabxaxaxaLLLLLLLLLLLLLLL221122222121112121111, 1,222221, 11212111 nnnnnnnnnnnaxaaxa

13、xaaxaxaxaLLLLLLLLLLLLL181 解：第二个方程乘以解：第二个方程乘以2 2，再与第一个方程对换次序，再与第一个方程对换次序得得第二个方程减去第一个方程的第二个方程减去第一个方程的2 2倍，倍，1.1 顺序顺序Gauss消元法消元法解线性方程组解线性方程组12311112322212323113250 xxxxxxxxx 12312312322313250 xxxxxxxxx 第三个方程减去第一个方程的第三个方程减去第一个方程的3 3倍，倍，得得 0-5231-1/2-1/21/21-312,bA用增广矩阵进行运算用增广矩阵进行运算 0-5231-3-122-111,bA19

14、第三个方程减去第二个方程的第三个方程减去第二个方程的5 5倍，得倍，得第三个方程乘以第三个方程乘以 ，得，得 13123232323526xxxxxxx -6-250-3-1102-111,bA1232332339xxxxxx 9300-3-1102-1-11,bA123233233xxxxxx 3100-3-1102-1-11,bA123(5,0,3).xxx123(5,0,3).xxx20这样，对于方程组这样，对于方程组或者或者 nnnnnnnnnnbxaxaxabxaxaxabxaxaxaLLLLLLLLLLLLLLL22112222212111212111bAx 我们用增广矩阵表示，并

15、给出我们用增广矩阵表示，并给出gauss消元法的具体算法消元法的具体算法 )1()1()1(3)1(2)1(1)1(3,n+1)1(3)1(33)1(32)1(31)1(2,n+1)1(2)1(23)1(22)1(21)1(1,n+1)1(1)1(13)1(12)1(11)1()1(,n,n+1nnnnnnnnaaaaaaaaaaaaaaaaaaaabAbALMMLMMMLLL21 )1()1()1(3)1(2)1(1)1(3,n+1)1(3)1(33)1(32)1(31)1(2,n+1)1(2)1(23)1(22)1(21)1(1,n+1)1(1)1(13)1(12)1(11)1()1(,n

16、,n+1nnnnnnnnaaaaaaaaaaaaaaaaaaaabALMMLMMMLLL）；（个个元元素素成成为为行行第第第第行行第第行行第第为为消消第第一一次次消消元元：用用1213213201111111121111111111 njniaaaaajiaainiaaijijijii,:,: ),()()()()()()()()()(LLLaii022 )2()2()2(3)2(2)2(3,n+1)2(3)2(33)2(32)2(2,n+1)2(2)2(23)2(22)1(1,n+1)1(1)1(13)1(12)1(11)2()2(000,n,n+1nnnnnnnaaaaaaaaaaaaaa

17、aaabALLLLLLLLL），；，（个个元元素素成成为为：行行第第，则则第第行行第第行行第第消消为为把把第第二二次次消消元元：用用）（）（）（）（）（）（）（）（）（13243 243 02222222232222222222 njniaaaaajiaainiaaijijijii,: ),(LLL23如此继续消元下去第如此继续消元下去第n1步结束后，得到矩阵步结束后，得到矩阵: )3()3()3(3)3(3,n+1)3(3)3(33)2(2,n+1)2(2)2(23)2(22)1(1,n+1)1(1)1(13)1(12)1(11)3()3(00000,n,n+1nnnnnnaaaaaaaaa

18、aaaaaabALLLLLLLLL）；（个个元元素素成成为为：行行第第，则则第第行行第第行行第第消消为为把把次次消消元元：用用第第）（）（）（）（）（）（）（）（）（11 k1 0k1kkkkk nkjnkiaaaaajiaainkiaakkkkikkkjkijkijkkkikik,: ),(LLL24这是与原线性方程组等价的方程组这是与原线性方程组等价的方程组. )()()3(3,n+1)3(3)3(33)2(2,n+1)2(2)2(23)2(22)1(1,n+1)1(1)1(13)1(12)1(11)()(000000,nn,n+1nnnnnnnnaaaaaaaaaaaaaabALLLLL

19、LLLL增广矩阵增广矩阵A（n），），b（n）对应上三角形方程组：对应上三角形方程组： )()()2(2)2(22)2(22)1(1)1(12)1(121)1(11nnnnnnnnnnbxabxaxabxaxaxaLLLLLLLLLLLLLLL25再进行回代求解，可以得到：再进行回代求解，可以得到： 1111211),.,(, nnkxaaaxaaxnkjjkjnkkkknnnnn261.消元过程：消元过程： 2.回代过程：回代过程： （k1，2，n1）第第k次消元执行计算：次消元执行计算：下面给出顺序下面给出顺序Gauss消元法的计算流程：消元法的计算流程：）；（）（）（）（）（）（11 1

20、 nkjnkiaaaaakkkkikkkjkijkij,LL对于对于kn1，n2，1，计算：，计算：先置：先置：nnnnnaax1 , nkjjkjnkkkkxaaax11127MATLAB程序程序 11. 消元消元A=2 -1 -3; 1/2, -1/2, -1/2; 3 2 -5; b=1 1 0;a=A,b; n=length(b);for k=1:n-1 % 第第k次消元次消元 for i=k+1:n % 第第i行行 m=a(i,k) /a(k,k); for j=k:n+1 % 第第j列列 a(i,j) = a(i,j)- a(k,j) *m; end endend 0-5231-1

21、/2-1/21/21-312,bA12311112322212323113250 xxxxxxxxx ）；（）（）（）（）（）（11 1 nkjnkiaaaaakkkkikkkjkijkij,LL28vx(n)=a(n,n+1)/a(n,n); %计计算算vfor k=n-1:-1:1 % 计算计算x(n-1:1)v s=0;v for j=k+1:nv s = s + a(k,j)*x(j);v endv x(k)=(a(k,n+1) - s )/a(k,k);vendvx2. 回代回代 3100-3-1102-1-11,bA)()(nnnnn，n+1naax , nkjjkjnkkkkxa

22、aax111(kn1，n2，1)29矩阵运算法程序矩阵运算法程序 2）；（）（）（）（）（）（11 1 nkjnkiaaaaakkkkikkkjkijkij,LL )1()1()1(3)1(2)1(1)1(3,n+1)1(3)1(33)1(32)1(31)1(2,n+1)1(2)1(23)1(22)1(21)1(1,n+1)1(1)1(13)1(12)1(11)1()1(,n,n+1nnnnnnnnaaaaaaaaaaaaaaaaaaaabALMMLMMMLLLvclearvA=2 -1 -3; 1/2 -1/2 -1/2; 3 2 -5; b=1 1 0;va=A,b; n=length(b

23、);vfor k=1:n-1v a(k+1):n, k:(n+1) = a(k+1):n, k:(n+1). v - a(k+1):n, k)*a(k, k:(n+1)/a(k,k); vendvx=zeros(n,1); x(n)=a(n,n+1)/a(n,n); vfor k=n-1:-1:1v x(k)=(a(k,n+1) - a(k, (k+1):n)*x(k+1):n) )/a(k,k);vendvx, nkjjkjnkkkkxaaax111)()(nnnnn，n+1naax 301.2 列主元列主元Gauss消元法消元法顺序顺序Gauss消去法计算过程中的消去法计算过程中的 akk

24、(k) 在第在第k步消元步消元时要用它作除数时要用它作除数,则可能会出现以下则可能会出现以下2种情况：种情况：若出现若出现 0，消元过程就不能进行下去。，消元过程就不能进行下去。 )(kkka 0，消元过程能够进行，但若 很小，也会造成舍入误差积累很大导致计算解的精度下降。 )(kkka)(kkka31 例例 4.1 4.1 解方程组解方程组精确解为：精确解为：x1=1/3x1=1/3，x2=2/3x2=2/3。数值计算取。数值计算取5 5位数位数字：字：解解 顺序消元：顺序消元： 交换方程的顺序：交换方程的顺序： 经高斯消元：经高斯消元：12120.00033.00002.00011.000

25、01.00001.0000 xxxx121220.00033.00002.0001 0 9999.06666.00.6667 xxxxx 12121.00001.00001.00000.00033.00002.0001xxxx122211.00001.00001.00000.6667 2.99971.99980.3333xxxxx32列主元列主元Gauss消元法思想消元法思想v在第在第k次消元前，将第次消元前，将第k列里的列里的n-k+1个元素中选取个元素中选取绝对值最大的一个作为主元素，并把此主元素所在绝对值最大的一个作为主元素，并把此主元素所在的行与第的行与第k行交换。行交换。v避免小数做

26、分母。避免小数做分母。1,22 1,22 1,222221, 11212111 nnnnnnninininnnnnnaxaxaaxaxaaxaxaaxaxaxaLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLL33下面将列主元下面将列主元Gauss消元法的计算步骤叙述如下：消元法的计算步骤叙述如下： 给定线性方程组 Axb, 记A（1）, b（1） A,b，列主元Gauss消去法的具体过程如下： 1. 首先在增广矩阵A（1）,b（1）第一列的n个元素中选取绝对值最大的一个作为主元素，并把此主元素所在的行与原第一行交换，即,max)1(11)1(1inikaa 2. 进行第

27、一步消元进行第一步消元,得到增广矩阵得到增广矩阵A（2）,b（2）。在矩阵。在矩阵A（2）,b（2） 第二列的第第二列的第 2:n个元素中选取绝对值最大的一个作个元素中选取绝对值最大的一个作为主元素，并把此主元素所在的行与原第二行交换，即为主元素，并把此主元素所在的行与原第二行交换，即,max)2(22)2(2inikaa )1(1)1()1(1)1(,bbaakjkj)2(2)2()2(2)2(,bbaakjkj34 3. 再进行第二步消元，得到增广矩阵A（3），b（3）。按此方法继续进行下去，经过n1步选主元和消元运算，得到增广矩阵A（n），b（n），它对应的方程组A（n）xb（n）是一个

28、与原方程组等价的上三角形方程组，可进行是一个与原方程组等价的上三角形方程组，可进行回代求解。回代求解。 容易证明，只要det（A）0，列主元Gauss消去法就可以顺利完成，即不会出现主元素为零或者绝对值太小的情形出现。35vA=2 -1 -3; 1/2 -1/2 -1/2; 3 2 -5; b=1 1 0;va=A,b; n=length(b);vfor k=1:n-1v val, p=max(abs(a(k:n, k); p=p+k-1;v a(k;p, :) = a(p;k, :);v a(k+1):n, k:(n+1) = a(k+1):n, k:(n+1). v - a(k+1):n,

29、 k)*a(k, k:(n+1)/a(k,k); vendvx=zeros(n,1); x(n)=a(n,n+1)/a(n,n); vfor k=n-1:-1:1v x(k)=(a(k,n+1) - a(k, (k+1):n)*x(k+1):n) )/a(k,k);vendvx列主元列主元Gauss消元程序消元程序 12311112322212323113250 xxxxxxxxx 36 由于这两种方法的精度差不多，且全主元Gauss消元法程序设计复杂占用机器时间较多，实际应用中一般采用列主元Gauss消元法，它既简单又能保证计算精度。 有时候在消元过程中可以在系数矩阵所有元素中选择绝对值最大

30、的元素作为主元素，这样的Gauss消元法叫做全主元Gauss消元法。 )1()1()1(3)1(2)1(1)1(3)1(3)1(33)1(32)1(31)1(2)1(2)1(23)1(22)1(21)1(1)1(1)1(13)1(12)1(11)1()1(,nnnnnnnnnbaaaabaaaabaaaabaaaabALMMLMMMLLL371.3. 1.3. 高斯高斯- -约当消元法约当消元法 高斯高斯-约当消元法约当消元法 就是将给定方程组通过加减就是将给定方程组通过加减消元化为对角形方程组的方法消元化为对角形方程组的方法. nnnnbxabxabxa22221111是无回代的消去法是无回

31、代的消去法, ,所进行的乘除运算次数比高斯消所进行的乘除运算次数比高斯消去法去法( (有回代的消去法有回代的消去法) )多多. .38解解 第第1步：将方程（步：将方程（2.1）乘上（）乘上（-3/2）加到方程（）加到方程（2.2） 将方程（将方程（2.1）乘上（）乘上（-1/2）加到方程（）加到方程（2.3） 则得到与原方程等价的方程组则得到与原方程等价的方程组2/5932/ 14231222321321321xxxxxxxxx（2.1）（2.3）（2.2）例例2用高斯用高斯-约当消元法解方程组约当消元法解方程组282112223232321xxxxxxx（2.4）（2.5）（2.6）（先消

32、下三角，后消上三角）（先消下三角，后消上三角）39第第2 步：消去（步：消去（2.4）和（）和（2.6）式中未知数）式中未知数 x2 ，得到等价的方程组，得到等价的方程组0, 1, 2/1321 xxx不用回代，即可求得原方程组的解为：不用回代，即可求得原方程组的解为：第第3 步：消去（步：消去（2.7）和（）和（2.8）式中未知数）式中未知数 x3 ，得等价的对角形方程组，得等价的对角形方程组 0114233231xxxxx（2.7）（2.8） 0112321xxx402 2 三角分解方法三角分解方法11112212311(1)(1)(1)(1)22222322(2)(2)(2)33333(

33、1)(1) nnnnnnnnnnnna xa xa xa xbaxaxaxbaxaxbaxbLLLLL 11112211211222221122 nnnnnnnnnna xa xa xba xa xa xba xa xa xbLLLLL）；（）（）（）（）（）（11 1 nkjnkiaaaaakkkkikkkjkijkij,LL利用初等变换消元：利用初等变换消元：4121111213141111213141(1)(1)(1)(1)112223242212223242(1)(1)(1)(1)31313233343323334311(1)(1)41424344442434411110 0 01 0

34、 000 1 0000 0 1aaaaabaaaabaaaabaaaabaaaaabaaabaaaaabaaaaa(1)(1)44b111213141111213141(1)(1)(1)(1)2223242212223242(1)(1)(1)(1)3132333433233343(1)(1)(1)(1)4142434444243444000aaaabaaaabaaabaaaabaaaabaaabaaaabaaab用矩阵相乘来解释高斯消元过程，取用矩阵相乘来解释高斯消元过程，取n=4。k1：2111111314110 0 01 0 0, where /,2,3,4.0 1 00 0 1iilLl

35、aaill(1)(1) ,AxbA xb等价于等价于 L1A b=A(1) b(1), 用矩阵运算的观点来看，消元的每一步计算等价于用用矩阵运算的观点来看，消元的每一步计算等价于用一个单位下三角矩阵左乘前一步约化得到的矩阵。一个单位下三角矩阵左乘前一步约化得到的矩阵。42k=2111213141111213141(1)(1)(1)(1)(1)(1)(1)(1)2223242222324232(1)(1)(1)(1)(2)(2)(2)22323334333343(1)(1)(1)(1)42424344422100 0010 00001 0000000 1aaaabaaaabaaabaaabaaa

36、aabaabaaaaba(2)(2)(2)4344400aab111213141111213141(1)(1)(1)(1)(1)(1)(1)(1)22232422223242(1)(1)(1)(1)(2)(2)(2)323334333343(1)(1)(1)(1)(2)(2)(2)42434444344400000000aaaabaaaabaaabaaabaaabaabaaabaab(1)(1)222223242100 0010 0, where /,3,4.01 000 1iiLlaaill(1)(1)(2)(2) ,A xbAxb等价于等价于 L2A(1) b(1)= A(2) b(2),

37、43k=3111213141111213141(1)(1)(1)(1)(1)(1)(1)(1)22232422223242(1)(1)(1)(2)(2)(2)333433334343(1)(1)(1)(2)(2)43444444331 0000 100000 01000000 0100000aaaabaaaabaaabaaabaabaabaaababa111213141111213141(1)(1)(1)(1)(1)(1)(1)(1)22232422223242(1)(1)(1)(2)(2)(2)3334333343(1)(1)(1)(2)(2)4344444400000000000aaaab

38、aaaabaaabaaabaabaabaabab(1)(1)33333431 0000 100, where /,4.0 0100 01iiLlaail(2)(2)(3)(3) ,AxbAxb等价于等价于 L3A(2) b(2) =A(3) b(3),4421111213141112122232423231223132333431143414243444423341221110 0 0100 01 0001 0 0010 00 10001 00 0100 1 00 0100 10 0 1aaaaabaaaaabaaaaaaabaaaaaabaaaaa 111213141(1)(1)(1)(1)

39、2223242(1)(1)(1)33343(1)(1)444000000aaaabaaabaabab111123 U=ALL L L1L2L3LAA*其中，其中，(2)(1)3323211111321123 , , AL AL L AL L L AAL L LAL L L ALUALU上三角矩阵上三角矩阵4511112311121313243414221313241421 00010 0 0100 00 1001 0 0010 0 0 0100 1 001 00 010 0 100 11 0 0 01 0 0 01 0 0 00 1 0 01 0 00 1 0 00 1 001 00 0 10

40、0 1LL L Llllllllllll 432131324142430 0 1 00 01100 010 0 1 01lllllll初等变换乘法因子初等变换乘法因子单位下三角矩阵单位下三角矩阵/ijijjjlaa 46由此可见，在顺序由此可见，在顺序Gauss消元法的过程中，系数矩阵消元法的过程中，系数矩阵A经过一系经过一系列单位下三角矩阵的左乘运算约化为上三角矩阵列单位下三角矩阵的左乘运算约化为上三角矩阵A，即，即 112111111121123 ., ., nnnALL L AALL LAL L LLALUALU 1111121323121111211nnnnnllllllLLLLLMM

41、L nnnnuuuuuuMLL22211211 ) 1() 1(2) 1(2211211nnnnaaaaaaUMLL从顺序从顺序Gauss消元法的矩阵运算表示式可知，系数矩阵消元法的矩阵运算表示式可知，系数矩阵A可分解可分解为一个单位下三角矩阵为一个单位下三角矩阵L和一个上三角矩阵和一个上三角矩阵U的乘积的乘积47第一个方程组的系数矩阵为下三角矩阵，第二个方程第一个方程组的系数矩阵为下三角矩阵，第二个方程组的系数矩阵为上三角矩阵，两个方程组都非常容易组的系数矩阵为上三角矩阵，两个方程组都非常容易求解，具体求解结果如下：求解，具体求解结果如下：我们将我们将称为矩阵称为矩阵A的三角分解的三角分解L

42、UA 对线性方程组：对线性方程组：bAX bLUX 令令YUX 则有则有 YUXbLY48bLY 对于对于 nnnnnnbbbbyyyyllllllMMLMM3213211213231211111由由 nkylbybykiikikk, 3 , 2,1111L解得解得49YUX 对于对于 nnnnnnyyyxxxuuuuuuMMMLL212122211211由由 1 , 2, 1,/1Lnniuxuyxuyxiinijjijiinnnn求得求得50三角分解法求解过程：三角分解法求解过程： nkylbybykiikikk, 3, 2,1111L 1 , 2, 1,/1Lnniuxuyxuyxiin

43、ijjijiinnnn YUXbLYbAX bLUX LUA 51由这个简单的计算过程可知，系数矩阵的三角分解很由这个简单的计算过程可知，系数矩阵的三角分解很关键，如何进行三角分解更容易？下面介绍几种方法。关键，如何进行三角分解更容易？下面介绍几种方法。 杜利特尔杜利特尔(Doolittle)分解：分解：L为单位下三角矩阵，为单位下三角矩阵， U为上三角矩阵；为上三角矩阵；库朗库朗(Courant)分解：分解： L为下三角矩阵，为下三角矩阵， U为单位上三角矩阵。为单位上三角矩阵。522.1 2.1 杜利特尔杜利特尔(Doolittle)(Doolittle)分解分解 设设 A A 的各阶主子

44、式不为零，的各阶主子式不为零，A A 分解为分解为 A=LU A=LU，即，即先计算先计算U U的行元素，后计算的行元素，后计算L L的列元素：的列元素：U U的第的第1 1行：行：111211112121222212221212100100 100nnnnnnnnnnnnaaauuuaaaluuaaalluL LLLLLLLL LL LLLL11111111101 00,0nrrrual uuLM53111211112121222212221212100100 100nnnnnnnnnnnnaaauuuaaaluuaaalluL LLLLLLLL LL LLLL121111111 00,0

45、,1,2,., .jnjjrrjjrjjuual uuuajnLMU的第的第1行：行：54111211112121222212221212100100 100nnnnnnnnnnnnaaauuuaaaluuaaalluL LLLLLLLL LL LLLL112121212111121211101 00,00/;nrrrual ull ulauLML的第的第1列：列： 11111,11111111101 00,0/,2,3,., .niirrii iiriiual ulll ulauinLLM55111211112121222212221212100100 100nnnnnnnnnnnnaaau

46、uuaaaluuaaalluL LLLLLLLL LL LLLLv再计算再计算U U的第的第2 2行元素；行元素；v计算计算L L的第的第2 2列元素；列元素；vv计算计算U U的第的第k k行元素：行元素：,1,.,kkk kk nuuu111111 (when ,0) (1) , ,1,., (4.17)nkkjkrrjkrrjkrrrkkrrjkjkkrkkjkjkrrjral ul urk ll uulual ujk kn56LU分解的运算量：分解的运算量：v计算计算L L的第的第k k列元素：列元素：1,2,1,2,.,1:kkkkn klllkn111111 ( ,0) /, 1,

47、2,., (4.18)nkkikirrkirrkirrkikkkrkrrrkikikirrkkkral ul ul ul urk ulal uuikkn321().3nO n111211112121222212221212100100 100nnnnnnnnnnnnaaauuuaaaluuaaalluL LLLLLLLL LL LLLL5711 , ,1,.,kkjkjkrrjrual ujk kn11/, 1,2,.,kikikirrkkkrlal uuikkn1) A=LU2) Ly=b3) Ux=y11 , 1,2, (4.19)iiiijjjybl yinL1/ , (4.20)/,

48、1,2,1nnnnniiijjiij ixb uxyu xuin L杜利特尔直接分解算法：杜利特尔直接分解算法：58vclearva=2 1 1; 1 3 2; 1 2 2;vn=length(a); vL=eye(n); U=zeros(n); vfor k=1:n %计算计算U的第的第 k 行元素和行元素和L的第的第 k 列元素列元素v for j=k:n %计算计算U的第的第 k 行元素行元素v s=0;v for r=1:k-1v s=s+L(k,r)*U(r,j);v endv U(k,j)=a(k,j)-s;v endv for i=k+1:n %计算计算L的第的第 k 列元素列元

49、素v s=0;v for r=1:k-1v s=s+L(i,r)*U(r,k);v endv L(i,k)=(a(i,k)-s)/U(k,k);v endvendvL,U11 , ,1,.,kkjkjkrrjrual ujk kn11/, 1,2,.,kikikirrkkkrlal uuikkn2111321221000.5100.5 0.6121102.51.5000.6杜利特尔杜利特尔MATLAB程序程序1：LU11 , ,1,.,kkjkjkrrjrual ujk kn11/, 1,2,.,kikikirrkkkrlal uuikknA=59vclearvA=2 1 1; 1 3 2;

50、1 2 2;vn=length(A); vU=zeros(n); L=eye(n);vU(1,:)=A(1,:); L(2:n, 1)=A(2:n, 1)/U(1,1);vfor k=2:nv U(k, k:n) = A(k, k:n) - L(k, 1:k-1)*U(1:k-1, k:n);v L(k+1:n, k)=(A(k+1:n, k) - L(k+1:n, 1:k-1)*U(1:k-1,k)/U(k,k);vendvL,U11 , ,1,.,kkjkjkrrjrual ujk kn11/, 1,2,.,kikikirrkkkrlal uuikkn杜利特尔杜利特尔MATLAB程序程序2

51、：11 , ,1,.,kkjkjkrrjrual ujk kn11/, 1,2,.,kikikirrkkkrlal uuikkn60vclearvA=2 1 1; 1 3 2; 1 2 2;vn=length(A); vA(2:n, 1)=A(2:n, 1)/A(1,1);vfor k=2:nv A(k, k:n) = A(k, k:n) - A(k, 1:k-1)*A(1:k-1, k:n);v A(k+1:n, k) = (A(k+1:n, k) - A(k+1:n, 1:k-1)*A(1:k-1,k)/A(k,k);vendvA 2110.52.51.50.50.60.621113212

52、2A=A=L+U杜利特尔杜利特尔MATLAB程序程序3：11 , ,1,.,kkjkjkrrjrual ujk kn11/, 1,2,.,kikikirrkkkrlal uuikkn61v A=2 1 1; 1 3 2; 1 2 2;v L, U=lu(A)vL =v 1.0000 0 0v 0.5000 1.0000 0v 0.5000 0.6000 1.0000vU =v 2.0000 1.0000 1.0000v 0 2.5000 1.5000v 0 0 0.6000MATLAB程序程序4：211132122A=62例3.2 用杜利特尔分解求解方程组：解 设 A=LU，即 1231231

53、23242 1 1326 1 3 2 .2251 2 2xxxxxxAxxx1112132122233132332 1 110 01 3 21 00,1 2 2100uuuluullu100211 0.5100 2.5 1.5 .0.5 0.6 1000.6LU63解下三角方程组 Ly = b，即解上三角方程组 Ux= y，即Thus . Let , then .AxbLUxbUxyLyb11223310044 0.5106 4 0.5 0.6 150.6yyyyyy 132231211410 2.5 1.54 1000.60.61xxxxxx642.2 库朗分解 先计算L的第一列，再计算U的

54、第一行，其余类此。类似于杜利特尔分解法，得到：对k=1，2，n111211112121222212221212001001 (4.21) 001nnnnnnnnnnnnaaaluaaaalluaaalllL LLLLLLLL LL LLLL1111, ,1,., (4.22) /, 1,., (4.23)kikikirrkrkkjkjkrrjkkrlal uik knual uljkn65vclear allvA=2 -1 1; 4 2 1; 2 1 2;vn=length(A); vA(1,2:n)=A(1,2:n)/A(1,1);vfor k=2:nv A(k:n,k) = A(k:n,k

55、) - A(k:n,1:k-1)*A(1:k-1,k);v A(k,k+1:n) = (A(k,k+1:n) - A(k,1:k-1)*A(1:k-1,k+1:n)/A(k,k);vendvA库朗程序库朗程序11111, ,1,., /, 1,.,kikikirrkrkkjkjkrrjkkrlal uik knual uljkn66注：为保证LU分解能够进行，要求 ，既要求A的所有顺序主子式不为零，而线性方程组有解，只须|A| 0即可。为取消此限制，采用类似列主元消元法处理。库朗列主元直接分解算法输入：n，A，b计算L的第 k 列元素和U的第 k 行元素0,0iiiiul11for 1 : f

56、or : ( -th rank of L) |=max| (choose main elements) change -th row witkikikirrkrukikk i nkniknklal ullk h -th row, ( , ) enduA b6711 for 1 : ( -th row of U) / endkkjkjkrrjkkrjknkual ul68库朗列主元程序库朗列主元程序2vclear allvA=2 -1 1; 4 2 1; 2 1 2;vn=length(A); L=zeros(n); U=eye(n);vL(:,1)=A(:,1); U(1,2:n)=A(1,2

57、:n)/L(1,1);vfor k=2:nv L(k:n,k) = A(k:n,k) - L(k:n,1:k-1)*U(1:k-1,k);v val, p=max(abs(L(k:n, k); p=p+k-1;v L(k;p, :) = L(p;k, :) v U(k,k+1:n) = (A(k,k+1:n) - L(k,1:k-1)*U(1:k-1,k+1:n)/L(k,k);vendvL,U1111, ,1,., /, 1,.,kikikirrkrkkjkjkrrjkkrlal uik knual uljkn69例2.3 用库朗分解求解方程组：解 设 A=LU，即 12321 154 2

58、14 .2 1 25xxx 11121321222331323321 10014 2 1001,2 1 2001luullulll2 0011/ 2 1/ 2 4 40, 011/ 4 .2 2 3/ 2001LU70解下三角方程组 Ly = b，即解(单位)上三角方程组 Ux= y，即Thus . Let , then .AxbLUxbUxyLyb1122332 0055/ 2 4 404 3/ 2 2 2 3/ 252 yyyyyy 13223111/ 2 1/ 25/ 22011/ 43/ 2 100121xxxxxx 71解线性代数方程组的迭代法4.3 4.3 雅可比迭代雅可比迭代4.

59、4 4.4 高斯高斯- -赛德尔赛德尔(Gauss-Seidel)(Gauss-Seidel)迭代迭代72非线性方程迭代法：非线性方程迭代法：f(x) = 0f(x) = 0转化为等价的转化为等价的x=(x), x=(x), 构造构造xk+1=(xk), xk+1=(xk), 由由x0 x0计算序列计算序列xkxk的极限即的极限即x x* *由初始向量由初始向量x(0)计算向量序列计算向量序列 x(k) ，若收敛到某个极限，若收敛到某个极限向量向量 x*，即，即 则则 x* 就是线性方程组的解。就是线性方程组的解。( )*,kxx我们可以对线性方程组我们可以对线性方程组Ax=y进行等价变换，构

60、造出等价方进行等价变换，构造出等价方程程组组 x = Mx+g,由此构造迭代关系式由此构造迭代关系式:(1)( ), 0,1,2,kkxMxgkL73若向量序列若向量序列 x(k) 收敛到某个极限向量收敛到某个极限向量 x*，即，即 则则 x* 就是线性方程组的解。就是线性方程组的解。( )*,kxx1111 () , (,),NP xyNxPxyxNPxNyMxgMNP gNy(0)(0)(0)(0)T12(1)(2) (,) , ,nxxxxxxLL取初始向量由迭代得到向量序列例如，分解例如，分解A=N-P，则，则(1)( ), 0,1,2,kkxMxgkLAx=y74 4.3 雅可比迭代