数值分析解线性方程组的直接方法演示文稿

数值分析解线性方程组的直接方法演示文稿目前一页\总数九十一页\编于二十一点优选数值分析解线性方程组的直接方法Ppt目前二页\总数九十一页\编于二十一点

在自然科学和工程技术中有很多问题的解决常常归结为解线性代数方程组．如三次样条函数问题，用最小二乘法求实验数据的曲线拟合问题，解非线性方程组问题，用差分法或者有限元方法解常微分方程、偏微分方程的边值问题等都导致求解线性代数方程组，而这些方程组的系数矩阵大致分为两种，一种是低阶稠密矩阵，另一种是大型稀疏矩阵。

关于线性方程组的数值解法一般有两类：1.直接法2.迭代法§5.1引言与预备知识5.1.1引言目前三页\总数九十一页\编于二十一点

本章讨论n元线性方程组

（5.1）的直接解法。方程组(5.1)的矩阵形式为

Ax=b其中

若矩阵A非奇异,即det(A)≠0,则方程组(2.1)有唯一解。目前四页\总数九十一页\编于二十一点

所谓直接解法是指，若不考虑计算过程中的舍入误差，经过有限次算术运算就能求出线性方程组的精确解的方法。但由于实际计算中舍入误差的存在，用直接解法一般也只能求出方程组的近似解。Cramer法则是一种不实用的直接法，本章将介绍几种实用的直接法。目前五页\总数九十一页\编于二十一点5.1.2预备知识M行n列矩阵.n维列向量.目前六页\总数九十一页\编于二十一点矩阵的基本运算:(1)矩阵的加法(7)矩阵的行列式行列式性质:(a)det(AB)=det(A)det(B)(6)非奇异矩阵(5)单位矩阵(4)转置矩阵(3)矩阵与矩阵的乘法(2)矩阵与标量的乘法目前七页\总数九十一页\编于二十一点矩阵特征值与谱半径定义1设若存在一个数λ(实数或复数)和非零向量使(1.1)则称λ为A的特征值,x为A对应λ的特征向量,A的全体特征值称为A的谱，记作称为A的谱半径.(1.2)目前八页\总数九十一页\编于二十一点由式(1.1)知,λ

可使齐次方程组有非零解,故系数行列式记称为矩阵Ａ的特征多项式，方程(1.3)称为Ａ的特征方程．(1.3)目前九页\总数九十一页\编于二十一点在复数域中有n个根故由行列式(1.3)展开可知：的n个特征值故是它的特征方程(1.3)的几个根,并有(1.4)(1.5)A的迹.目前十页\总数九十一页\编于二十一点A的特征值λ和特征向量x还有以下性质:(1)AT与A有相同的特征值λ及相同的特征向量x.(2)若A非奇异，则A-1的特征值为λ-1,特征向量为x.

(3)相似矩阵B=S-1AS有相同的特征多项式.目前十一页\总数九十一页\编于二十一点例1求的特征值及谱半径.解:A的特征方程为故A的特征值为A的谱半径为目前十二页\总数九十一页\编于二十一点5.1.4特殊矩阵目前十三页\总数九十一页\编于二十一点目前十四页\总数九十一页\编于二十一点定理1.目前十五页\总数九十一页\编于二十一点定理2.目前十六页\总数九十一页\编于二十一点定理3.定理4(Jordan标准型)设A为n阶矩阵,则存在一个非奇异矩阵P使得目前十七页\总数九十一页\编于二十一点其中:(1)当A的若当标准型中所有若当块Ji均为一阶时,此标准型变成对角矩阵.返回主页目前十八页\总数九十一页\编于二十一点求解高斯消去法(逐次消去法)及消去法和矩阵三角分解之间的关系：§5.2高斯消去法目前十九页\总数九十一页\编于二十一点例2

用消去法解方程组解第1步.将方程(2.2)乘上-2加到方程(2.4)上,消去未知数x1,得到(2.2)(2.3)(2.4)第2步.将方程(2.3)加到方程(2.5)上去,消去方程(2.5)中的x2,得到与原方程组等价的三角形方程组解为:首先举一个简单的例子来说明消去法的基本思想.目前二十页\总数九十一页\编于二十一点上述过程相当于思路首先将A化为上三角阵/*upper-triangularmatrix*/，再回代求解/*backwardsubstitution*/。=目前二十一页\总数九十一页\编于二十一点消元记Step1：设，计算因子将增广矩阵/*augmentedmatrix*/第i行mi1

