线性代数方程组的直接法 (2)幻灯片.ppt

线性代数方程组的直接法第1页，共54页，编辑于2022年，星期一问题驱动：投入产出分析问题驱动：投入产出分析 投入产出分析是投入产出分析是2020世纪世纪3030年代由美国经济学家首先提出的，年代由美国经济学家首先提出的，它是研究整个经济系统各部门之间它是研究整个经济系统各部门之间“投入投入”与与“产出产出”关系的线性关系的线性模型，一般称为投入产出模型。国民经济各个部门之间存在着模型，一般称为投入产出模型。国民经济各个部门之间存在着相互依存的关系，每个部门在运转中将其它部门的成品或半成相互依存的关系，每个部门在运转中将其它部门的成品或半成品经过加工（称为投入）变为自己的产品（称为产出），如何品经过加工（称为投入）变为自己的产品（称为产出），如何根据各部门之间的投入根据各部门之间的投入-产出关系，确定各部门的产出水平，以产出关系，确定各部门的产出水平，以满足社会的需求，是投入产出综合平衡模型研究的问题，试讨论满足社会的需求，是投入产出综合平衡模型研究的问题，试讨论如下简化问题。如下简化问题。设国民经济仅由农业、制造业和服务业三个部门构成，已知设国民经济仅由农业、制造业和服务业三个部门构成，已知某年它们之间的投入和产出关系、外部需求、初始投某年它们之间的投入和产出关系、外部需求、初始投入等如表入等如表6.1.16.1.1所示（数字表示产值，单位为亿元）所示（数字表示产值，单位为亿元）。第2页，共54页，编辑于2022年，星期一表表6.1.1 6.1.1 国民经济各个部门间的关系国民经济各个部门间的关系表中第一行数字表示农业总产出为表中第一行数字表示农业总产出为100100亿元，其中亿元，其中1515亿元农产品用亿元农产品用于农业生产本身，于农业生产本身，2020亿元用于制造业，亿元用于制造业，3030亿元用于服务业，剩下的亿元用于服务业，剩下的3535亿元农产品用于满足外部需求。类似地可以解释第二、三行数字。亿元农产品用于满足外部需求。类似地可以解释第二、三行数字。第一列数字中，第一列数字中，1515亿元如前所述，亿元如前所述，3030亿元是制造业对农业的投入，亿元是制造业对农业的投入，2020亿元是服务业对农业的投入，亿元是服务业对农业的投入，3535亿元的初始投入包括工资、税收、亿元的初始投入包括工资、税收、进口等，总投入进口等，总投入100100亿元和总产出相等。假定每个部门的产出和投亿元和总产出相等。假定每个部门的产出和投入是成正比的，由表入是成正比的，由表6.1.16.1.1能够确定这三个部门的投入产出表，如能够确定这三个部门的投入产出表，如表表6.1.26.1.2所示。所示。第3页，共54页，编辑于2022年，星期一表表6.1.2 6.1.2 投入产出表投入产出表表中的第一行，第二列的数字表示生产表中的第一行，第二列的数字表示生产1 1个单位产值的制造业个单位产值的制造业产品需要投入产品需要投入0.100.10个单位的产值的农产品，同样第三行、第个单位的产值的农产品，同样第三行、第一列的数字表示，生产一列的数字表示，生产1 1个单位产值的农产品需要个单位产值的农产品需要0.200.20个单位个单位的服务业产值。表的服务业产值。表6.1.26.1.2的数字称为投入系数和消耗系数，如的数字称为投入系数和消耗系数，如果技术水平没有变化，可以假设投入系数是常数。已知投入果技术水平没有变化，可以假设投入系数是常数。已知投入系数如表系数如表2.1.22.1.2所示，若今年对农业、制造业和服务业的外部所示，若今年对农业、制造业和服务业的外部需求分别为需求分别为5050、150150、100100亿元，试计算三个部门的总产出分亿元，试计算三个部门的总产出分别为多少？别为多少？第4页，共54页，编辑于2022年，星期一若共有若共有n个部门，记一定时期内第个部门，记一定时期内第i个部门的总产出为个部门的总产出为 xi,其中对第其中对第 j个部门的投入为个部门的投入为xij，满足的外部需求为，满足的外部需求为 di，则，则 (6.1.1)记第记第j个部门的单位产出需要第个部门的单位产出需要第i个部门的投入为个部门的投入为aij，在每个部门的产出，在每个部门的产出与投入成正比的假定下，有与投入成正比的假定下，有(6.1.2)投入系数即为投入系数即为aij，将，将(6.1.2)式代入式代入(6.1.1)）式得方程组）式得方程组 用矩阵表示为用矩阵表示为或(6.1.3)因此投入产出模型最终可归结为求解线性方程组的问题，下面因此投入产出模型最终可归结为求解线性方程组的问题，下面介绍求解线性方程组数值方法。介绍求解线性方程组数值方法。第5页，共54页，编辑于2022年，星期一AXAX=b b（3.1）第6页，共54页，编辑于2022年，星期一线性方程组数值解法的分类线性方程组数值解法的分类 直接法直接法（适用于中等规模的（适用于中等规模的n n阶线性方程组）阶线性方程组）GaussGauss消去法及其变形消去法及其变形 矩阵的三角分解法矩阵的三角分解法 迭代法迭代法（适用于高阶线性方程组）（适用于高阶线性方程组）JacobiJacobi迭代法迭代法 Gauss-SeidelGauss-Seidel迭代法迭代法 逐次超松弛法逐次超松弛法 共轭斜量法共轭斜量法第7页，共54页，编辑于2022年，星期一1 1 高斯消去法高斯消去法1 1三角形方程组的解法三角形方程组的解法-回代法回代法（3.2）（3.3）第8页，共54页，编辑于2022年，星期一2 2顺序高斯消去法顺序高斯消去法 基本思想：基本思想：通过消元将上述方程组化为上三角形通过消元将上述方程组化为上三角形 方程组，再利用回代法进行求解。方程组，再利用回代法进行求解。第9页，共54页，编辑于2022年，星期一记记第10页，共54页，编辑于2022年，星期一消元公式消元公式消元公式消元公式回代公式回代公式回代公式回代公式对对k=1,2,k=1,2,n-1,n-1,若若则可得则可得第11页，共54页，编辑于2022年，星期一顺序顺序GaussGauss消去法可执行的前提消去法可执行的前提定理定理 1 1 给定线性方程组给定线性方程组 ，如果，如果n n阶方阵阶方阵 的所有顺序主子式都不为零，即的所有顺序主子式都不为零，即 则按顺序则按顺序GaussGauss消去法所形成的各主元素消去法所形成的各主元素 均不为零，从而顺序均不为零，从而顺序Gauss 消去法可顺利执行。消去法可顺利执行。注：注：当线性方程组的系数矩阵为当线性方程组的系数矩阵为对称正定对称正定或或严格对角占优阵严格对角占优阵时，时，按顺序按顺序GaussGauss消去法计算是稳定的。消去法计算是稳定的。第12页，共54页，编辑于2022年，星期一3、列主元、列主元Gauss消去法计算步骤：消去法计算步骤：1、输入矩阵阶数输入矩阵阶数n，增广矩阵增广矩阵 A(n,n+1);2、对于对于(1)按列选主元：选取按列选主元：选取 l 使使(2)如果如果 ，交换，交换 A(n,n+1)的第的第k行与第行与第l 行元素行元素(3)消元计算消元计算 ：3、回代计算回代计算第13页，共54页，编辑于2022年，星期一4 4无回代过程的主元消去法无回代过程的主元消去法(Gauss-Jordan(Gauss-Jordan消去法消去法)第一步：选主元，在第一列中选绝对值最大的元素，设第第一步：选主元，在第一列中选绝对值最大的元素，设第k行为行为主元行，将主元行换至第一行，将第一个方程中主元行，将主元行换至第一行，将第一个方程中x1的系数变为的系数变为1，并从其余并从其余n 1n 1个方程中消去个方程中消去x1。第二步：在第二列后第二步：在第二列后n 1个元素中选主元，将第二个方程中个元素中选主元，将第二个方程中 x2的系数变为的系数变为1，并从其它，并从其它n n 1 1个方程中消去个方程中消去x2。第第k步：在第步：在第k列后列后n k个元素中选主元，换行，将第个元素中选主元，换行，将第k个方程个方程 xk的系数变为的系数变为1，从其它，从其它n-1n-1个方程中消去变量个方程中消去变量xk第14页，共54页，编辑于2022年，星期一消元公式为：消元公式为：对对k=1,2,按上述步骤进行到第按上述步骤进行到第n步后，方程组变为：步后，方程组变为：即为所求的解即为所求的解.第15页，共54页，编辑于2022年，星期一注：注：无回代的无回代的GaussGauss消元法实际上就是将消元法实际上就是将 方方程组的系数矩阵化为行最简形矩阵。程组的系数矩阵化为行最简形矩阵。第16页，共54页，编辑于2022年，星期一5无回代消去法的应用无回代消去法的应用（1）解线性方程组系解线性方程组系设要解的线性方程组系为：设要解的线性方程组系为：AX=b1,AX=b2,AX=bm上述方程组系可以写为上述方程组系可以写为AX=B=(b1,bm)第17页，共54页，编辑于2022年，星期一因此因此X=A-1B 即为线性方程组系的解。即为线性方程组系的解。在计算机上只需要增加几组右端常数项的存贮单元，在计算机上只需要增加几组右端常数项的存贮单元，其结构和解一个方程组时一样。其结构和解一个方程组时一样。行行系数系数右端右端第18页，共54页，编辑于2022年，星期一（2）求逆矩阵求逆矩阵设设A=(aij)n n是非奇矩阵是非奇矩阵，A 0，且令且令由于由于 AA-1=AX=I因此，求因此，求A-1的问题相当于解下列线性方程组的问题相当于解下列线性方程组相当于相当于（1）中中m=n,B=I 的情形。的情形。第19页，共54页，编辑于2022年，星期一（3）求行列式的值求行列式的值用高斯消去法将用高斯消去法将 A 化成上三角形化成上三角形第20页，共54页，编辑于2022年，星期一例例 用用Gauss-Jordan消去法解方程组消去法解方程组，并求出并求出 其中其中 解：把系数矩阵、单位矩阵和右端项组成增广矩阵，对增广矩解：把系数矩阵、单位矩阵和右端项组成增广矩阵，对增广矩 阵实行阵实行Gauss-Jordan消元过程消元过程。第21页，共54页，编辑于2022年，星期一故故，系数矩阵的逆为系数矩阵的逆为 第22页，共54页，编辑于2022年，星期一2 2 解三对角方程组的追赶法解三对角方程组的追赶法第23页，共54页，编辑于2022年，星期一第24页，共54页，编辑于2022年，星期一33 矩阵的矩阵的三角分解法三角分解法 高斯消元法的矩阵形式高斯消元法的矩阵形式 每一步消去过程相当于左乘初等下三角矩阵每一步消去过程相当于左乘初等下三角矩阵L Lk k记记其中，其中，第25页，共54页，编辑于2022年，星期一记记其中，其中，一般的，对一般的，对第26页，共54页，编辑于2022年，星期一A 的的 LU 分解分解(LU factorization)记记则则其中，其中，于是于是第27页，共54页，编辑于2022年，星期一 矩阵三角分解的定义矩阵三角分解的定义定义定义1 1设设A A为为n n n n实矩阵，如果存在下三角矩阵实矩阵，如果存在下三角矩阵L L与与上三角矩阵上三角矩阵U U，使得，使得 A=LU A=LU，则称之为矩阵，则称之为矩阵A A的的三角分解；三角分解；若存在单位下三角矩阵若存在单位下三角矩阵L L，对角，对角矩阵矩阵D D及单位上三角矩阵及单位上三角矩阵R R，使得，使得A=LDR,A=LDR,则称之则称之为矩阵为矩阵A A的的LDRLDR分解分解。第28页，共54页，编辑于2022年，星期一注注注注：（1 1）L 为为单位单位下三角阵而下三角阵而 U 为为一般一般上三角上三角阵的分解称为阵的分解称为Doolittle 分解；分解；（2）L 为为一般一般下三角阵而下三角阵而 U 为为单位单位上三角阵上三角阵的分解称为的分解称为Crout 分解分解。第29页，共54页，编辑于2022年，星期一定理定理2：设设A为为n n实矩阵，如果求解实矩阵，如果求解AX=b用顺序高斯消用顺序高斯消去法能够完成（即去法能够完成（即 ），则矩阵），则矩阵A可唯一分解为一个单位下三角矩阵可唯一分解为一个单位下三角矩阵L与上三角矩阵与上三角矩阵U的乘积。的乘积。A=LU 矩阵三角分解的存在性矩阵三角分解的存在性第30页，共54页，编辑于2022年，星期一通过比较法直接导出通过比较法直接导出L 和和 U 的计算公式。的计算公式。思思路路 矩阵三角分解的实现矩阵三角分解的实现DoolittleDoolittle三角三角分解：分解：第31页，共54页，编辑于2022年，星期一求解线性方程组的（求解线性方程组的（Doolittle）三角分解法）三角分解法第32页，共54页，编辑于2022年，星期一直接三角分解法解直接三角分解法解AX=b的计算公式的计算公式对于对于r=2,3,n（2）计算计算U的第的第r行元素行元素（3）计算计算L的第的第r 列元素列元素(r n)(1)第33页，共54页，编辑于2022年，星期一(4)(5)第34页，共54页，编辑于2022年，星期一例例 用矩阵的用矩阵的Doolittle Doolittle 三角分解法解方程组三角分解法解方程组第35页，共54页，编辑于2022年，星期一 Doolittle分解法的变形分解法的变形 紧凑格式的紧凑格式的DoolittleDoolittle三角分解法三角分解法例例所以所以第36页，共54页，编辑于2022年，星期一紧凑格式的紧凑格式的列主元列主元DoolittleDoolittle分解法分解法例例第37页，共54页，编辑于2022年，星期一4 4 平方根法平方根法1 1矩阵的矩阵的LDRLDR分解分解定理定理3：如果如果n阶矩阵阶矩阵A的所有顺序主子式均不等于零，的所有顺序主子式均不等于零，则矩阵则矩阵A存在唯一的分解式存在唯一的分解式A=LDR，其中其中L和和R分别是分别是n阶单位下三角阵和单位上三角阵，阶单位下三角阵和单位上三角阵，D是对角元素不为零的是对角元素不为零的n阶对角阵，上述分解称为阶对角阵，上述分解称为A的的LDR分解分解。第38页，共54页，编辑于2022年，星期一 2平方根法平方根法 如果如果A为对称正定矩阵为对称正定矩阵，则存在一个实的非奇异下三则存在一个实的非奇异下三角矩阵角矩阵，使使A=LLT ，且当限定的对角元素为正时且当限定的对角元素为正时，这种分解是唯一的，称为矩阵这种分解是唯一的，称为矩阵A的的cholesky分解分解。定理定理4：（对称正定矩阵的三角分解）（对称正定矩阵的三角分解）第39页，共54页，编辑于2022年，星期一将将对称对称 正定阵正定阵 A 做做 LU 分解分解U=uij=u11uij/uii111u22unn记为记为 A 对称对称即即记记 D1/2=则则 仍是下三角阵，且有仍是下三角阵，且有对称正定阵对称正定阵choleskycholesky分解的实现分解的实现第40页，共54页，编辑于2022年，星期一用平方根法解对称正定线性代数方程组的算法用平方根法解对称正定线性代数方程组的算法（1）对矩阵对矩阵A进行进行Cholesky分解分解，即即A=LLT，由矩阵乘法由矩阵乘法：对于对于 i=1,2,n 计算计算第41页，共54页，编辑于2022年，星期一（2）求解下三角形方程组求解下三角形方程组 （3）求解求解LTX=y第42页，共54页，编辑于2022年，星期一3改进平方根法改进平方根法 其中其中第43页，共54页，编辑于2022年，星期一改进平方根法解对称正定方程组的算法改进平方根法解对称正定方程组的算法 第44页，共54页，编辑于2022年，星期一令令LTX=y，先解下三角形方程组先解下三角形方程组 LDY=b 得得解上三角形方程组解上三角形方程组 LTX=Y 得得 第45页，共54页，编辑于2022年，星期一求解求解 时，时，A 和和 的误差对解的误差对解 有何影响？有何影响？设设 A 精确，精确，有误差有误差 ，得到的解为，得到的解为 ，即，即绝对误差放大因子绝对误差放大因子又又相对误差放大因子相对误差放大因子5 线性方程组的性态和解的误差分析线性方程组的性态和解的误差分析第46页，共54页，编辑于2022年，星期一第47页，共54页，编辑于2022年，星期一6 Error Analysis for .设设 精确，精确，A有误差有误差 ，得到的解为，得到的解为 ，即，即 Wait a minute Who said that(I+A 1 A)is invertible?(只要只要 A充分小，使得充分小，使得 是关键是关键的误差放大因子，称为的误差放大因子，称为A的的条件数条件数，记为，记为cond(A),越越 则则 A 越病态，越病态，难得准确解。难得准确解。大大第48页，共54页，编辑于2022年，星期一例：例：Hilbert 阵阵cond(H2)=27cond(H3)748cond(H6)=2.9 106cond(Hn)as n 注：注：现在用现在用MatlabMatlab数学软件可以很方便求矩阵的状态数数学软件可以很方便求矩阵的状态数!定义定义2:2:设线性方程组的系数矩阵是非奇异的设线性方程组的系数矩阵是非奇异的,如果如果condcond(A A)越大越大,就称这个方程组就称这个方程组越病态越病态.反之反之,condcond(A A)越小越小,就称这个方程组就称这个方程组越良态越良态.第49页，共54页，编辑于2022年，星期一一般判断矩阵是否病态，并不计算一般判断矩阵是否病态，并不计算A A 1 1，而由经验得出。，而由经验得出。行列式很大或很小（如某些行、列近似相关）；行列式很大或很小（如某些行、列近似相关）；元素间相差大数量级，且无规则；元素间相差大数量级，且无规则；主元消去过程中出现小主元；主元消去过程中出现小主元；特征值相差大数量级。特征值相差大数量级。第50页，共54页，编辑于2022年，星期一 近似解的误差估计及改善：近似解的误差估计及改善：设设 的近似解为的近似解为，则一般有，则一般有cond(A)误差上限误差上限 改善方法改善方法(1)：Step 1:近似解近似解Step 2:Step 3:Step 4:若若 可被精确解出，则有可被精确解出，则有 就是精确解了。就是精确解了。经验表明经验表明：若：若 A 不是非常病态（例如：不是非常病态（例如：），则如），则如此迭代可达到机器精度；但若此迭代可达到机器精度；但若 A 病态，则此算法也不能改病态，则此算法也不能改进。进。第51页，共54页，编辑于2022年，星期一 改善方法改善方法(2)(2)对方程组进行预处理，即适当选择非奇异对角阵对方程组进行预处理，即适当选择非奇异对角阵D D，C,C,使求解使求解 Ax=b Ax=b 的问题转化为求解等价方程组的问题转化为求解等价方程组 DACCDACC-1-1x=Db,x=Db,且使且使DACDAC 的条件数得到改善。（的条件数得到改善。（P88P88，例，例3.103.10）用双精度进行计算，以便改善和减轻病态矩阵的影响用双精度进行计算，以便改善和减轻病态矩阵的影响。第52页，共54页，编辑于2022年，星期一历史与注记历史与注记 高斯高斯(Carl Friedrich Gauss,1777-1855)高斯是德国高斯是德国数学家、物理学家、天文学家。数学家、物理学家、天文学家。1777年生于德国布伦瑞年生于德国布伦瑞克，克，1855年在哥廷根逝世。高斯是近代数学奠基者之一，年在哥廷根逝世。高斯是近代数学奠基者之一，在历史上影响之大，在历史上影响之大，可以和阿基米德、牛顿、欧拉并列，可以和阿基米德、牛顿、欧拉并列，有有“数学王子数学王子”之称。之称。1795 年进入格丁根大学学习。年进入格丁根大学学习。1798 年转入黑尔姆施泰特大学，年转入黑尔姆施泰特大学，1799 年获博士学位。年获博士学位。1807 1807 年以后一直在格丁根大学任教授和哥廷根天文台台长，一直到逝世。年以后一直在格丁根大学任教授和哥廷根天文台台长，一直到逝世。1833年和物理学家年和物理学家W.E.韦伯共同建立地磁观测台，组织磁学学会以联系全韦伯共同建立地磁观测台，组织磁学学会以联系全世界的地磁台站网。世界的地磁台站网。高斯的数学研究几乎遍及所有领域，在数论、代数学、非欧几何、复高斯的数学研究几乎遍及所有领域，在数论、代数学、非欧几何、复变函数和微分几何等方面都做出了开创性的贡献。高斯还把数学应用于天变函数和微分几何等方面都做出了开创性的贡献。高斯还把数学应用于天文学、大地测量学和磁学的研究，他在数值计算方面的贡献主要有：提出文学、大地测量学和磁学的研究，他在数值计算方面的贡献主要有：提出了高斯消去法和高斯迭代法、高斯积分，发明了最小二乘法原理等。了高斯消去法和高斯迭代法、高斯积分，发明了最小二乘法原理等。1810年，高斯提出了消去法，矩阵是凯莱在年，高斯提出了消去法，矩阵是凯莱在1855年提出的，而将高斯消年提出的，而将高斯消去表示成矩阵分解是在去表示成矩阵分解是在20世纪世纪40年代由年代由Dwyer、冯、冯.诺伊曼等人提出的。诺伊曼等人提出的。第53页，共54页，编辑于2022年，星期一参考文献参考文献 1 J.H.Wilkinson.Error analysis of direct methods of matrix inversion.J.ACM,8:281-330,1961.2 J.H.Wilkinson.Rounding Errors in Algebraic Processs.Prentice Hall,Englewood Cliffs,NJ,1963.3 J.H.Wilkinson.The Algebraic Eigenvalue Problem.Oxford University Press,New York,1965.4 G.E.Forsythe and C.B.Moler.Computer Solution of Linear Algebraic Systems.Prentice Hall,Englewood Cliffs,NJ,1967.第54页，共54页，编辑于2022年，星期一