111-5-3平方根法、追赶法

合集下载

追赶法

追赶法

追赶法/平方根法 例2.4.1 设4阶方程组AX=B 为⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡-=⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡------121653342231124321x x x x这就是一个三对角方程组,既系数矩阵除了对角线的“三斜线”以外的元素均为0。

用追赶法求解三对角方程组的一种做法是把系数矩阵A 写成下列形式的LU 分解(这里采用Doolittle 分解,类似地也可以采用Crout 分解):()1.4.2321111153342231124321432⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡---⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡=⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡------u u u u l l l 即L 为单位上三角阵,两斜行,主对角线元素为1,其下方的斜行元素特定;U 为上三角阵,也是两斜行,主对角线元素特定,其上方斜行的元素与A 对应的斜行元素相同(直接验算可知道)。

利用矩阵乘法规则,按顺序依次考虑A 的11a ->21a ->22a ->32a ->33a ->43a ->44a ,并对比(2.4.1)式两端可得2=1u → 1u =2 -1=12u l → 2l =-1/1u =-1/23=-2l +2u →2u =3+2l =5/2-2=23u l → 3l =-2/2u =-4/5 4=-23l +3u → 3u =4+23l =12/5-3=34u l →4l =-3/3u =-5/45=-34l +4u 4u =5+3⨯(-5/4)=5/4即得分解⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡---⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡---=⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡------4/535/1222/51214/515/412/115334223112于是用前推过程求解下三角方程组⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡-=⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡---121614/515/412/114321y y y y 得⎪⎪⎩⎪⎪⎨⎧=+==+-==+==2/54/515/65/4242/1163423121y y y y y y y再用回代过程求解上三角方程组⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡=⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡---2/55/6464/535/1222/5124321x x x x 得⎪⎪⎩⎪⎪⎨⎧=+==+==+===52/)6(4)2/5/()24(3)5/12/()35/6(2)4/5/()2/5(2132434x x x x x x x 即的方程组的解()T x 2,3,4,5=.从实例看到,三对角方程组的追赶法是三角分解发的一种特殊应用,因此,一般地,如果对三角矩阵n n R A ⨯∈非奇异,其顺序主子式)1,...,2,1(0-=≠∆n i i ,则解三对角方程组Ax=d:()2.4.2...............12112111122211⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡=⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡-----n n n n n nn n n d d d d x x x x b a c b a c b a c b 的追赶法可描述如下:令A=LU ,则⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡=⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡-----n n n nn nn n n u c u c u c u l l l b a c b a c b a c b 1122113211122211......1 (1)11.........利用矩阵乘法规则,可求L 和U 的计算公式:⎩⎨⎧=-===--)),...3,2((/1111n i c l b u u a l b u i i i i i i i (2.4.3)于是,求解LY=D 得⎩⎨⎧=-==-)),...3,2(111n i y l d y d y i i i i (2.4.4)再求解UX=Y ,得三对方程组的解⎩⎨⎧=-==+)),...3,2(/)(/1n i u x c y x u y x i i i i i n n n (2.4.5)上述3个公式便组成解三角方程组的追赶法,国外称Thomas 算法。

(完整版)2.6追赶法

(完整版)2.6追赶法

第16讲 追赶法、误差分析在实际应用问题中,经常会遇到解三对角线方程组。

例如:用三次样条函数的插值问题中得到的三转弯及三弯矩方程组,当时说可用追赶法来求解。

还有用差分法解二阶线性常微分方程边值问题,若用三点插值格式也得到解三对角线方程组,本节介绍该类方程组中的特例及该种方程组的解法:追赶法。

优点:1.计算量小。

2.方法简单,存贮量小。

3.数值稳定的(对舍入误差来说)。

1 追赶法三对角线方程组的一般表示方法:可见,对A 的分解只需求i i u l ,且按n n n l u l u l u l −→−−→−−→−−→−−→−−→−−→−--112211.....的递推过程进行,形象地称为“追”的过程⎩⎨⎧=-==-),....2(/)(/1111n i l y a f y l f y i i i i⎩⎨⎧-=-==+)1,2,.....1(1n i x u y x y x i i i inn 形象地称回代求解过程为“赶”的过程追赶法的计算量为5n-4次乘除法,可用4个 一 维数组存放{}{}{}{}i i i i f c b a ,,,。

共占用4n-2个单元,在计算过程中{}{}{}i i i y u l ,,依次覆盖掉{}{}{}i i i f c b ,,最后,{}i x 覆盖掉{}i y ,所以,追赶法具有计算量小,占用内存单元少的特点。

2、误差分析⎪⎪⎪⎪⎪⎪⎪⎪⎭⎫⎝⎛=-n n l u u u U 121....111⎪⎪⎪⎪⎪⎪⎪⎪⎭⎫⎝⎛=n nl a l a l a l L ....33221)1,...,3,2,1(-=n i ⎪⎩⎪⎨⎧+===+++11111i i i i ii i lu a b u l c l b ⎪⎩⎪⎨⎧-===+++ii i i ii i u a b l l c u b l 11111/)1,...,3,2,1(-=n i病态方程组与条件数一个线性方程组Ax=b 是由它的系数矩阵A 和它的右端项b 所确定,在实际问题中,由于各种原因,A 或b 往往有误差,从而使得解也产生误差。

(整理)线性方程组的直接法

(整理)线性方程组的直接法

第二章线性方程组的直接法在近代数学数值计算和工程应用中,求解线性方程组是重要的课题。

例如,样条插值中形成的关系式,曲线拟合形成的法方程等,都落实到解一个元线性方程组,尤其是大型方程组的求解,即求线性方程组(2.1)的未知量的数值。

(2.1)其中ai j,bi为常数。

上式可写成矩阵形式Ax = b,即(2.2)其中,为系数矩阵,为解向量,为常数向量。

当detA=D0时,由线性代数中的克莱姆法则,方程组的解存在且惟一,且有为系数矩阵的第列元素以代替的矩阵的行列式的值。

克莱姆法则在建立线性方程组解的理论基础中功不可没,但是在实际计算中,我们难以承受它的计算量。

例如,解一个100阶的线性方程组,乘除法次数约为(101·100!·99),即使以每秒的运算速度,也需要近年的时间。

在石油勘探、天气预报等问题中常常出现成百上千阶的方程组,也就产生了各种形式方程组数值解法的需求。

研究大型方程组的解是目前计算数学中的一个重要方向和课题。

解方程组的方法可归纳为直接解法和迭代解法。

从理论上来说,直接法经过有限次四则运算,假定每一步运算过程中没有舍入误差,那么,最后得到方程组的解就是精确解。

但是,这只是理想化的假定,在计算过程中,完全杜绝舍入误差是不可能的,只能控制和约束由有限位算术运算带来的舍入误差的增长和危害,这样直接法得到的解也不一定是绝对精确的。

迭代法是将方程组的解看作某种极限过程的向量极限的值,像第2章中非线性方程求解一样,计算极限过程是用迭代过程完成的,只不过将迭代式中单变量换成向量而已。

在用迭代算法时,我们不可能将极限过程算到底,只能将迭代进行有限多次,得到满足一定精度要求的方程组的近似解。

在数值计算历史上,直接解法和迭代解法交替生辉。

一种解法的兴旺与计算机的硬件环境和问题规模是密切相关的。

一般说来,对同等规模的线性方程组,直接法对计算机的要求高于迭代法。

对于中等规模的线性方程组,由于直接法的准确性和可靠性高,一般都用直接法求解。

线性方程组的解法及其应用

线性方程组的解法及其应用

线性方程组的解法及其应用摘要:线性方程组是线性代数的核心内容之一,其解法研究是代数学中经典且重要的研究课题.本文综述了几种不同类型的线性方程组的解法,如消元法、克拉默法则、广义逆矩阵法、直接三角形法、平方根法、追赶法,并以具体例子介绍不同解法的应用技巧. 在这些解法中,广义逆矩阵方法,具有表达式清晰,使用范围广的特点.另外,这些方法利于快速有效地解决线性方程组的求解问题,为解线性方程组提供一个简易平台,促进了理论与实际的结合.关键词:线性方程组解法广义逆矩阵应用实例1. 引言线性方程组理论是高等数学中十分重要的内容,而线性方程组的解法是利用线性方程组理论解决问题的关键.本文主要介绍线性方程组的广义逆矩阵法、追赶法、平方根法等求解方法,为求解线性方程组提供一个平台.文章也给出线性方程组在其他领域中的应用实例,揭示了各学科之间的内通性.首先,我们讨论一般线性方程组.这里所指的一般线性方程组形式为11112211211222221122,,.n n n n s s sn n s a x a x a x b a x a x a x b a x a x a x b +++=⎧⎪+++=⎪⎨⎪⎪+++=⎩ ()i()i 式中(1,2,,)i x i n =代表未知量,(1,2,,;1,2,,)ij a i s j n ==称为方程组的系数,(1,2,,)j b j n =称为常数项.线性方程组)(i 称为齐次线性方程组,如果常数项全为零,即120s b b b ====.令111212122212n n s s sn a a a a a a A a a a ⎡⎤⎢⎥⎢⎥=⎢⎥⎢⎥⎣⎦,12n x x X x ⎡⎤⎢⎥⎢⎥=⎢⎥⎢⎥⎣⎦, 12s b b B b ⎡⎤⎢⎥⎢⎥=⎢⎥⎢⎥⎣⎦,则()i 可用矩阵乘法表示为AX B =,,,.m n n m A C X C B C ⨯∈∈∈2. 线性方程组的解法2.1 消元法在初等代数里,我们已经学过用代入消元法和加减消元法解简单的二元、三元线性方程组.实际上,这个方法比用行列式解方程组更具有普遍性.但对于那些高元的线性方程组来说,消元法是比较繁琐的,不易使用.例 1 解线性方程组123123123123324,32511,23,237.x x x x x x x x x x x x +-=⎧⎪+-=⎪⎨++=⎪⎪-++=-⎩ 解 分别将第一个方程的(-3)倍,(-2)倍和2倍加到第二、三、四个方程上,整理得123232323324,71,555,7 1.x x x x x x x x x +-=⎧⎪-+=-⎪⎨-+=-⎪⎪-=⎩将此方程组第二个方程加到第四个方程上,使该方程两边全为零,并将第三个方程的两边乘以15-,得1232323324,71,1.x x x x x x x +-=⎧⎪-+=-⎨⎪-=⎩再将第三个方程的7倍加到第二个方程上,消去第二个方程中的未知量2x ,整理得123233324,1,6 6.x x x x x x +-=⎧⎪-=⎨⎪-=⎩最后解得123(,,)(2,0,1)T T x x x =--.正如消元法是我们接触比较早的,被我们所熟悉的一种方法,在此只给出三元线性方程组的解法,三元以上的方程组的具体理论、性质和解题过程详见参考文献[1]. 2.2 应用克莱姆法则对于未知个数与方程个数相等的情形,我们有定理1[1] 如果含有n 个方程的n 元线性方程组11112211211222221122,,.n n n n n n nn n n a x a x a x b a x a x a x b a x a x a x b +++=⎧⎪+++=⎪⎨⎪⎪+++=⎩ ()ii的系数矩阵111212122212n n n n nn a a a a a a A a a a ⎡⎤⎢⎥⎢⎥=⎢⎥⎢⎥⎣⎦的行列式111212122212det 0n n n n nna a a a a a A a a a =≠,那么线性方程组()ii 有唯一解:det (1,2,,),det j j B x j n A==其中det j B 是把矩阵中第j 列换成线性方程组的常数项12,,,n b b b 所成的矩阵的行列式,即111,111,11222,122,121,1,1det,1,2,,.j j n j j n j n n j n n j nna ab a a a a b a a B j n a a b a a -+-+-+==此外,还可以叙述为,如果含有n 个未知数、n 个方程的线性方程组Ax b =的系数矩阵的行列式det 0A ≠,则线性方程组Ax b =一定有解,且解是唯一的. 例2 解线性方程组12342341242342344,3,31,73 3.x x x x x x x x x x x x x -+-=⎧⎪-+=-⎪⎨++=⎪⎪-++=-⎩ 解 由已知可得系数行列式12341234123401110111111det 16013015352073173148A ---------====≠----,因此线性方程组有唯一解.又因124234143431110311det 128,det 48,1301110137310331B B -------==-==-341244123401310113det 96,det 0.1311130107310733B B ------====--故线性方程组的解为1234(,,,)(8,3,6,0)T T x x x x =-.克莱姆法则主要给出了解与系数的明显关系,但只能应用于系数矩阵的行列式不为零的线性方程组,并且它进行计算是不方便的. 2.5 直接三角分解法[5]设有线性方程组11112211211222221122,,,n n n n n n nn n n a x a x a x b a x a x a x b a x a x a x b +++=⎧⎪+++=⎪⎨⎪⎪+++=⎩或写成矩阵形式Ax b =,其中111212122212n n n n nn a a a a a a A a a a ⎡⎤⎢⎥⎢⎥=⎢⎥⎢⎥⎣⎦,12n x x x x ⎡⎤⎢⎥⎢⎥=⎢⎥⎢⎥⎣⎦,12n b b b b ⎡⎤⎢⎥⎢⎥=⎢⎥⎢⎥⎣⎦.若A 为非奇异矩阵,且有分解式A LU =,其中U 为上三角矩阵,L 为单位下三角矩阵,即11121212221,1111n n n n n nn u u u l u u A LU l l u -⎡⎤⎡⎤⎢⎥⎢⎥⎢⎥⎢⎥==⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦⎣⎦, 则线性方程组Ax b =的求解等价于 解以下两个三角方程组:(1)Ly b =,求y ; (2)Ux y =,求x .直接三角形分解法求解线性方程组,基本步骤如下: 第一步: 11,(1,2,,),i i u a i n == 1111,(2,3,,)i i l a u i n ==,计算U 的第r 行,L 的第r 列元素,2,3,,r n =.第二步: 11,(,1,,)r ri ri rk ki k u a l u i r r n -==-=+∑.第三步: 11,(1,,;)r ir ir ik kr rr k l a l u u i r n r n -==(-)=+≠∑.求解Ly b =,Ux y =的计算公式如下:第四步: ()1111,,2,3,.i i i ik k k y b y b l y i n -==⎧⎪⎨=-=⎪⎩∑第五步: 1,(),(1,2,,1).n n nn n i i ik k ii k i x y u x y u x u i n n =+=⎧⎪⎨=-=--⎪⎩∑例5 求解线性方程组1231212321,42,227.x x x x x x x x ++=⎧⎪+=-⎨⎪-++=⎩解 由直接三角分解法第二、三步可得211100211410210012221131004A LU ⎡⎤⎡⎤⎡⎤⎢⎥⎢⎥⎢⎥==--=⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥---⎣⎦⎣⎦⎣⎦. 于是线性方程组变为LUx b =,求解线性方程组(1,2,7)T Ly =-,得(1,4,4)T y =--;求解线性方程组(1,4,4)T Ux =--,得(1,2,1)T x =-.2.6 平方根法[7]在许多应用中,欲求解的线性方程组的系数矩阵是对称正定的.所谓平方根法,就是利用对称正定矩阵的三角分解而得到的求解具有对称正定矩阵的线性方程组的一中有效方法,目前在计算机上广泛应用平方根法解此类方程组.定理6[12] 若A 的各阶顺序主子式非零,则A 可以分解为A LDU =,其中L 是单位下三角矩阵,U 是单位上三角矩阵,D 是对角矩阵,且这种分解是唯一的.定理7[12] 设A 为对称正定矩阵,则存在三角分解T A LL =,其中L 是非奇异下三角形矩阵,且当限定L 的对角线元素为正时,这种分解是唯一的.应用对称正定矩阵的平方根法,可以解具有对称正定系数矩阵的线性方程组Ax b =,具体算法如下:1) 对j =1,2,,n ,计算11221()j jj jj jkk l a l -==-∑,11j ij ij ik jk k l a l l -==-∑(1,,)i j n =+.2) 求解线性方程组Ax b =等价于解两个三角方程组,.TLy b L x y =⎧⎨=⎩ 计算11()i i i ik k ii k y b l y l -==-∑,(i =1,2,,n ), 1()ni i ki kii k i x b lx l =+=-∑,(i n =,1n -,,2,1),即可.例6 求解线性方程组12341161 4.25 2.750.5.1 2.75 3.5 1.25x x x -⎡⎤⎡⎤⎡⎤⎢⎥⎢⎥⎢⎥-=-⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦⎣⎦⎣⎦ 解 设1111213121222232313233334111 4.25 2.751 2.75 3.5l l l l l l l l l l l l -⎡⎤⎡⎤⎡⎤⎢⎥⎢⎥⎢⎥-=⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦⎣⎦⎣⎦, 由矩阵乘法得1121223132332,0.5,2,0.5, 1.5, 1.l l l l l l ==-====解下三角方程组123260.520.50.5 1.51 1.25y y y ⎡⎤⎡⎤⎡⎤⎢⎥⎢⎥⎢⎥-=-⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦⎣⎦⎣⎦, 得1233,0.5,1,y y y ===-再由123230.520.50.5 1.511Tx x x ⎡⎤⎡⎤⎡⎤⎢⎥⎢⎥⎢⎥-=⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥-⎣⎦⎣⎦⎣⎦, 得线性方程组的解为123(,,)(2,1,1)T T x x x =-.可以用消元法解此方程组,但发现此方程组的系数矩阵为正定矩阵,运用平方根法解这个方程组比较容易,而且理论分析指出,解对称正定方程组的平方根法是一个稳定的算法,其在工程计算中使用比较广泛. 2.7 追赶法[5]在许多实际问题中,都会要求解系数矩阵为对角占优的三对角方程组11112222211111iiii i n n n n n nn n n x k b c x k a b c a b c x k a b c x k a b x k -----⎡⎤⎡⎤⎡⎤⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥=⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦⎣⎦⎣⎦, 简记作 Ax k =, 其中A 满足下列对角占优条件:(1) 110b c >>;(2) i i i b a c ≥+, i a ,i c 0≠(i =2,3, ,1n -);(3) 0n n b c >>.由系数矩阵A 的特点,可以将A 分解为两个三角矩阵的乘积,即A LU =,其中L 为下三角矩阵,U 为单位上三角矩阵.求解线性方程组Ax k =等价于解两个三角方程组Ly k =与Ux y =,先后求y 与x ,从而得到以下解三角方程组的追赶法公式:第一步:计算的递推公式111c b β=,1()i i i i i c b a ββ-=-,(2i =,3,,1)n -;第二步:解Ly k =:111y k b =,11()()i i i i i i i y k a y b a β--=--,(2,3,,)i n =;第三步:解Ux y =:n n x y =,1i i i i x y x β+=-,(1,2,,2,1)i n n =--.例7 求解三对角线性方程组123421001131020111200210x x x x ⎡⎤⎡⎤⎡⎤⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥=⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦⎣⎦⎣⎦.解 设有三角分解111122222233333344441111b c p q a b c a p q a b c a p q a b a p ⎡⎤⎡⎤⎡⎤⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥=⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦⎣⎦⎣⎦, 由矩阵乘法易得111,,1,2,3.,2,3,4.i i i ii i i p b q c p i p b a q i -=⎧⎪==⎨⎪=-=⎩ 将已知系数矩阵的元素代人上式有11223342,12,52,25,35,53,73.p q p q p q p ==⎧⎪==⎪⎨==⎪⎪=⎩ 解线性方程组112233441121220p y p y p y p y ⎡⎤⎡⎤⎡⎤⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥=⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦⎣⎦⎣⎦, 得123412,35,73, 2.y y y y ====再解线性方程组111222333441111x y q x y q x y q x y ⎡⎤⎡⎤⎡⎤⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥=⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦⎣⎦⎣⎦,得原线性方程组的为1234(,,,)(0,1,1,2)T T x x x x =-.追赶法是以LU 分解为基础的求解方法,因此它的不足之处是当某个0=k u 时,就不能进行.但是当方程组的系数矩阵A 中有很多零元素时,利用三对角方程组系数矩阵的稀疏性,使零元素不参加运算,可以类似于追赶法来简化计算过程,从而极大地节省了计算量和存储量.这也是追赶法的最大特点.3. 应用举例3.1 线性方程组在解析几何中的应用例8 已知平面上三条不同直线的方程分别为1L :230ax by c ++=,2L :230bx cy a ++=,3L :230cx ay b ++=,试证:这三条直线交于一点的充分必要条件为0a b c ++=.证 必要性 设三直线1L ,2L ,3L 交于一点,则线性方程组232323ax by cbx cy a cx ay b +=-⎧⎪+=-⎨⎪+=-⎩ ()iii有惟一解,故系数矩阵222a b A b c c a ⎡⎤⎢⎥=⎢⎥⎢⎥⎣⎦与增广矩阵232323a b c A b c a c a b --⎡⎤⎢⎥=-⎢⎥⎢⎥-⎣⎦的秩均为2,于是0A -=,即22223236()()23a bcA bc a a b c a b c ab ac bc ca b--=-=++++----=0,所以0a b c ++=.充分性 由0a b c ++=,则从必要性的证明可知,0A -=,故()3r A -<.由于22222132()2[()]2[()]0224a b ac b a a b b a b b b c =-=-++=-++≠, 故()()2r A r A -==.因此线性方程组()iii 有惟一解,即三直线1L ,2L ,3L 交于一点. 3.2 线性方程组在产品生产量中的应用例9 设有一个经济系统包括3个部门,在某一个生产周期内各部门间的消耗及最终产品如表所示:求各部门的总产品.解 设i x 表示第i 部门的总产品.由已知可以得到线性方程组()I A x y -=,其中0.250.10.1()0.20.20.10.10.10.2ij A a ⎡⎤⎢⎥==⎢⎥⎢⎥⎣⎦,0.750.10.10.20.80.10.10.10.8I A --⎡⎤⎢⎥-=--⎢⎥⎢⎥--⎣⎦,(245,90,175)T y =. 利用矩阵的初等变换可以求得1126181810()34118198912017116I A -⎡⎤⎢⎥-=⎢⎥⎢⎥⎣⎦, 所以线性方程组()I A x y -=的解为消耗系数 消耗部门 生产部门123最终产品1 0.25 0.1 0.1 2452 0.2 0.2 0.1 90 30.10.10.21751126181824540010()3411819902508912017116175300x I A y -⎡⎤⎡⎤⎡⎤⎢⎥⎢⎥⎢⎥=-==⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦⎣⎦⎣⎦. 4. 结束语本文针对不同的线性方程组给出了一些计算方法,及线性方程组的应用实例.根据线性方程组自身所具有的特点,可以选择相应合适的方法,而对于那些特殊类型的线性方程组的解法,有待进一步的讨论与研究.参考文献:[1] 北京大学数学系几何与代数教研室前代数小组编. 高等代数[M].3版.北京:高等教育出版社,2003.105-112.[2] 白梅花. 线性方程组若干应用实例举例[J].科技资讯,2011,(27):200-201.[3] 康道坤,陈劲. 广义逆下线性方程组的解结构及其推广[J].大理学院学报,2011,10(4):7-9. [4] 卢刚.线性代数[M]. 北京:高等教育出版社,2002.64-72.[5] 李庆扬,王能超,易大义. 数值分析[M].4版.武汉:华中科技大学出版社,2006.177-185. [6] 苏育才,姜翠波,张跃辉. 矩阵理论[M].北京:科学出版社,2006.200-206. [7] 首都师范大学数学系组编. 数值分析[M].北京:科学出版社,2000.28-32.[8] 徐仲,张凯院,陆全,等. 矩阵论简明教程[M].2版.北京:科学出版社,2005.141-147. [9] 谢寿才,陈渊. 大学数学[M].北京:科学出版社,2010.37-40.[10] 徐仲,张凯院,陆全. 矩阵论[M].西安:西北工业大学出版社,2002.228-245.[11] 尹钊,钟卫民,赵丽君. 线性方程组的广义逆矩阵解法[J].哈尔滨师范大学自然科学学 报,1999,15(5):21-22. [12] 张明淳. 工程矩阵理论[M].1版.南京:东南大学出版社,1995.172-173.[13] 赵树嫄. 线性代数(经济应用数学基础)[M].4版.北京:中国人民大学出版社,2008.150-157.。

