数学系毕业论文.docx

东北师范大学本科论文 M-矩阵西尔维基斯特方精确解的数值模拟 摘要本文主要讨论扰动理论及其关于M-矩阵西尔维基特方程相关分量的精确解，其中的对角元都是正的，非对角元都是非正的，C是非负矩阵是非奇异M-矩阵，能够证明A,B,C分量小的扰动可以引起解X分量小的相关误差.因此，相对较大的分量X中较小的分量不会引起跟更大的误差，它对一般的西尔维基特方程是不成立的.接下来我们将对3种数值方法进行小的却十分重要的修正使他们可以用来精确计算X.关键词：M-矩阵；西尔维基特方程；扰动；相关误差AbstractThis paper is concerned with a relative perturbation theory and its entrywise relatively accurate numerical solutions of an M-matrix Sylvester equationby which we mean both A and B have positive diagonal entries and nonposive off-diagonal entries and is a nonsingular M-matrix ,and C is entrywise nonnegative. It is proved that small relative perturbations to the entries of A,B,and C introduce small relative errors to the entries of the solution X .Thus the smaller entries of X do not suffer bigger relative errors than its larger entries , unlikely the existing perturbation theory for (general)Sylvester equations . We then discuss some minor but crucial implementation changes to three existing numerical methods so that they can be used to compute X as accurately as the input data deserve. Keywords :M-matrix ;the Sylvester equation ;perturbation ;relative accurate error目录第一章背景介绍- 4 -1.1M-矩阵及西尔维基特方程介绍- 4 -1.2研究MSE的意义及目标- 4 -1.3 论文的结构介绍和相关规定- 5 -第二章 M-矩阵的逆阵- 5 -2.1 扰动理论- 5 -2.2 如何精确计算可逆矩阵- 7 -第三章 分量扰动分析- 8 -3.1 主要结果- 8 -第四章 西尔维基特方程的算法及数值案例- 11 -4.1 直接法- 13 -4.2 固定点迭代法- 13 -4.3 Smith 算法- 14 -4.4 数值案例- 16 -第五章 结论- 18 -参 考 文 献- 20 -致 谢- 22 -附录- 23 - 第一章 背景介绍1.1 M-矩阵及西尔维基特方程介绍 对于实数域上阶实矩阵,若存在阶非负矩阵使,则称是M-矩阵，其中表示阶单位矩阵,表示求一个矩阵的谱半径。 若,称是非奇异M-矩阵;若,称是奇异的M-矩阵.必然地，一个非奇异M-矩阵有非正的非对角元和非负的对角元.对于一个非奇异或不可约的M-矩阵，它的对角元一定是正的.本文我们主要关心的是下面的西尔维基特方程 (1.1)有正的对角元和非正的非对角元，是一个非奇异M-矩阵, ,是非负矩阵，是矩阵（或向量）的克罗内克积.我们称（1.1）这种类型的西尔维基特方程为M -矩阵西尔维基特方程（简称MSE）. 1.2研究MSE的意义及目标MSE经常出现现M-矩阵代数里卡蒂方程的迭代方法中.M-矩阵代数里卡蒂方程 （简称MARE）定义如下: （1.2）其中都是矩阵.它们的阶数有下面矩阵分割确定 ， (1.3)它是一个非奇异或不可约奇异的M-矩阵.在(1.3)中令，(1.3)变成MSE（1.1）.但是根据MSE的定义，并不是所有的MSES都有这种形式，例如那些和不全是M-矩阵的MSES就不会有（1.3）这样的形式.MARE(1.3)有唯一最小非负解是显而易见的14,15.MARE(1.3)一个非常重要的应用是随机流体模型23,24，在随机流体模型中解的每一个分量都有一个物理意义，精确计算解的每个微小分元变得十分重要.为了保证这点，我们就需要解决出现在某些迭代方法中的MSE . 本文我们的第一个目标就是针对MSE（1.1）提出分量是扰动分析.特别地我们要确定由系数矩阵A,B,C分量扰动引起的解X每个分量误差范围结果显示每个解及解的分量（无论它们多小） 都决定与系数矩阵分量精度可对比的相关精度.我们的第二个目标是介绍一些能够精确计算MSE解的算法，他们包括有（无论它们多么小）都决定与系数矩阵分量精度可对比的相关精度.我们的第二个目标是介绍一些 能够精确计算MSE解的算法.1.3论文的结构介绍和相关规定论文的第二部分主要讨论了M-矩阵的可逆阵的扰动理论以及如何通过理论精确计算可逆阵.利用第二部分和第三部分的结果建立MSE的一个相关扰动理论.第四部分解释了3种经典方法：GTH-like算法，固定点迭代法以及Smith 算法.第五部分给出了一些数值实例.第六部分陈述了论文的结论.在随后的文章中我们做了一下规定:表示n阶单位矩阵的第j个分量;表示所有分元都是1的向量;；表示一个矩阵,它的分元为;当时表示与矩阵Z对角元相同的对角矩阵：令;表示将矩阵的列向量组成一个向量，称为矩阵的外部运算；指的是所有的类似定义等,特别即X是非负的.是一个矩阵，它的(i,j)分元是并规定；对任意矩阵或向量，用带浪线的相同字母表示它们的扰动.第二章 M-矩阵的逆阵在这部分，我们会介绍关于M-矩阵逆阵的一些结果，它们会在后面分量式扰动分析以及MSE(1.1)的算法中发挥作用。下面的结果非常出名7.定理 2.1 有非正的非对角元，下面几条是等价的.(1)是一个非奇异M-矩阵；(2);(3)存在一个，使得(4)的所有特征值都有正的实部.实数域上阶方阵,如果存在一个置换矩阵，使得，其中是方阵，则称是不可约的；否则，称是可约的.在接下来的小部分中规定是一个非奇异M-矩阵，它的扰动矩阵为，令 (2.1)2.1扰动理论定理2.2 2 如果存在和，使得 (2.2)则是一个非奇异M-矩阵，且有 (2.3)附注2.1（1）在（2.2）中(2)（2.2）中可由 替换.由（2.2）可知A的对角元是确定的.因为 (2.4)其中，因此（2.2）强于. 定理2.3 假设,如果， 则是非奇异M-矩阵 且有 （2.5）证明：假设在（2.1）中是不可约的，令是贝龙特征向量，即有 我们知道,易知，其中 ，则由于及定理2.2知是非奇异M-矩阵，因此也是非奇异M-矩阵，同时有现在考虑是可约的情况.对于充分小的，是不可约M-矩阵，应用刚刚证明的结果，再令即证得结论.附注2.2（2.5）中不等式可推出，对于充分小的有 (2.6)附注2.3 在定理2.3条件下，通常的一阶误差分析如下， 展开 可得到,由，则 （2.7）= = (2.8)又由于 (2.9)对比(2.5)， (2.9)仅是一阶误差范围，(2.9)没有指出为什么和什么时候是微小的，不同于（2.5）和成正比，（2.9）式在实际操作中容易实现,只需计算足够的就可以得到解中每个分量至少一位正确的小数.尽管有时这个序列 收敛缓慢，但它却仍是非常有价值的. 2.2 如何精确计算可逆矩阵根据1,A的三元数表示具有数值上的优势，其中 .A的对角元可以通过（2.4）得到，下面我们将探讨识别一个M-矩阵以及它有效的参数化的三元数 表示 其中 (2.10)需要指出的是在数值上并不精确的等于，但它是分量精确近似值.当矩阵的三元数表示已知，由GTH-like算法1是向前稳定的，它计算出的，它计算出的具有分量相对精度；当矩阵三元数表示未知，则我们需要找到一个正向量使 ,它可以通过求解中的实现.理论上，除非是几乎奇异的，相对于1n来说剩余的对的计算影响是微小的，这意味着 实际上对于几乎奇异的情况= 也是成立的.假设这种情况出现，不难证明存在使得，则有,其中是机器舍入精度.分解.然后我们发现.应用定理2.3,有 (2.11)其中. 但是因为是不知道的,在实际计算中我们不得不用M-矩阵的逆阵代替它.由于GTH-like算法是向前稳定的，则计算得到的不同于，且有.定理2.2说明，则 ，其中是低阶多项式.这暗示作为的近似值保证了定理2.3要求的精度.当或不再适用通过来计算时，是几乎奇异的.解决是反迭代非常基础的一步，它意味着非常接近的贝龙特征向量.有两种子情况：是不可约的，则的贝龙特征向量分量都是正的，因此我们需要迭代更多步 直到和,解 得;令其中.表示求一个矩阵的无穷范数；是可约的. 更多的迭代也没有帮助，因为贝龙特征向量有的分量可能是零.对于这种情况我们首先找一个置换矩阵P使得其中是非负M-矩阵，它能够精确计算可以直接应用上面不可约的结果，因为是分块上三角矩阵且它的对角矩块是 ，而非对角块没有单独的减法可以很容易计算出来.第三章 分量扰动分析在这部分，,按如下方式分解： (3.1a) (3.1b)且有和有非负的对角元， 和 是非负的. 令 (3.2)3.1 主要结果MSE(1.1)主要的方式扰动结果是定理3.1和定理3.2 .定理3.1假设和都是非奇异M-矩阵,而定理3.2没有要这个但是有条件暗示是一个非奇异M-矩阵.MSE(1.1)等价于,是一个非奇异M-矩阵，关于一个非奇异M-矩阵逆阵的分量扰动结果可以直接应用到的分量相对扰动范围上.事实上这种情况就是定理3.1.对于和空间维数特别大的时候，很难估计，为此我们将采取一个不直接的方法建立涉及不涉的范围. (3.3)其中对于充分小的我们近似有,这部分主要的结果有下面的形式： (3.4a)或近似 (3.4b)其中之后会详细说明.定理3.1 假设和是非奇异M-矩阵，存在,使得 (3.5a) (3.5b) (3.5c)则和也非奇异M-矩阵，（3.4）成立，且.证明：MSE等价地写成和.因为和是非奇异M-矩阵，则是非奇异M-矩阵.3.5a和3.5b中不等式保证了和也是非奇异M-矩阵；因此也是非奇异M-矩阵.由3.5和3.6有类似地因此 另一方面,对于和相同位置的非对角元必是下面3种情况之一:0和0，和，和，其中因此由（3.5a）和（3.5b）有.由于满足定理2.2的条件，因此同理附注3.1在（3.5a）和（3.5b）中最后面的不等式可以由替代，结论仍成立.实际上求解没有必要是非奇异M-矩阵，只要求是非奇异M-矩阵即可下面引理要求的条件保证了的非奇异性.引理3.1 如果 (3.6) 和 (3.7)其中有一个成立，则是非奇异M-矩阵. 证明：首先考虑和都是不可约的，我们只需考虑3.6式，令分别是贝龙特征向量即若易知P是一个非奇异M-矩阵；假设，我们有注意,则 （3.8）由定理2.1知是一个非奇异M-矩阵.由3.9及7，P.28知 .设是具有最小实部的特征值,则有18,189页问题19知 (3.9)现在考虑或是可约的，令,类似定义,我们有 令，由特征值的收敛性，因此3.10对于可约的情况也成立.综上是一个非奇异M-矩阵.定理3.2 在引理3.1条件下，假设, (3.10) 令 如果则（3.4）成立.证明：引理3.1说明是一个非奇异M-矩阵，假设和都是不可约的，设分别是的贝龙特征向量，所以3.8成立. (3.13)(3.14) (3.15)由定理条件知不等式的右边均是正的.现在我们观察的分量比，由(3.13)和(3.14) （3.16） （3.17）令,对（3.16）右侧关于t求导，又由（3.6）,有 则（3.16）中右侧是一个关于t的递增函数，（3.17）不等号成立.令因此 ，同理可证，易证应用定理2.2完成证明.若A和B是可约的，对应用刚刚证明的结果，再令，得出结论.附注3.2 类似附注2.3，常用的一阶误差分析也可以执行.在P是一个非奇异M-矩阵及（3.11）成立假设下，,其中Z可由A,B,C,X等替换，表示对Z的扰动.把它们替换中，则有定义线性算子,其中它的矩阵表示为且,因此 (3.18)因为和有相同的零或非零模式，所以我们 (3.19)其中且是由下面的式子定义的， , (3.20)第四章 西尔维基特方程的算法及数值案例关于MSE (1.1)的扰动理论在第三部分中包含两种情况：A和B是非奇异M-矩阵且;和 中至少有一个是非奇异M-矩阵（另一个可能不是M-矩阵）且是一个非奇异M-矩阵，.实际上这两种情况放到一起就能包含我们所以定义的MSE的所有可能结果.下面我们来证明一下这个结论.证明：我们注意到的特征值为和和的可能结果，其中和分别是和的特征值,是一个非奇异M-矩阵，由定理2.1知，其中表示取一个复数的实部，则 或，其中表示一个矩阵所有特征值组成的集合. 因此和至少有一个是非奇异M-矩阵.此外,令 （4.1）则和的所有特征值都有正的实部，因此和都是非奇异M-矩阵.不是一般性，我们可以假设和都是不可约的；否则，等式（1.1）可以被分解成比较小的西尔维基特方程序列，这些比较小的方程都有（1.1）这样的形式且其系数矩阵都是不可约的.事实上，令是置换矩阵，存在向量，使得即其中，所有的分块都是不可约的，是非奇异M-矩阵，(). 把分成个分块，每个子块则 有 分别取等式两边i,j块有 ,它暗示 的列快可以从最后一列依次向前计算，每一列块又可以从最低端依次向上 计算.当A和B都是矩阵时，我们考虑这种情况（4.2）如果向量,都是预先给出的，则显然成立；如果向量 和给出的是一般矩阵形式，和的近似三元数表示可以通过2.2部分描述的方法计算出.4.1 直接法直接法是以高斯消去法为基础解决MSE等价问题,由（4.2）， （4.3）因此应用GTH-like算法证明结论.当A和B中有一个不是M-矩阵，但是M-矩阵，我们不会有（4.2）式，而P的三元数表示仍可找到.假设，是一个非奇异M-矩阵，我们利用2.2部分的想法找到的贝龙向量且假设是不可约的.因为不是M-矩阵，所以的分量有可能是负的.即便如此，若我们任由（4.3）；如果 我们重新找的三元数表示，其中是的贝龙向量且.如果定理3.2条件满足，则，则我们又有(4.3).直接法计算量大为，即使对适度的也很困难,但对较小的它是一个非常理想的计算精确解的方法.4.2 固定点迭代法对于的任一分解 ， (4.4)给出的一个分解 , , (4.5)对于MSE(1)相应的一个迭代法：, (4.6)为了方便称（4.5）为正则分解，即 在这样的分解下(4.6)变得很容易解决.下面给出五种明显的策略： , (4.7a) , (4.7b) , (4.7c) , (4.7d) , (4.7e) 其中和是MATLAB-like中符号，分别表示求下三角矩阵和上三角矩阵., . 因为是非奇异M-矩阵，所以 ，则（4.6）中迭代法是收敛的，而且 （4.8）虽然西尔维基特方程本身很难执行，但对于（4.7）中每一个分解执行（4.6）是容易的. 对于（4.7b）-(4.7e), 和都是三角矩阵且（4.6）可以被分解成一个三角矩阵序列，（4.6）的简单实现总可以给出但它不再保持（4.8）中的单调性.下面给出一个比较好的方法，令算法 4.1（固定点迭代法）第一步：解方程 得;第二步：;第三步：直到收敛；第四步：解方程 得；第五步：；第六步：停止.对于（4.7）中每一种分解，算法4.1都能保证产生一个的线性收敛单调序列，它的收敛速率为.在我们数值实验中，如果有，则停止迭代.为了以后方便，我们用定义在（4.7a）-(4.7e)各自分解下的算法 Smith 算法这个解决（1.1）的迭代方法来自于Smith 算法24,它能计算出有（4.2)形式的A和B, A和B中至少有一个是非奇异M-矩阵，而另一个只需是M-矩阵就可以.对任意一个数，有.如果，则和都是非奇异M-矩阵，从而,其中 此外，对由于非奇异M-矩阵所有特征值实部都是正的，则（1.1）的解X承认下面级数展开式由24知，它 是快速得到近似的， ，实际上.为了我们的目的，我们将找一个使得 （4.9）所有的和都可以被高精度地算出来.不等式（4.9）能确保下面式子成立，且近似有 （4.10）算法 4.2(Smith算法)第一步：找到一个，这样在计算所有的和 时不会有巨量消失；第二步：第三步：用GTH-like算法计算第四步：第五步：直到收敛；第六步： 第七步： 第八步： 停止. 附注4.1 1用三元数表示和是为了应用GTH-like 算法.不同于解决（1.1），那里没有要求一定要有三元数表示，本算法 在利用三元数表示时效果会更好; 2 前四步都没有涉及单独的运算，因此所有的计算都是向前稳定的;3 需要解释什么时候停止迭代使最后一个具有的分量对精度可以作为一个近似值.我们借鉴Prof.W.Kahan的想法，考虑一个非负单调上升序列即有 （4.11）令如果当是递减的，且，则 （4.12）为了得到（4.12）式，令.对，有则，从而有 .为了计算有相对精度为的，只需要 （4.13）就可以停止迭代了，这个迭代标准可以在第五步-第八步循环圈中很容易得到.4.3 数值案例在这部分，我们主要给出一个数值案例检验我们第四部分给出的固定点迭代方法计算分量相关精确解的能力.接下来，我们将用两个误差尺度去估计算的解的精度：标准化的残差（简称 (5.1)这是估计精度的一般测度因为它能够很容易地被计算出来.分量相关误差（简称）， （5.2） 看下面的例子采用五种策略并比较这两种误差尺度例1 令其中，如果，则是一个M-矩阵.事实上我们已经证明了取很多值时观察到都有类似的结果，接下来我们仅列举的情景. 这三个图分别是=10, =20,=40时五种固定点迭代的ERErr曲线，观察发现图像发现只出现三条曲线，因为FPa曲线和FPd 曲线重合,FPb曲线和FPC曲线重合,三条曲线中FPa(FPd)拟合效果最好.随着n增加FPa(FPd)和FPe变化不大，但FPb(FPC)有上下波动，随着迭代步骤增加效果会比FPe好，却始终没有FPa(FPd)效果好.而下面三个图分别是=10, =20, =30时五种固定点迭代的NRes曲线观察发现图像发现只出现三条曲线，因为FPa曲线和FPd 曲线重合,FPb曲线和FPC曲线重合,三条曲线中FPe拟合效果最差. =10时FPb(FPc)在迭代40次左右突然变大,在50次左右就变为1（由于FPb(FPc)的谱半径不严格小于1）.随着n增加大，FPb(FPC)效果会比FPa(FPd)差，且随着迭代步骤增加FPb(FPC)效果会比FPa(FPd)差. 第五章 结论如果是M-矩阵，是非奇异的且, 我们能更自然地称（1.1）为一个MSE，我们关于MSE的目的会更加显著，但实际上（1.1）等价于这个自然的结果因为对于任意一个数和（4.1）给出的来说,是一样的.我们已经提出了一个关于MSE(1.1)的分量扰动分析,能够证明A,B,C分量相关扰动只能够引起解X分量的相关改变，无论量级多大.我们认为范围中线性部分是最好的近似结果.一般的一阶误差分析也可以很容易实现，正如我们在附注3.2中概述，目的是得到一个明显的且容易实现的一阶误差界.但是我们新的关于对解敏感度的影响的分析给出一个深刻的误差范围. 根据我们的分析，对三种方法进行小的但十分重要的修正后它们能够计算出符合分量相关精度要求的解，我们的数值结果能够说明算法的分析及精度理论是对的.本论文的不足之处是没有对Smith算法进行相应的数值模拟，Smith算法初值是由GTH-like算法给出十分复杂，GTH-like及Smith算法将作为后续学习的重点.参 考 文 献1.1Alfa,AS,XueJ.,Ye,Q.:Accurate computation of the smallest eigenvalue of a dia- gonally M -matrix.Math.Comput.71,217-236(2002)2.Alfa,A.S.,Xue,J.,Ye,Q:Entry wise perturbation theory for diagonally dominant M- matrices with applications .Numer.Math.90(3),401-414(2002)3.American National Standards Institute and Institute of Electrical and Electronic Engineers; IEEE standard for binary floatingpoint arithmetic .ANSI/IEEE Stan- dard ,Std 754-1985,New York (1985)3.American Nation Standards Institute and Institute of Electrical and Electrical Engineers :IEEE standard for radix independent floating-point arithmetic .ANSI/IEEE Standard ,Std 854-1987,New York（1987） 5.Bartels,R.H.,Stewart G.W.: Algorithm 432:The solution of the matrix equation. Commun .ACM 8,820-826(1972)6.Benner,P.,Li,R.-C,Truhar,N.:On ADI method for Sylvester equations. J. Comput. APPL .Math. 2334(4),1035-1045(2009)7.Berman,A.,Plemmons,R.J.: Nonnegative Matrices in the Mathematical Sciences .SIAM ,Philadelphia, 1994.This SIAM edition is a corrected reproduction of the work first published in 1979 by Academic Press, San Diego ,CA8.Elsner,L .,Koltracht ,I .:Mixed ,componentwise ,and structured condition numb- ers .SIAM J.Matrix Anal.Appl.14(3) ,688-704(1991)9.Goldberg,D.:What every computer scientist should know about floating-point ar- ithmetic .ACM Comput .Surv .23(1),5-47(1991)10.Gohberg, I .,Koltracht ,I. :Mixed componentwise ,and structured condition num- bers. SIAM J. Matrix Anal. Appl .14(3),688-704(1993)11.Golub,G.H.,Nash,S.,van Loan, C.F.: Hessenberg Schur method for the problem.IEEE Trans .Autom .Control AC-24,909-913(1979)12.Grassmann,W.K.,Taksar,M.J.,Heyman,D.P.:Regenerative analysis and stady state distributions for Markov chains .Oper. Res. 33,1107-1116(1985)13.Guo,C.,Higham,N.:Iterative solution of a nonsymmetric algebraic Riccati equation SIAMJ. Matrix Anal.APPL.29,396-412(2007)14.Guo,C.-H.:Nonsymmetric algebraic Riccati equation and Wiener- Hopffactorization for M-matrices .SIAM J. Matrix Anal.Appl.23,225-242(2001)15.Higham,N.J.:Accuracy and Stability of Numerical Algorithms.2nd edn.SIAM ,Philadephia (2002)16.Horn,R.A.,Johnson,C.R.:Topicsin Matrix Analysis .Cambridge University Press ,Cambridge(1991)17.Juang,J.:Existence of algebraic matrix Riccati equations arising in transport theory .Linear Algebra Appl.230,89-100(1995)18.Li,R.-C.:Solving secular equations stably and efficiently. Technical Report UCB/CSD-94-851,Computer Science Division, Department of EECS ,University of California at Berkeley(1993)19.Ramaswami, V .:Matrix analytic methods for stochastic fluid flows. In :Key, P . Smith, D.(eds.) Teletraffic Engineering in a Competitive World ,vol.3a of Teletraffic Science and Engineering. Elsevier Science ,Amsterdam, pp.1019-1030(1999)20.Smith,R.A.:Matrix equation.SIAM J .Appl .Math. 16(1),198-201(1968)21.Stoer,J.,Bulirsch,R.:Introduction to Numerical Analysis.2nd edn.Springer- Ver lag ,Berlin(1992)22.Varga, R.S .:Matrix Iterative Analysis .Englewood Cliffs, NJ(1962)23.Virnik, E .:Analysis. PhD thesis ,Technischen Universitat Berlin(2008)24.Xue,J.:Computing the smallest eigenvalue of an-matrix .SIAM J .Matrix Anal .Appl.17(4),748-762(1996)25.Xue,J.,Jiang, E .:Entrywise relative perturbation theory for nonsingular M-matrices and applications .BIT 35(3),417-427(1995)26.Xue,J.,Xu,S.,Li,R.-C.:Accurate solution of M-matrix algebraic Riccati equations.Numer.Math.doi:10.1007/s00211-011-0421-0 致 谢时间匆匆，四年的大学时光转瞬即逝.回想这四年的学习生活，收获很多。在学院老师和领导的悉心指引下，我的学习成绩有很大提升！在本文的选题、选材及写作的整个过程中，我的导师刁怀安副教授给了我很大的启迪与帮助，论文的修改更是有老师的不少心血.他严谨的工作态度，积极地工作热情时时感染着我，激励着我。他的谆谆教诲必将使我受益一生.感谢我的同学在我查找资料和论文写作过程中给与的帮助和支持. 在论文的写作过程中参阅了许多研究文献，在此也向他们表示感谢! 附录function X,des_0,iteration,nres_0=Point_1(A,B,C,option) %固定点迭代的程序Maxiter=100;e_0=1.0e-1;des_0=;nres_0=;iteration=;switch option case FPa M_1=diag(diag(A);M_2=diag(diag(B); case FPb M_1=tril(A);M_2=triu(B); case FPc M_1=triu(A);M_2=tril(B); case FPd M_1=tril(A);M_2=tril(B); case FPe M_1=triu(A);M_2=triu(B);endK_1=M_1-A;K_2=M_2-B;X_1=lyap(M_1,M_2,C);X=X_1;for i=1:Maxiter C_1=K_1*X_1+X_1*K_2; X_2=lyap(M_1,M_2,C_1); X=X+X_2; temp=X_2-X_1; des=norm(pinv(diag(X_2(:)*(temp(:),inf); nres=norm(A*X+X*B-C,1)/(norm(X,1)*(norm(A,1)+norm(B,1)+norm(C,1); des_0=des_0;des; nres_0=nres_0;nres; iteration=iteration;i; if des=e_0 disp(X); break end X_1=X_2;endend function A,B=struction(n,w) %构造A和B的程序q=(1/3)/w(n-1);m=eye(n-1);p=eye(n);e_n=p(:,n);e_1=p(:,1);U=zeros(n-1,1) m;zeros(1,n-1) 0;A=eye(n)-w*U-q*e_n*e_1;B=A;end function plot_0(n,w) %五种策略关