第1行，得到其中Stepk：设，计算因子且计算目前二十二页\总数九十一页\编于二十一点回代注意1：只要A

非奇异，即A1

存在，则可通过逐次消元及行交换，将方程组化为三角形方程组，求出唯一解。共进行？步n

1目前二十三页\总数九十一页\编于二十一点注意2:设Ax=b,其中A为非奇异矩阵,如果由于A为非奇异矩阵,所以A的第一列一定有元素不等于零.例如目前二十四页\总数九十一页\编于二十一点定理5

设Ax=b,其中(1)如果则可通过高斯消去法将Ax=b约化为等价的三角形线性方程组(2.10),且计算公式为:①消元计算(k=1,2,…,n-1)②回代计算目前二十五页\总数九十一页\编于二十一点(2)如果A为非奇异矩阵,则可通过高斯消去法(及交换两行的初等变换)将方程组Ax=b约化为方程组(2.10).定理6

约化的主元素aii(i)

≠0(i=1,2,…,k)的充要条件是矩阵A的所有顺序主子式

/*determinantofleadingprincipalsubmatrices*/Di≠0(i=1,2,…,k).即(2.12)推论如果A的顺序主子式Dk≠0(k=1,2,…,n-1),则目前二十六页\总数九十一页\编于二十一点§5.2.2三角分解法

/*MatrixFactorization*/高斯消元法的矩阵形式

/*MatrixFormofG.E.*/：Step1:记L1=，则Stepn

1:其中

Lk=目前二十七页\总数九十一页\编于二十一点记为L单位下三角阵/*unitarylower-triangularmatrix*/记

U=目前二十八页\总数九十一页\编于二十一点定理7

若A的所有顺序主子式/*determinantofleadingprincipalsubmatrices*/

均不为0，则A

的

LU

分解唯一（其中L

为单位下三角阵）。证明：由§1中定理可知，LU分解存在。下面证明唯一性。若不唯一，则可设A=L1U1=L2U2

，推出Upper-triangularLower-triangularWithdiagonalentries1注：L

为一般下三角阵而U

为单位上三角阵的分解称为Crout分解。实际上只要考虑A*的LU

分解，即

，则即是A的Crout分解。目前二十九页\总数九十一页\编于二十一点例3

对于例2,系数矩阵由高斯消去法,返回主页目前三十页\总数九十一页\编于二十一点5.3.1列主元消去法:在顺序消元过程中,只要即可进行计算,但如果很小,则将导致舍入误差增长,使解的误差很大.例4用Gauss消去法求解方程组§5.3高斯主元素消去法目前三十一页\总数九十一页\编于二十一点解:因故方程有唯一解,且精确解为若用Gauss消去法取四位有效数字计算,可得解比较,误差很大,若将两个方程互换为用Gauss消去法取四位有效数字计算,可得解目前三十二页\总数九十一页\编于二十一点

本例表明通过行交换可避免舍入误差增长,这就是列主元消去法的基本思想.其计算步骤如下:第1步，在中的第1列选主元，即行为主元,若将

的第i1行与第1行互换，再按消元公式计算得到假定上述过程已进行(k-1)步，得到第k步，在中的第k列选主元，即若则在中将ik行与第k行互换，再按消元公式计算得到对k=1,2,…,n-1,重复以上过程则得如果某个k出现主元方程没有唯一解或严重病态，否则可由(3.2.4)求得解.则表明detA=0，目前三十三页\总数九十一页\编于二十一点它也表明当Ａ非奇异时，存在排列矩阵Ｐ(若干初等排列矩阵的乘积),使PA=LU,其中L为单位下三角矩阵,其元素|lij|<=1,U为上三角矩阵.上述每步行交换后再消元相当于其中是指标为k的初等下三角矩阵,为初等排列矩阵时,表示不换行,经过(n-1)步换行与消元,A化为上三角矩阵.即:目前三十四页\总数九十一页\编于二十一点解:例5用列主元消去法解Ａx=b,其中消元消元目前三十五页\总数九十一页\编于二十一点消元结束．由回代公式求得解此例的精确解为可见结果精度较高．若不选列主元Ｇauss消去法，求得解，误差较大．除列主元消去法外，还有一种消去法，是在Ａ的所有元素aij中选主元，称为全主元消去法．因计算量较大且应用列主元已能满足实际要求，故不再讨论．目前很多数学软件库都有列主元消去法，可直接调用．目前三十六页\总数九十一页\编于二十一点注:为了减少计算的舍入误差,使用消去法通常都要选主元.目前最常用的是列主元消去法,也就是每步消元之前选主元,当A=(aij)第一步选A中第1列的主元,即max|ai1|=ai1.然后将i1行与第１行互换，再进行消元，以后每步消元做法类似，先选主元，再消元．目前三十七页\总数九十一页\编于二十一点5.3.2高斯若当消去法消去对角线上方和下方的元素.假设已经完成k-1步,得到与方程Ax=b等价的方程组返回主页目前三十八页\总数九十一页\编于二十一点