Ch3.2.3 追赶法

Ch3.2.3  追赶法

1 b1 ( 1 ) , c b 1 1 1
解(7.1)的追赶法计算公式 1 r (1)分解计算公式( A LU): 2 2 1 c1 b1 a ai i i )) , ( i 2, , n 1) i ri i 11 i ci (bii (2)求解Ly f 逆推公式
三次样条插值问题中得到的三转弯方程组(8.9)、 (8.10)、(8.11)及三弯矩方程组(8.19)或(8.20) 三对角线 方程组
用差分法解二阶线性常微分方程边值问题, 解法:追赶法 若用三点插值格式也得到三对角线方程组 本节介绍该类方程组中的特例及该种方程组的解法:追赶法。 优点: 1.计算量小(仅5n-4次乘除法运算)。 2.方法简单,存贮量小。 3.数值稳定(对舍入误差来说)。
LU
1 b1 c1 1 a b c r 2 2 2 2 2 a b c r aii bii cii rii i i an bn rn n n a l u 由矩阵乘法 ij ik kj , 得:
b3 a 4
c3 b4
b1 c1 x1 f 1 x3 f 3 x f 与 a b x f c x 。 2 3 2 2 2 2 4 4
§7 解三对交线方程组的追赶法
k 1
1 1 2 1 i2 11 ii 11 i 1 i 1 1 n 1 1
1 b1 , (1) b1 1 ,c1 1 1 , ( 1 c1 1 ) 1 c1 b1

数值分析列主元高斯消去顺序高斯平方根法追赶法

数值分析列主元高斯消去顺序高斯平方根法追赶法

课题名称:课题一解线性方程组的直接方法解决的问题:给定三个不同类型的线性方程组,用适当的直接法求解。

采用的数值方法:对第一个普通的线性方程组,采用了高斯顺序消去法和高斯列主元消去法。

对第二个正定线性方程组,采用了平方根法。

对第三个三对角线性方程组,采用了追赶法。

算法程序:(1)普通的线性方程组①顺序消去法#include<stdio.h>#include<math.h>int main(void){float A[10][10]= {{4,2,-3,-1,2,1,0,0,0,0},{8,6,-5,-3,6,5,0,1,0,0},{4,2,-2,-1,3,2,-1,0,3,1},{0,-2,1,5,-1,3,-1,1,9,4},{-4,2,6,-1,6,7,-3,3,2,3},{8,6,-8,5,7,17,2,6,-3,5},{0,2,-1,3,-4,2,5,3,0,1},{16,10,-11,-9,17,34,2,-1,2,2},{4,6,2,-7,13,9,2,0,12,4},{0,0,-1,8,-3,-24,-8,6,3,-1}};float b[10]= {5,12,3,2,3,46,13,38,19,-21}; float x[10]= {0};float Aik,S,temp;int i,j,k;int size=10;for(k=0; k<size-1; k++){if(!A[k][k])return -1;for(i=k+1; i<size; i++){Aik=A[i][k]/A[k][k];for(j=k; j<size; j++){A[i][j]=A[i][j]-Aik*A[k][j]; }b[i]=b[i]-Aik*b[k];}}printf("A[]\n");for(i=0; i<size; i++){for(j=0; j<size; j++)printf("%f ",A[i][j]);printf("\n");}printf("b[]\n");for(i=0; i<size; i++)printf("%f ",b[i]);printf("\n\n");x[size-1]=b[size-1]/A[size-1][size-1]; for(k=size-2; k>=0; k--){S=b[k];for(j=k+1; j<size; j++){S=S-A[k][j]*x[j];}x[k]=S/A[k][k];}printf("x[]=\n");for(i=0; i<size; i++)printf("%f ",x[i]);return 0;}②列主元消去法#include<stdio.h>#include<math.h>int main(void){float A[10][10]= {{4,2,-3,-1,2,1,0,0,0,0}, {8,6,-5,-3,6,5,0,1,0,0},{4,2,-2,-1,3,2,-1,0,3,1},{0,-2,1,5,-1,3,-1,1,9,4},{-4,2,6,-1,6,7,-3,3,2,3},{8,6,-8,5,7,17,2,6,-3,5},{0,2,-1,3,-4,2,5,3,0,1},{16,10,-11,-9,17,34,2,-1,2,2},{4,6,2,-7,13,9,2,0,12,4},{0,0,-1,8,-3,-24,-8,6,3,-1}};float b[10]= {5,12,3,2,3,46,13,38,19,-21}; float x[10]= {0};float Aik,S,temp;int i,j,k;float max;int col;int size=10;for(k=0; k<size-1; k++){max=fabs(A[k][k]);col=k;for(i=k; i<size; i++){if(max<fabs(A[i][k])) {max=fabs(A[i][k]); col=i;}}for(j=k; j<size; j++){temp=A[col][j];A[col][j]=A[k][j];A[k][j]=temp;}temp=b[col];b[col]=b[k];b[k]=temp;if(!A[k][k])return -1;for(i=k+1; i<size; i++){Aik=A[i][k]/A[k][k]; for(j=k; j<size; j++){A[i][j]=A[i][j]-Aik*A[k][j]; }b[i]=b[i]-Aik*b[k];}}printf("A[]\n");for(i=0; i<size; i++){for(j=0; j<size; j++)printf("%f ",A[i][j]);printf("\n");}printf("b[]\n");for(i=0; i<size; i++)printf("%f ",b[i]);printf("\n\n");x[size-1]=b[size-1]/A[size-1][size-1]; for(k=size-2; k>=0; k--){S=b[k];for(j=k+1; j<size; j++){S=S-A[k][j]*x[j]; }x[k]=S/A[k][k];}printf("x[]=\n");for(i=0; i<size; i++)printf("%f ",x[i]);return 0;}(2)对称正定线性方程组平方根法:#include <stdio.h>#include <math.h>#define n 8int main(void){float A[8][8]={{4,2,-4,0,2,4,0,0},{2,2,-1,-2,1,3,2,0},{-4,-1,14,1,-8,-3,5,6},{0,-2,1,6,-1,-4,-3,3},{2,1,-8,-1,22,4,-10,-3},{4,3,-3,-4,4,11,1,-4},{0,2,5,-3,-10,1,14,2},{0,0,6,3,-3,-4,2,19}};float g[8][8]= {0};float b[8]= {0,-6,6,23,11,-22,-15,45}; float x[8]= {0};float y[8]= {0};int k,m,i,sq;for(k=0; k<n; k++){float p=0,q=0,s=0;for(m=0; m<=k-1; m++){p=p+A[k][m]*A[k][m];}g[k][k]=sqrt(A[k][k]-p);A[k][k]=g[k][k];for(i=k+1; i<n; i++){q=0;for(m=0; m<=k-1; m++){q=q+A[i][m]*A[k][m];}g[i][k]=(A[i][k]-q)/A[k][k]; A[i][k]=g[i][k];}s=0;for(m=0; m<=k-1; m++){s=s+A[k][m]*y[m];}y[k]=(b[k]-s)/A[k][k];}x[n-1]=y[n-1]/A[n-1][n-1];for(k=n-2; k>=0; k--){float sum=0;for(m=k+1; m<n; m++){sum=sum+A[m][k]*x[m];}x[k]=(y[k]-sum)/A[k][k];}for(sq=0; sq<n; sq++){printf("%f ",x[sq]);}return 0;}(3)三对角线性方程组追赶法#include <stdio.h>#include <math.h>#define n 10int main(void){float a[10]={4,4,4,4,4,4,4,4,4,4};float c[9]={-1,-1,-1,-1,-1,-1,-1,-1,-1};float d[9]={-1,-1,-1,-1,-1,-1,-1,-1,-1}; float b[10]={7,5,-13,2,6,-12,14,-4,5,-5}; float x[10]={0};float y[10]={0};float arf[10]={0};float bt[9]={0};arf[0]=a[0];int i;for(i=0;i<n-1;i++){bt[i]=c[i]/arf[i];arf[i+1]=a[i+1]-d[i+1]*bt[i];//printf("%f %f \n",bt[i],arf[i+1]); }y[0]=b[0]/arf[0];//printf("%f\n",y[0]);for(i=1;i<n;i++){y[i]=(b[i]-d[i]*y[i-1])/arf[i];}//printf("%f\n",y[1]);x[n-1]=y[n-1];for(i=n-2;i>=0;i--){x[i]=y[i]-bt[i]*x[i+1]; }for(i=0;i<n;i++)printf("%lf ",x[i]);return 0;}数值结果:(1) 普通的线性方程组①顺序消去法②列主元消去法(2) 对称正定线性方程组平方根法:(3)三对角线性方程组追赶法:对实验计算结果的讨论和分析:(1) 普通的线性方程组①顺序消去法x1~x10的绝对误差:0.000001,-0.000001,0.000001,0,0.000001,0,0.000002,0,0,0x1~x10的相对误差:0.000001,0.000001,-1,0,0.0000005,0,0.00000067,0,0,0误差很小,基本可以忽略。

追赶法解三对角方程组

追赶法解三对角方程组

《数值分析》课程设计追赶法解三对角方程组院(系)名称信息工程学院专业班级10普本信计学号100111014学生姓名刘银朋指导教师张荣艳2013 年05 月31日数值分析课程设计评阅书题目追赶法解三对角方程组学生姓名刘银朋学号100111014 指导教师评语及成绩指导教师签名:年月日答辩评语及成绩答辩教师签名:年月日教研室意见总成绩:教研室主任签名:年月日课程设计任务书2012—2013学年第二学期专业班级:10普本信息与计算科学学号:100111014 姓名:刘银朋课程设计名称:数值分析Ⅰ、Ⅱ设计题目:追赶法解三对角方程组完成期限:自2013 年05月21 日至2013年05 月31日共10天设计依据、要求及主要内容:一、设计目的理解追赶法,掌握追赶法的算法设计以及关于追赶法的分析和综合应用,能够较熟练的应用Matlab软件编写求解追赶法的程序和应用Matlab软件数据库软件.二、设计内容(1)认真挑选有代表性的三对角方程组.(2)认真梳理解三对角方程组的解题思路.(3)比较追赶法和高斯消去法的计算精度.三、设计要求1.先用Matlab数据库中的相应的函数对选定的方程,求出具有一定精度的解. 2.然后使用所用的方法编写Matlab程序求解.3.对于使用多个方程解同意问题的,在界面上要设计成菜单的形式.计划答辩时间:2013年06 月5 日工作任务鱼工作量要求:查阅文献资料不少于3篇,课程设计报告1篇不少于3000字.指导教师(签字):教研室主任(签字):批准日期:2013 年05 月20 日追赶法解三对角方程组摘要本文主要通过运用追赶法来求解三对角方程组的问题.追赶法是用来求解三对角方程组的专用方法,对于三对角方程组,追赶法比Gauss消去法的计算量要小的多,本文主要介绍了追赶法的原理,并用Matlab编写求解程序,以实现对三对角方程组的求解,进一步解决实际中的问题.并且根据所得出的结果分析追赶法算法和高斯消去的法的计算精度.关键词:追赶法,三对角方程组,追赶法的Matlab程序目录1.前言 (1)2.解题思想和方法 (1)2.1 追赶法解题思想 (2)2.2追赶法解题原理 (4)3.对追赶法的MATLAB求解 (5)3.1实验程序 (5)3.2 应用举例 (6)4.与高斯消去法的精度比较 (7)课程设计总结 (9)参考文献 (10)1.前言当今很多科学与工程计算问题大都可以化为线性方程组的形式,所以有效的求解线性方程组在科学和工程计算中是非常重要的.虽然线性代数方程的求解方法和数值计算软件包均很成熟,但随着并行计算机的发展,问题的求解速度和解题规模都大大提高,因而使数值计算方法和响应的数学软件包都产生了变化,相应的线性方程组的有效并行求解也引起了人们的普遍关注.追赶法是用来求解三对角方程组的专用方法,生活中很多实际问题,都归结为求解线性方程组.例如解常微分方程边值问题,解热传导方程以及船体数学放样中建立三次样条函数等,都会要求就解系数矩阵成三对角线性的方程组.而解三对角方程组的最简单方法是追赶法,公式简单,计算量小,所占用的存储单元少,所以在小机器上也能求解.追赶大事用来求解三对角方程组的专用方法,对于三对角方程组,追赶法比Gauss消去法的计算量要小得多.应用追赶法求解三对角线性方程,追赶法仍然保持LU分解特性,它是一种特殊的LU分解.充分利用了系数矩阵的特点,而且使之分解更简单,得到对三对角线性方程组的快速解法. 本文讨论使用追赶法解线性方程组.介绍追赶法的理论,求解线性方程组的追赶法的实现以及追赶法的应用.2.解题思想和方法三对角矩阵是一种具有特殊意义的带状矩阵.用差分法求解二阶常微分方程边值问题时,最后常规解为求解具有三对角系数矩阵的线性方程组.对三对角矩阵实行Doolittle (或Crout)分解,便得到求解三对角方程组的最有效方法---追赶法.设矩阵A非奇异,A有Crout分解A=LU,其中L为下三角矩阵,U为单位上三角矩阵,可先依次求出L,U中的元素后,令Ux=y,先求解下三角方程组Ly=f得出y,再求解上三角方程组Ux=y,事实上求解三对角方程组的2追赶法将矩阵三角分解的计算与求解两个三角方程组的计算放在一起,使算法更为紧凑.本文主要介绍追赶法的原理,以及编写Matlab 程序实现在计算机上的应用并分析他们的计算精度,比较解三对角方程组的最优解的问题. 2.1 追赶法解题思想在实际问题中,经常遇到以下形式的方程组⎪⎪⎪⎪⎩⎪⎪⎪⎪⎨⎧=+=++=++=++=+-------+-n n n n n n n n n n n n k k k k k k k d x b x a d x c x b x a d x c x b x a d x c x b x a d x c x b111112111232221212111 (2.11)这种方程组的系数矩阵A 为三对角矩阵,即⎪⎪⎪⎪⎪⎪⎪⎪⎭⎫⎝⎛=---n nn n n k k k b a c b a c b a c b a c b A 11122211 以下针对这种方程组的特点提供一种简便有效的算法—追赶法.追赶法实际上是高斯消去法的一种简化形式,它同样分消元与回代两个过程.先将(2.11)第一个方程中x 1的系数化为1112111b dx b c x =+111111b d y b c r ==(2.12)有1211y x r x =+注意到剩下的方程中,实际上只有第二个方程中含有变量x 1,因此消元手续可以简化.利用(2.12)可将第二个方程化为2312y x r x =+,这样一步一步地顺序加工(2.11)的每个方程,设第k – 1个方程已经变成111---=+k k k k y x r x(2.13)再利用(2.13)从第k 个方程中消去x k -1,得:k k k k k k k k k a y d x c x a r b 111)(-+--=+-同除()k k k a r b 1--,得n k a r b a y d x a r b c x kk k kk kk k k k k k ,,3,21111 =--=-+--+-记kk k kk k k kk k kk a r b a y d y a r b c r 111-----=-=则有 k k k k y x r x =++1 这样做n – 1步以后,便得到:111---=+n n n n y x r x将上式与(2.11)中第11个方程联立,即可解出x n = y n ,这里nn n nn n n a r b a y d y 11----=于是,通过消元过程,所给方程组(2.11)可归结为以下更为简单的形式:⎪⎪⎪⎩⎪⎪⎪⎨⎧==+=++nn k k k k y x y x r x y x r x 11211 (2.14)这种方程组称作二对角型方程组,其系数矩阵中的非零元素集中分步在主对角线和一条次主对角线上⎪⎪⎪⎪⎪⎪⎭⎫⎝⎛-11111121n k r r r r 对加工得到的方程组(2.14)自下而上逐步回代,即可依次求出x n ,x n -1,…,x 1,计算公式为:⎩⎨⎧--=-==+1,,2,11n n k x r y x y x k k k k nn (2.15)上述算法就是追赶法,它的消元过程与回代过程分别称作“追”过程与“赶”过程.综合追与赶的过程,得如下计算公式:⎪⎪⎪⎩⎪⎪⎪⎨⎧=--=-===---nk a r b a y d y a r b c r b d y b c r k k k k k k k kk k k k ,,3,2111111111(2.16)⎩⎨⎧--=-==+1,,2,11n n k x r y x y x k k k k nn (2.17)2.2追赶法解题原理 (1) 计算{}i β的递推公式111/c b β=,/(),2,3,,1i i i i i c b a i n ββ=-=-.(2) 解Ly=f111/y f b =, 1()/(),2,3,,;i i i i i i i y f a y b a i n β-=--=(3) 解Ux=y11,,1,2,2,1n n i i i i x y x y x i n n β-+==-=--.我们将计算系数121n βββ-→→→及12n y y y →→→的过程称为追的过程.将计算方程组的解11n n x x x -→→→的过程称为追赶的过程.3.对追赶法的MATLAB 求解3.1实验程序function x=chase(a,b,c,f) %定义函数chase n=length(b); if n-1==length(a) for i=n-1:-1:1 a(i+1)=a(i); end end%将a 设置为n 维向量 c(1)=c(1)/b(1);f(1)=f(1)/b(1); for i=2:n-1b(i)=b(i)-a(i)*f(i-1)/b(i); c(i)=c(i)/b(i);f(i)=(f(i)-a(i)*f(i-1))/b(i); endf(n)=(f(n)-a(n)*f(n-1))/(b(n)-a(n)*c(n-1)); for i=n-1:-1:1f(i)=f(i)-c(i)*f(i+1);endx=f;3.2 应用举例例1 用追赶法解三对角方程组设270329045A⎛⎫⎪= ⎪⎪⎝⎭567f⎛⎫⎪= ⎪⎪⎝⎭求解线性方程组Ax=F解:(1)在Matlab中编写一个名为chase.m的M文件,依次输入数据如下;>> A=[2 7 0;3 2 9;0 4 5];>> a=[3 4];>> b=[2 2 5];>>c=[7 9];>> f=[5 6 7];>> x=chase (a,b,c,f)x= -3.0140 1.5754 0.1397得到输出结果x= -3.0140 1.5754 0.1397即为原线性方程组的解.结果验证:>>A*x’ans=5.00006.00007.0000得到ans=f 即结果正确4. 与高斯消去法的精度比较事实上,追赶法的求解过程就是将系数矩阵分解两个简单的二对角线矩阵,从而归结为求解两个简单三角形方程组的过程. 高斯消去法是求解线性方程组的最基本方法之一.Gauss 消去法是针对一般的线性方程组,与线性代数中的初等变换解线性方程组方法类似.追赶法只是针对系数矩阵为三对角阵的方程组,因此是一种特殊的方程组.此方法效率较高,不过不适用于一般的线性方程组.下面对追赶法与高斯消去法计算量进行比较.例2 编写高斯消去法解线性方程组的Matlab 程序,并求解下面的线性方程组1231231235282821361x x x x x x x x x ++=⎧⎪+-=⎨⎪--=⎩ 解:编程如下function [RA,RB,n,X]=gaus(A,b) B=[A b];n=length(b);RA=rank(A); RB=rank(B);zhica=RB-RA; if zhica>0,disp('请注意:因为RA~=RB,所以此方程组无解.') return end if RA==RBif RA==ndisp('请注意:因为RA=RB=n,所以此方程组有唯一解.')X=zeros(n,1); C=zeros(1,n+1);for p= 1:n-1for k=p+1:nm= B(k,p)/ B(p,p); B(k,p:n+1)= B(k,p:n+1)-m* B(p,p:n+1);endendb=B(1:n,n+1);A=B(1:n,1:n); X(n)=b(n)/A(n,n);for q=n-1:-1:1X(q)=(b(q)-sum(A(q,q+1:n)*X(q+1:n)))/A(q,q);endelse disp('请注意:因为RA=RB<n,所以此方程组有无穷多解.')endend用Matlab调用程序:在Matlab中编写一个名为gauss.m的M文件,依次输入数据如下:A=[5,2,1;2,8,-3;1,-3,-6];b=[8;21;1]; [RA,RB,n,X]=gauss(A,b)运行结果RA =3RB =3n = 3X =1 2 -1ans = 1.00002.0000-1.0000综上所述,追赶法的原理和高斯消去法相同,但考虑到方程组的特点,计算时会把大量零元素撇开,从而大大节省计算量.追赶法实际上就是把高斯消去法用到求解三对角方程组上的结果.这时由于A特别简单,因此使得求解的计算公式非常简单,而且计算量仅为2(n-1)+(n-1)+[1+2(n-1)]=5n-4次乘除法,而另外解一个方程组Ax=f仅增加3n-2次乘除运算,易见追赶法的计算量是比较小的.课程设计总结首先,这次的课程设计让我把以前学习到的非线性方程的求解知识得到巩固和进一步的提高,对已有知识有了更进一步的理解和认识.同时也发现了自身有许多的不足之处.再者,在这次课程设计中碰到了很多的问题,自己通过上网,查阅了相关的书籍以及和我组其他的成员讨论的方式来解决的.另外编写程序是我的一大弱项,在以前做实验时总是按照老师给出的程序不加思考的敲进去,得出来结果就行,但在做这次课程设计时,需要自己真正的去理解和编程.通过这次课程设计,我从中学会了很多,也发现自己真的还有很多不足以及很多东西需要去学习.所以在以后的生活学习中要不断的扩大自己的视野,多学习一些与专业有关的知识,不能只满足于课本上的知识.所以在完成本专业的基础上,要不断涉猎,完善自我,希望自己在以后的课程中会得到更好的锻练.总的来说这次课程设计还是有很多的收获的,并且特别感谢我们组的成员在做课程设计的过程中对我的帮助.参考文献[1] 冯国忱,黄明游.数值分析(上册)[M].北京:高等教育出版社,2007.[2] 宋叶志,贾东志.MATLAB数值分析与应用[M].北京:机械工业出版社,2009.[3] 张德丰.MATLAB数值分析与应用[M].北京:国防工业出版社,2007.[4] 周煦,计算机数值方法及程序设计,机械工业出版社,2004.10.01.[5] 张静,杨文平,石傅强等,MATLAB程序设计与实例应用,中国铁道出版社,2003.Welcome To Download !!!欢迎您的下载,资料仅供参考!。

ch2.3平方根法和追赶法

ch2.3平方根法和追赶法

则线性方程组可化为两个三角形方程组
Ly b T 1 L x D y
方程组求解公式
y1 b1 k 1 yk bk lkj y j j 1 yn xn d n n yk xk d l jk x j j k 1 k
k 2, , n
k n 1, ,1
举例
二、追赶法 • 追赶法仍然保持LU分解特性,它是一种特殊的 LU分解。充分利用了系数矩阵的特点,而且使 之分解更简单,得到三对角线性方程组的快速 解法。
三对角线性方程组:
b1 x1 c1 x2 a x b x c x 2 1 2 2 2 3 an1 xn1 bn1 xn1 cn1 xn an xn1 bn xn d1 d2 d n 1 dn
x k ( yk ck x k 1 ) uk

“赶”的过程
k n 1 , n 2 , ,1
举例
(**)
若记d (d1 , d 2 , , d n )T ,则三对角方程的矩阵 表示 Ax d, 当 A LU时,可由Ly d及Ux y解出。
三对角矩阵计算公式为 :
x n yn un
y1 d 1 “追”的过程 y k d k l k y k 1 k 2 ,3 , , n
d1ln1l31 d 2ln 2l32 d 3ln3
d1l21ln1 d 2ln 2 d1ln1l31 d 2ln 2l32 d 3ln3 ... ... d1ln1
关于对角矩阵D对角元素及下三角矩阵元素的计算 公式如下(可以课后自己推导)
2 d k akk lkm d m (k 1,2,, n) m 1 k 1
  1. 1、下载文档前请自行甄别文档内容的完整性,平台不提供额外的编辑、内容补充、找答案等附加服务。
  2. 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
  3. 3、如文档侵犯您的权益,请联系客服反馈,我们会尽快为您处理(人工客服工作时间:9:00-18:30)。

x1 2x2
1
2 1
1
A
1
2
1 2 1
x2 2x3 x4 0
x3 2x4 1
1 2
u1 b1, lk ak uk1 , uk bk lk ck1
y1
d1 ,
yk dk lk yk1
x
n
yn
un,xk ( yk ck xk1 )
uk
同解三角方程组为UX=Y,即
x3 1, x2 1 (2)1 1
u32 0 1 2 2, l32 2
x1 4 2111 1
d3 14 11 (2)(2) 9
方程组的解为x1 1, x2 1, x3 1
典型例题
练习 用改进的平方根法解方程组
4 1 1 0 x1 7

1
1 0
3 1 0
1 5 2
2 78 / 25
6 156
/ 11 / 25
典型例题
例5.5 用改进的平方根法解方程组
4 1 1 0 x1 7

1
1 0
3 1 0
1 5 2
0
2 4
x2 x3 x4
8
4 6
方 根
系数矩阵为对称正定矩阵

4
1
1
0
7
1/4 1/4
0
11/ 4 3 / 11
0
3/4 50 / 11 11/ 25
b1
c1
0
A
an1
bn1
cn1
an bn
bi ai ci aici 0(i 2,3, , n 1) (*)
追 赶 法
bn
an
0
则它可分解为
1 l2 1
u1 c1
O
u2 c2 O
A LU
l3 O 1
O
cn1
ln 1
un
分解法的实质是高斯消元 法、LU分解法的应用。事实 上,将系数矩阵A做LU分解:
程序设计
P5={I0[[1]],I0[[2]],{0,0,1/A4[[3,3]],0},I0[[4]]};
A={{-2,1,0,0},{1,-2,0,0},{0,1,-2,1},{0,0,1,-2}}; A5=P5.A4;
b={1,1,0,-1};
MatrixForm[%];
I0=IdentityMatrix[4];
0 0
2 1
Байду номын сангаас0 0
1 2
1 3
2
0
1
2 1
2 0
0
1 3
2 2 3
0
0 0
2 1
2
0 0
1 3
2
1 3
2
1
1 2
1
1 1* 0.3333
x4 3 0.3333 x3 2 0.3333
3 0*(0.3333)
x2 2
3
1
2
1 1*(1) x1 2 1
程序设计
A LDLT 其 中L为 单 位 下 三 角 阵 ,D为 对 角 阵 。
比较
平方根法(Cholesky):A LLT 改进平方根法:A L1 DL1T
优点:1)不必选主元 2)算法稳定
改进的平方根法
A LDLT

