求解线性方程组直接解法.docVIP

求解线性方程组的直接解法 5.2?????LU分解 ①?? Gauss消去法实现了LU分解 顺序消元结束时的上三角矩阵U和所用的乘数,严格下三角矩阵。将下三角矩阵的对角元改成1,记为L,则有A=LU, 这事实是一般的,我们不难从消去的第k个元素时的矩阵k行及k列元素的历史得到这一点.因为从消元的历史有 ukj=akj-mk1u1j- mk2u2j -…- mk,k-1uk-1,j, j=k,k+1,…,n mik=(aik-mi1u1k- mi2u2k -…-mi,k-1uk-1,k)/ukk i=k+1,k+2,…,n 于是 akj=mk1u1j+mk2u2j+…+mk,k-1uk-1,j+ukj, j=k,k+1,…,n aik=mi1u1k+mi2u2k+…+mi,k-1uk-1,k+mikukk i=k+1,k+2,…,n 从前面两个式子我们可以直接计算L和U(见下段).将矩阵分解为单位下三角矩阵和上三角矩阵之积称为矩阵的LU分解.顺序消元实现了LU分解,同时还求出了g, Lg=b的解. ② ? 直接LU分解 上段我们得到(lij=mij) ukj=akj-lk1u1j-lk2u2j -…- lk,k-1uk-1,j, j=k,k+1,…,n lik=(aik-li1u1k-li2u2k -…-li,k-1uk-1,k)/ukk i=k+1,k+2,…,n 二式也可从A=LU的n2个等式解出.下面以n=3为例说明. u11=a11 u12=a12 u13=a13 l21=a21/u11 u22=a22-l21u12 u23=a23-l21u13 l31=a31/u11 l32=(a32-l31u12)/ u22 u33=a33-l31u13- l32u23 从表上看到每个元素由所在位置的元素减去同行L左边诸元素与上方U诸元素对应乘积,只不过算L的元素时还要除以同列对角元.这一规律很容易记住.可写成算法(L和U可存放于A): for k=1:n-1 for j=k:n ukj=akj-lk1u1j-lk2u2j -…- lk,k-1uk-1,j end for i=k+1:n lik=(aik-li1u1k-li2u2k -…-li,k-1uk-1,k)/ukk end end 这一算法也叫Gauss消去法的紧凑格式,可一次算得L,U的元素,不需逐步计算存储. 考察上面的表格会发现还可安排其它计算次序,只要在这一次序下每个元素左边的L的元素与上方的U的元素已计算在先。例如, 逐行自左而右的次序, 逐列自上而下的次序, … 易知g的计算规律同U. 利用LU分解解Ax=b分三步: 1．分解A=LU 2．解 Lg=b 求g 3．解 Ux=y 求x 例3.????????????? 用直接LU分解法解 解 用分解公式计算得 求解 ③??? 其它分解 我们用顺序消元和直接分解两种方法实现了LU分解.还有更一般的三角分解,比如,下三角矩阵和单位上三角矩阵之积,又如单位下三角矩阵,对角矩阵,单位上三角矩阵之积,等等.下面给出第二种分解形式的算法LDR分解法。A=LDR,L是单位下三角矩阵,D是对角矩阵,R是单位上三角矩阵.逐列计算(逐列作LU分解,再用U的对角元素除各行),结果存入A。 for j=1:n for i=2:j aij=aij-ai1a1j-ai2a2j -…-ai,i-1ai-1,j end for i=j+1:n aij=(aij-ai1a1j-ai2a2j -…-ai,j-1aj-1,j)/ajj end for i=1:j-1 aij= aij/aii end end ④???? 列主元素的LU分解 对照顺序消元和LU分解,列主元素法也可得列主元素的LU分解: PA=LU P是行交换结果的排列阵,L和U同前. 例4.????????????? 列主元素法解方程组并写出系数矩阵的LU分解. 括号内是乘数,k=2时2,3行交换.因而有 直接作列主元素LU分解,因为在k步要先选主元素,所以作如下改变: for k=1:n-1 for i=k:n aik=aik-li1u1k-li2u2k -…-li,k-1uk-1,k end 找p: p行k行 ik=p for j=k+1:n ukj=akj-lk1u1j-lk2u2j -…- lk,k-1uk-1,j end for i=k+1:n lik= aik/ukk end end 可将lik存于aik,ukj存于akj. 二 实验部分 本章实验内容： 实验题目：