高斯消去法有很多变形,有的是高斯消去法的改进、改写，有的是用于某一类特殊性质矩阵的高斯消去法的简化。5.3.1直接三角分解法.

将高斯消去法改写为紧凑形式,可以直接从A的元素得到计算L,U元素的递推公式,而不需要任何中间步骤,这就是所谓直接三角分解法.

一旦实现A的LU分解,那么求解Ax=b的问题就等价于求解两个三角形方程组

(1)Ly=b,求y;(2)Ux=y,求x.§５.4矩阵三角分解法

/*MatrixFactorization*/目前三十九页\总数九十一页\编于二十一点１.不选主元的三角分解法设A为非奇异矩阵,且有分解式A=LU,其中L为单位下三角,U为上三角即L,U元素可以由n步直接计算定出,其中第r步定出U的第r行和L的第r列元素.由上式有:目前四十页\总数九十一页\编于二十一点故同样有:

设已经定出U的第1行到第r-1行元素与L的第1列到第r-1列元素.利用矩阵乘法(注意当r<k时,lrk=0),有得目前四十一页\总数九十一页\编于二十一点通过比较法直接导出L和

U的计算公式。思路目前四十二页\总数九十一页\编于二十一点固定r:对i=r,r+1,…,n

有lii=1a将r

，i