1 0 0
l
21
1
0
l31
l32
1
ln1 ln2 ln3
0
0
d1
“赶”的过程
典型例题
例5.5
2x1 x2
x1 2x2
x2 2x3 x4
x3 2x4
1
1 A 0 1
2 1
1
2
1
2 1
1 2
u1 b1, lk ak uk1 , uk bk lkck1
y1
d1 ,
yk dk lk yk1
xn
yn
un,xk ( yk ck xk1)
P6={I0[[1]],I0[[2]],I0[[3]],{0,0,-A5[[4,3]],1}};
MatrixForm[%];
A6=P6.A5;
P1={{1/A[[1,1]],0,0,0},I0[[2]],I0[[3]],I0[[4]]}; MatrixForm[%];
A1=P1.A;
P7={I0[[1]],I0[[2]],I0[[3]],{0,0,0,1/A6[[4,4]]}};

y1 b1
yk
bk
k 1
lkj y j
(k 2, , n)
j1
在平方根法中,需进行 n 次开方运算, 为了避免这一点,在更多情况下是将A作
xn yn / dn
xk
yk dk
n
l jk x j
j k 1
(k
n 1,
,1)
A LDLT 分解。 这一方法称为改进的平方根法。
典型例题
Table[Switch[i-j,-1,a,0,b,1,c,_,0],{i,6},{j,6}]; MatrixForm[%]