对换，对r=i,i+1,…,n有b目前四十三页\总数九十一页\编于二十一点结论:用直接三角分解法解Ax=b(要求A的所有顺序主子式都不等于0)的计算公式如下.Step1:u1i=a1i;li1=ai1/u11;(i=2,…,n)计算U的第r行,L的第r列元素(r=2,3,…，n）Step2:求解Ly=b,Ux=y

的计算公式：Step3:Step4:Step5:目前四十四页\总数九十一页\编于二十一点例5

用直接三角分解法解解用分解公式计算得求解目前四十五页\总数九十一页\编于二十一点2.选主元的三角分解法当urr=0时计算中断，或者当urr绝对值很小时，按分解公式计算可能引起舍入误差的积累。但如果A非奇异，可以通过交换A的行实现矩阵A的LU分解，因此可采用与列主元消去法类似的方法，将直接三角分解法修改为（部分）选主元的三角分解法。设第r-1步分解已完成，这时有目前四十六页\总数九十一页\编于二十一点第r步分解需用到（3.2）及（3.3）式，为避免用小的数urr做除数，引进量于是有：取交换A的r行与ir行元素，将调到(r,r)位置(将（i，j）位置的新元素仍记为lij

及

aii),于是有|lir|<=1(i=r+1,…，n）.由此再进行第r步分解计算。目前四十七页\总数九十一页\编于二十一点5.3.2平方根法当A对称正定时，A的顺序主子式故由定理知，A=LU的分解存在且唯一，其中L为单位下三角为了A利用对称性其中D为对角阵,U0为单位上三角阵,于是又代入到上式,就得到对称矩阵A的分解式矩阵，U为上三角矩阵，且目前四十八页\总数九十一页\编于二十一点定理9(对称阵的三角分解定理)设A为n阶对称阵,且A的顺序主子式则A可唯一分解为，其中L为单位下三角矩阵，.

D为对角矩阵定理10(对称正定矩阵的三角分解或Cholesky分解)

如果A为n阶对称正定矩阵,则存在一个实的非奇异下三角阵L使A=LLT,当限定L的对角元素为正时,这种分解是唯一的.将求方程组的解转化为求方程的解LLTx=b.令

，求得方程的解由根据矩阵乘法，由目前四十九页\总数九十一页\编于二十一点得

i=j有

当i＞j,得目前五十页\总数九十一页\编于二十一点例6用平方根法求以下方程组的解.

解先验证系数矩阵A对称正定，对称显然，

故A对称正定，可用Cholesky分解计算，求得

求解

得再求得目前五十一页\总数九十一页\编于二十一点5.3.3追赶法解三对角方程组

/*CroutReductionforTridiagonalLinearSystem*/在一些实际问题中,例如解常微分方程边值问题,解热传导方程以及船体数学放样中建立三次样条函数等,都会要求解系数矩阵为对角占优的三对角线方程组.简记为Ax=f.其中,当|i-j|>1时,aij=0,且：目前五十二页\总数九十一页\编于二十一点Step1:对A作Crout分解直接比较等式两边的元素，可得到计算公式。Ｌ为下三角，Ｕ为单位上三角目前五十三页\总数九十一页\编于二十一点注意当j=1时有对j=2，3，…，n求得L的元素，这就是A的Cholesky分解，然后再解两个三角方程组，得这就是对称正定方程组的平方根法

另外，由于

故有这表明分解过程中矩阵L中元素因此平方根法计算是数值稳定的.

的数量级不增长，目前五十四页\总数九十一页\编于二十一点Step2:追——即解Step3:赶——即解：求解Ａx＝f等价于Step１:计算的递推公式将计算系数为追的过程将计算方程组的解为赶的过程目前五十五页\总数九十一页\编于二十一点定理：设有三对角线方程组Ａx＝f，其中Ａ满足对角占优的条件，则Ａ为非奇异矩阵且追赶法计算公式中的满足：目前五十六页\总数九十一页\编于二十一点定理

若A

为对角占优

/*diagonallydominant*/的三对角阵，且满足，则追赶法可解以A

为系数矩阵的方程组。注：

如果A是严格对角占优阵，则不要求三对角线上的所有元素非零。

根据不等式可知：分解过程中，矩阵元素不会过分增大，算法保证稳定。

运算量为O(6n)。返回主页目前五十七页\总数九十一页\编于二十一点5.5.1内积与向量范数为了研究方程组Ax=b解的误差和迭代法收敛性，需对向量及矩阵的"大小"引进一种度量，就要定义范数，它是向量"长度"概念的直接推广，通常用表示n维实向量空间，表示n维复向量空间.定义2

设将实数称为向量x,y的数量积.非负实数称为向量x的欧氏范数或2-范数.§5.5向量和矩阵范数目前五十八页\总数九十一页\编于二十一点定理12设则内积有以下性质：(1)(2)(3)(4)(5)(柯西-施瓦茨不等式)等式当且仅当x与y线性相关时成立;(6)三角不等式目前五十九页\总数九十一页\编于二十一点定义3(向量范数)如果向量的某个实值函数满足条件:则称目前六十页\总数九十一页\编于二十一点对于由内积性质可知它满足定义2的三个条件，故它是一种向量范数.此外还有以下几种常用的向量范数.容易验证均满足定义2的三个条件.更一般的还可定义但只有p=1,2,∞时的三种范数是常用的向量范数.例如给定则可求出目前六十一页\总数九十一页\编于二十一点定理14设是则N(x)是向量x的分量上任一种向量范数，的连续函数.定理15(向量范数的等价性)设是上任意两种向量范数，则存在常数使目前六十二页\总数九十一页\编于二十一点5.5.2矩阵范数

矩阵可看成n×n维向量，如果直接将向量的2-范数用于矩阵A，则可定义称为矩阵A的Frobenius范数，简称F-范数.它显然满足向量范数的三条性质，但由于矩阵还有乘法运算，因此矩阵范数的定义中应增加新条件.目前六十三页\总数九十一页\编于二十一点定义4

如果的某个非负实函数N(A)，记作‖A‖，满足条件：则称目前六十四页\总数九十一页\编于二十一点显然满足定义中的四个条件,(3)，(4)两条均可由Cauchy-Schwarz不等式证明，故是一种矩阵范数.

除矩阵自身的运算外，在解方程中矩阵乘向量的运算即Ax，也是必不可少的.因此要求所引进的范数应满足条件：上式称为相容性条件.为使引进的矩阵范数满足条件(4.5)，我们给出以下定义.(4.5)目前六十五页\总数九十一页\编于二十一点定义6(矩阵的算子范数)设当给定向量范数时可定义称为矩阵的算子范数或从属范数.(4.6)定理17

设上的一种向量范数，则由(4.6)定义的是一种矩阵范数，且满足相容性条件目前六十六页\总数九十一页\编于二十一点证明因中有界闭集上的连续函数，故在D上有最大值，即使而对故所以从而当成立，而x=0时显然也成立.目前六十七页\总数九十一页\编于二十一点定理17

设则这里为矩阵的谱半径.目前六十八页\总数九十一页\编于二十一点例7

已知解从定理可以看出，计算较容易，而计算

时因为要求的特征值，所以较为困难.但当A对称时，有目前六十九页\总数九十一页\编于二十一点定理19定理18

对任何为任一种从属范数则反之，对任意ε＞0，至少存在一种从属范数使证明:设为A的特征值，则由得目前七十页\总数九十一页\编于二十一点非奇异，且证明用反证法.假定(I+B)奇异，则齐次方程有非零解而与‖B‖＜1的假设矛盾，故(I+B)非奇异.

又得取范数得定理20设

返回主页目前七十一页\总数九十一页\编于二十一点矩阵条件数与扰动方程组误差界

在解方程组Ax=b时，由于各种原因，A或b往往有误差，从而使得解也产生误差.例8

方程组

的准确解为

,当A与b有微小变化时，如变为方程则准确解为

它表明A,b的微小扰动引起方程解x的很大变化，这就是病态方程.§5.6误差分析与病态方程组目前七十二页\总数九十一页\编于二十一点定义7

求解线性方程组Ax=b时，若A或b有微小扰动

解x的误差

很大,，则称此方程组为病态方程组，相应的系数矩阵A称为病态矩阵.

反之，若此时

很小,，则称此方程组为良态方程组，相应的系数矩阵A称为良态矩阵.

注意方程组是否病态与用什么数值方法无关，它是由方程自身性质决定的.

在例8中因为行列式

因此出现病态.但有时A从表面上看性质很好，也可能是病态的.目前七十三页\总数九十一页\编于二十一点那么如何判断A是否病态?先给出如下定义.例9

方程组Ax=b表示为它的准确解

A对称正定且

表面看性质"较好"，但若对右端b作微小变化，如方程改为

则解变为

这里b的相对误差大约只有

但解的相对误差却很大，故A也是病态矩阵.目前七十四页\总数九十一页\编于二十一点定义8

设

非奇异，‖·‖v为矩阵的任一种从属范数，则

称为矩阵A的条件数.

从定义看到矩阵条件数依赖于范数的选取，如范数为2-范数，

则记为

同理有

等等.目前七十五页\总数九十一页\编于二十一点条件数有以下性质：

(1)(2)(3)U为正交矩阵，则

(4)若

与为A的按模最大与最小特征值，则若A对称，则

目前七十六页\总数九十一页\编于二十一点下面给出扰动方程组解的误差分析.先考察b有扰动

则扰动方程为由于Ax=b,故得

于是再由Ax=b,有

即故得目前七十七页\总数九十一页\编于二十一点下面再研究方程Ax=b，当A有扰动

时，其解

的误差分析.

此时扰动方程为

因Ax=b,故有

因存在,若假定

则由定理20可知

非奇异,并有:(5.6)由(5.6)可得

因此(5.7)目前七十八页\总数九十一页\编于二十一点定理22

设A为非奇异矩阵,Ax=b≠0,且如果则(5.7)式成立.从(5.7)看到，当A的条件数Cond(A)很大时，解的相对误差

也很大，故方程组为病态.在例9中

而于是条件数很大，故方程是严重病态的.

目前七十九页\总数九十一页\编于二十一点例10

Hibert矩阵是一个著名的病态矩阵，记作

它是一个对称正定矩阵，当n≥3时它是病态矩阵.例如

故另外还有

等等.因此Hn是严重的病态矩阵，且n越大

Cond(Hn)越大.目前八十页\总数九十一页\编于二十一点例11

在例10的方程组中可算出A的特征值

故例中实际相对误差是而根据(5.6)的误差估计为这与实际相差不大，即相对误差放大了将近3000倍.故方程为病态方程组.目前八十一页\总数九十一页\编于二十一点定理23(事后误差估计)设方程组

，则若实际求得解为证明记剩余则它表明如果方程组病态，即使剩余‖r‖很小，解的相对误差仍可能很大.目前八十二页\总数九十一页\编于二十一点5.6.2病态方程组的解法

如果A的条件数Cond(A)>>1，则Ax=b为病态方程，但计算Cond(A)时需要求A-1,计算量很大，相当于解方程组，在实际中常可通过求解过程直观地判断方程组的病态性质，如果解方程时出现下述情况之一，则可能是"病态"方程组.(1)在列主元消去法中出现小主元；(2)在计算过程中行或列几乎线性相关或三角分解中对角元出现近似零的元素；(3)矩阵A的元素数量级相差很大且无规律；(4)剩余很小，而解很大，又达不到精度要求.目前八十三页\总数九十一页\编于二十一点

(1)采用高精度运算，减轻病态影响，例如用双倍字长运算.对病态方程组求解可采用以下措施：(2)用预处理方法改善A的条件数，即选择非奇矩阵与Ax=b等价,而(3)平衡方法，当A中元素的数量级相差很大，可采用行均衡或列均衡的方法改善A的条件数.设非奇异，计算于是求Ax=b等价于求的条件数可得到改善，这就是行均衡法.目前八十四页\总数九十一页\编于二十一点例12

给定方程组Ax=b为A的条件数若用行均衡法可取则平衡后的方程用三位有效数字的列主元消去法求解得目前八十五页\总数九十一页\编于二十一点functionx=threedia(a,b,c,f)N=length(f);x=zeros(1,N);y=zeros(1,N);beta=zeros(1,N);gramma=zeros(1,N);beta(1)=b(1);fori=1:N-1gramma(i)=c(i)/beta(i);beta(i+1)=b(i+1)-a(i+1)*gramma(i);end%追的过程y(1)=f(1)/beta(1);fori=2:Ny(i)=(f(i)-a(i)*y(i-1))/beta(i);end%赶的过程x(N)=y(N);fori=N-1:-1:1x(i)=y(i)-gramma(i)*x(i+1);end

a=[0,-1,-1,-3];>>b=[2,3,2,5];>>c=[-1,-2,-1,0];>>f=[6,1,0,1]';>>x=threedia(a,b,c,f)追赶法求解三对角方程组目前八十六页\总数九十一页\编于二十一点Cholesky方法：

A=[4,-2,4;-2,17,10;4,10,9];

b=[8.7,13.7,-0.7]';

[x,L,D]=Chol_decompose(A,b)L=1.000000-0.50001.000001.00000.75001.0000D=416-4x=-5.1457-3.17275.7344%用Cholesky分解求解%A是对称矩阵%L是单位下三角阵%D是对角阵%对称阵A进行三角分解：A=LDL'目前八十七页\总数九十一页\编于二十一点function[x,L,D]=Chol_decompose(A,b)N=length(A);L=zeros(N,N);D=zeros(1,N);fori=1:NL(i,i)=1;endD(1)=A(1,1);fori=2:Nforj=1:i-1ifj==1L(i,j)=A(i,j)/D(j);elsesum1=0;fork=1:j-1sum1=sum1+L(i,k)*D(k)*L(j,k);endL(i,j)=(A(i,j)-sum1)/D(j);endendsum2=0;fork=1:i-1sum2=sum2+L(i,k)^2*D(k);endD(i)=A(i,i)-sum2;end%分别求解线性方程组Ly=b;L'x=y/Dy=zeros(1,N);y(1)=b(1);fori=2:Nsumi=0;fork=1:i-1sumi=sumi+L(i,k)*y(k);endy(i)=b(i)-sumi;endx=zeros(1,N);x(N)=y(N)/D(N);fori=N-1:-1:1sumi=0;fork=i+1:Nsumi=sumi+L(k,i)*x(k);endx(i)=y(i)/D(i)-sumi;end目前八十八页\总数九十一页\编于二十一点用Dollittle三角分解法求解方程组：>>A=[0.001,2,3;-1,3.712,4.623;-2,1.072,5.643];>>b=[1,2,3]';>>[x,L,U]=lu_decompose(A,b)x=-0.4904-0.05100.3675L=1.0e+003*0.001000-1.00000.00100-2.00000.00200.0010U=1.0e+003*0.00000.00200.003002.00373.0046000.0059目前八十九页\总数九十一页\编于二十一点function[x,L,U]=lu_decompose(A,b)%用Dollittle三角分解%