b1 c1
a2
b2
c2
对应的系数矩阵A
an1
bn1
cn1
an bn
分解定理
分 解 定 理 Add Your Title
b1 c1
a2
b2
c2
设矩阵A满足下列条件:
0
1
0 d2
0
1 0
l21 1
l31 l32
0 0 1
d
n
0
0
0
ln1
ln
2
ln
3
dk
akk
k 1
lkm 2dm
m1
lik
aik
k 1
l im d m l km
m1
dk
1
i k 1, k 2, , n
方 根
求解 Ly b, LT x D 1 y计算公式为
b
ln 1
un
这种把三对角方程组
追 赶 法
三 对 角 矩 阵 计 算 公 式 为:
y1 d1
“ 追 ” 的 过 程
yk dk lk yk1
k 2,3, , n
xn yn un
的解用递推公式表示出来 的方法,被形象化地叫做 追赶法。
xk ( yk ck xk1 ) uk
k n 1, n 2, ,1
uk

u1 2 y1 1
赶 法
l2
1 2
u2
3 2
y2
3 2
x4
1 3
0.3333
x3
1
1* 0.3333 2
0.3333
l3
2 3
l4
1 2
u3 2
u4
3 2
y3 1
y4
1 2
3 0*(0.3333)
x2 2
3
1
2
1 1*(1) x1 2 1
典型例题
例5.5
2x1 x2
A3=P3.A2;
MatrixForm[%]
MatrixForm[%];
L.U;
P4={I0[[1]],I0[[2]],{0,-A3[[3,2]],1,0},I0[[4]]}; MatrixForm[%]
答案:
A4=P4.A3; MatrixForm[%];
LinearSolve[L,b]; y=% LinearSolve[U,y]
三对角矩阵 分解定理 “追赶”过程 典型例题 通用程序设计
三对角矩阵
三对角线性方程组:
程序设计
b1x1 c1x2
a2 x1 b2 x2 c2 x3
d1 d2
Clear[i,j,a,b,c]
追 赶
an1 xn1 bn1 xn1 cn1 xn dn1 an xn1 bn xn dn
例5.4
用改进的平方根法解方程组
x1 2x2 x3 4
2x1 5x2 7
x1
14x3 15
求解 Ly b, LT x D 1计y 算公式为
dlxyyikn1kk
k 1
b1akk lkm 2d m
k 1 m 1
bk
lkjky1j
a
j 1
ik
l
im
(k
dml
yn / dn m1
0
3/4 5 2
0 25 / 4
2 4
4 6
4
1 1 0 7
1/4 1/4
0
11/ 4 3 / 11
0
3/4 50 / 11 11/ 25
0 25 / 4
2 4
6/ 6
11
4
1/4 1/4
0
1 11/ 4 3 / 11
0
1 3/4 50 / 11 11/ 25
0
7
0 25 / 4
L Table 0, i, n , j, n ;
U Table 0, i, n , j, n ;
For k 1, k n, k , For j k, j n, j , U k, j A k, j
For i k 1, i n, i , L i, k
A i, k Sum L i, t
L L IdentityMatrix n ;
平 方 (对称正定矩阵的LU分解形式更加简单, 根 平方根法就是针对正定矩阵的LU分解法) 法 定理 设 A 是对称正定矩阵,则存在惟一的
非奇异下三角阵L,使得
A LLT
且 L 的对角元素皆为正数。
相关文档
最新文档