利用Mathematica软件实现线性方程组求解

利用Mathematica 软件实现线性方程组求解摘要：Mathematica 作为一款独特的数学软件,有强大的数值计算和符号计算能力,具有很强的实用性,而线性方程组在线性代数中具有重要的地位.针对线性方程组的解具有零解、唯一解、无穷解的情况,文章借助Mathematica 数学软件分别给出求解的算法以及如何直接调用内部函数求解线性方程组的方法.文章中主要利用高斯-约当（Gauss-Jordan ）消元法、矩阵的LU 分解法求解线性方程组,并用Mathematica 软件对这两种算法予以实现.同时给出了带有可读性好的求解过程、能用于各种逆矩阵和线性方程组求解的通用程序,利用Mathematica 数学软件实现对求逆矩阵和线性方程组的可读性计算.关键字：Mathematica ;线性方程组;可读性计算;算法引言1.Mathematica 软件功能简介Mathematica 作为一款独特的数学软件,有强大的数值计算和符号计算能力.同时,线性方程组是连接矩阵和向量组的纽带,也是矩阵和向量组的直接应用.在我们的实际工作中出现的大量数学模型,他们最后都可以化为解线性方程组,所以对线性方程组的解的研究是非常重要的.在数学的许多分支中都涉及到线性方程组的问题,我们通常分齐次与非齐次两大类以及未知量个数与方程个数相等或不等对其进行求解,在系数矩阵为方阵时,还可以按照行列式是否为零进行分类.对于线性方程组的求解,其算法是非常清楚的,就是利用矩阵的行初等变换将对应的增广矩阵化为行最简矩阵,然后就可写出通解.但是由于整个过程涉及到大量数字运算,往往会因计算过程中不小心出现一些计算错误而导致错误的结论,有时甚至出现对没解的方程求出有解,或相反.然而,如果利用数学软件Mathematica ,可以将大家从繁琐的数字运算中解脱出来,而将注意力集中在基本概念和基本算法的学习上,从而增加学习兴趣,提高学习效率.笔者从事了一年多大学生创新实验《基于数学软件的高等代数解题实践研究》,通过实践发现,如果能够通过计算机编程让计算机自动求解各种线性方程组,这将是非常快速、方便的.因此本文以Mathematica 数学软件为平台,通过举例、演示与实验来理解线性方程组中的一些抽象概念和理论,并简捷直观地用计算机来解决复杂的线性方程组的求解问题,化简过难、过繁的运算技巧.本文根据Mathematica 软件强大的数值计算无误差的特点和符号计算功能,总结Mathematica 的关于线性方程组计算的功能、命令,总结了求解线性方程组的解的算法,并给出带有可读性好的求解过程、能用于各种线性方程组求解的通用程序,实现线性方程组的可读性计算.由于本章用到线性代数的许多基本知识,有必要对相关的概念和性质做一下说明. 定义1 如两个方程组具有完全相同的解,则称两个方程组为等价方程组.定义2 关于方程组的初等运算有三种类型：1） 交换方程组中两个方程：i j R R ↔,称为交换运算.2） 用一个非零的实数乘以某一方程：i i R R λ→,称为倍法运算.3） 某方程加上另一个方程的倍数：i j i R R R λ+→,称为消法运算.定理1 若一个方程由另一个方程经过有限次初等运算得到的,则这两个方程组等价.定义3 矩阵的初等变换,用i A 表示矩阵A 的第i 行,那么有：1） 交换A 的两行：i j A A ↔,称为交换变换.2） 用一个非零的实数λ乘以A 的某一行：i j A A λ→,称为倍法变换.3） 用A 某行加上另一行的λ倍：i j i A A A λ+→,称为消法表换.定义4 对矩阵A 施行初等变换,得到的矩阵称为相应的初等阵.例如23100100010001001010A A ↔⎛⎫⎛⎫ ⎪ ⎪−−−→ ⎪ ⎪ ⎪ ⎪⎝⎭⎝⎭记作23E2210010001000001001A A λλ⨯↔⎛⎫⎛⎫ ⎪ ⎪−−−−→ ⎪ ⎪ ⎪ ⎪⎝⎭⎝⎭记作()2E λ32310010001001000101A A A λλ+⨯↔⎛⎫⎛⎫ ⎪ ⎪−−−−−→ ⎪ ⎪ ⎪ ⎪⎝⎭⎝⎭记作()23E λ定理2 对矩阵施行初等变换相当于将矩阵左乘相应的初等阵.212+-210210221-101-5001001A A A ⨯→⎛⎫⎛⎫ ⎪ ⎪−−−−−−→ ⎪ ⎪ ⎪ ⎪⎝⎭⎝⎭（）相当于12(-2)E A100102102-21021-1=01-5001001001⎛⎫⎛⎫⎛⎫ ⎪⎪ ⎪ ⎪⎪ ⎪ ⎪⎪ ⎪⎝⎭⎝⎭⎝⎭定理3 对任何n 阶矩阵A,以下性质等价：1） A 存在逆矩阵,即A 是可逆矩阵.2） 0A ≠,即A 是非奇异矩阵.3） A 的行向量组(列向量组)线性无关.4） 方程组AX b =有唯一解.定理4 如果矩阵A 经过一系列的初等变换化为单位阵E,则单位阵E 通过相同的初等变换化为1A -.由此得出用初等变换求矩阵逆的方法下例所示. 例1：求矩阵123133247A ⎛⎫ ⎪= ⎪ ⎪⎝⎭的逆. 解：()(2)(1)(1)(2)(3)(1)13(3)(2)13(1)(3)(3)(1)(1)(2)211231*********33010010-110247001011-201123100100010-110010001-1-110A E +-⨯→+-⨯→+-⨯→+-⨯→+-⨯→⎛⎫⎛⎫ ⎪ ⎪=−−−−−−→−−−−−−→ ⎪ ⎪ ⎪ ⎪⎝⎭⎝⎭⎛⎫ ⎪−−−−−−→ ⎪ ⎪⎝⎭（）（）（）（）（）（）()-161-3-110=01-1-11E A ⎛⎫ ⎪ ⎪ ⎪⎝⎭ 2 Mathematica 软件常用命令及其功能为了便于读者更好地利用Mathematica 软件解决数学问题,首先给出本文程序涉及的Mathematica 数学软件的命令及其功能.Append[expr,elem]——生成一个在表expr 的最后添加元素elem 后的表;Complement[listall, list1,list2,…]——求全集listall 对listi 的差集; Delete[expr,{i}]——删除expr 中第i 个元素后剩下的表;Do[expr,{imax}]——重复执行expr 共imax 次;Table[expr,{imax}]-—生成一个由表达式expr 产生的元素构成的共imax 个元素的表; Do[expr,{i,min,max}]——循环变量i 按步长1从min 取至max,重复执行表达式expr; Do[expr,{i,imin,imax},{j,jmin,jmax},…]——多重循环;Dot[a,b]——矩阵、向量、张量a 与b 的内积;Transpose[list]——对矩阵list 进行转置; IdentityMatrix[n]——自动生成一个n 阶单位矩阵;Insert[expr,elem,n]——在表expr 的第n 个元素前插入元素elem 产生一个新表; If[condition,t,f]——如果条件condition 为True,执行t 段,否则执行f 段;Length[expr]——表expr 中第一层元素的个数( 通常称为表的长度);LinearSolve[A,b],给出向量x,使得Ax b =成立.——求解线性方程组Ax b =的一个特解; LUBackSubstitution[数据,b]——利用LUDecomposition[矩阵Ａ]输出的结果求解方程组矩阵.x=bLUDecomposition[A]——求出矩阵A 的LU 分解;LUMatrices[lu]——返回一个形式为{l,u}的列表;MatrixForm[expr]——按矩阵形式将表expr 输出;NullSpace[A]——求齐次线性方程组AX=0的基础解系;Print[expr1,expr2,…]——顺次输出表达式expri 的值;RowReduce[A]——将矩阵A 化为简化行梯形形式;Simplify[expr]——对表达式expr 进行化简;Solve[eqns,vars]——从方程组eqns 中解出变量vars;(*information*)--注释语句,其中information 表示对程序作说明;{a,b,c}--Mathematica 中表的结构,表示元素分别为a,b,c 的表;3线性方程组解的算法及在Mathematica 软件中的实现3.1线性方程组的解考虑线性方程组：,,,()n m m n A X b X R b R R A m ⨯=∈∈= （1）其对应的齐次线性方程组为：0AX = （2）讨论线性方程组AX b =的求解问题,式中A 为m n ⨯阶系数矩阵,b 为1m ⨯阶右端列向量,X 为待求的1n ⨯阶列向量.这里首先假定m n ≤,如果m n >留待稍后讨论.当m n =且行列式0A ≠时,称AX b =为恰定方程组;当m n <时,称AX b =为不定方程组;当m n >时,称AX b =为超定方程组;在上述方程组求解的讨论中,常常要用到系数矩阵A 的秩()R A ,与增广矩阵的秩(,)R A b ,有时还要计算对应齐次方程组0AX =的基础解系.线性方程组（1）的解有以下3个结论:线性方程组的解的基本问题：当()(,)R A R A b n ==时,（1）有唯一解;当()(,)R A R A b n =<时,（1）有无穷多组解;当()(,)R A R A b <时,（1）无解.线性方程组的解的结构问题：线性方程组（1）对应的齐次方程组（2）的基础解系由()n R A -个向量组成;（2）的通解为基础解系的任意线性组合;（1）的通解为（1）的特解与（2）的通解之和.Mathematica 软件中的命令Solve 可以直接求出方程组（不限于线性）的解,一般格式为：Solve[eqns,vars],其功能是从方程组eqns 中解出变量vars.（但对于大的方程组而言,这种方法就显得有点儿麻烦,而且效率不高.）也有专门用于求解线性方程组的命令,一般格式为：LinearSolve[A,b],给出向量x,使得AX b =成立.其功能是求解线性方程组AX b =的一个特解,此命令的缺陷是当方程组的解不唯一时,只能求出一个特解而不能求通解.其中A 是方程组的系数矩阵,而b 为线性方程组的“右边项”.如果A 为可逆阵,那么LinearSolve 就给出线性方程组的唯一解.如果A 为奇异阵,那么方程组就可能没有解,也有可能存在无穷组解.如果方程组具有唯一解,Mathematica 就会给出这个解（如例3）.如果方程组没有解,Mathematica 就会返回出错信息（*Linear equation encountered which has no solution.*）,即遇到了没有解的线性方程组.（如例4）例2：求解线性方程组.⎪⎩⎪⎨⎧=++=++=++147338523532321321321x x x x x x x x xA1={{2,1,3},{3,2,5},{3,3,7}};b1={5,8,14};Det[A1]=10≠,故知方程组有唯一解LinearSolve[A1,b1] （*利用Solve 函数求解11A X b =*）={-3,-4,5} （*求得解543321=-=-=x x x *）例3:求解线性方程组⎪⎩⎪⎨⎧=++=+-=++124313427233z y x z y x z y x⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡=⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡-=12171;4313422331b a a1={{3,3,2},{2,-4,3},{1,3,4}};b1={{7},{1},{12}};运行LinearSolve[a1,b1]得：231776,,701435⎧⎫⎧⎫⎧⎫⎧⎫-⎨⎨⎬⎨⎬⎨⎬⎬⎩⎭⎩⎭⎩⎭⎩⎭例4：方程组⎪⎩⎪⎨⎧=+-=+-=++1143313462z y x z y x z y x 没有解.⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡=⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡--=11162;4333411122b a a2={{2,1,1},{1,-4,3},{3,-3,4}};b2={{6},{1},{11}};运行LinearSolve[a2,b2]得：LinearSolve::nosol: Linear equation encountered that has no solution.（遇到了没有解的线性方程组）LinearSolve[{{2,1,1},{1,−4,3},{3,−3,4}},{{6},{1},{11}}]如果方程组没有解,Mathematica 就会返回出错信息（*Linear equation encountered which has no solution.*）,即这是一个没有解的线性方程组.如果方程组Ax=b 有无穷组解,问题就有一些复杂了.对于这种情况,Mathematica 就会返回一个解,我们称其为特解.而所有解是由特解加上对应的齐次方程组Ax=0的所有解构成的.所有满足0AX =的向量X 全体称为A 的零空间,可以很容易地用NullSpace 命令确定出零空间的构造.NullSpace[A ]返回A 的零空间的基向量,即齐次线性方程组AX=0的基础解系.A 的零度,即A 的零空间中最大线性无关向量组的向量个数,可以通过运行Length[NullSpace[a]]得到. A的秩可以用n-Length[NullSpace[A]]算出,其中n表示A 的列数.例5：方程组2y74323349x zx y zx y z++=⎧⎪-+=⎨⎪-+=⎩有无穷组解.a={{2,1,1},{1,-4,3},{3,-3,4}};b={{7},{2},{9}};nullspacebasis=NullSpace[a]{{−7,5,9}}（*由于零向量非空,因此没有唯一解*）particular = LinearSolve[a, b]{{10/3}, {1/3}, {0}}（*这是一个特解*）方程组的完全解具有形式t*nullspacebasis+particular,其中t是任意参数.然而,为了表示成单个列表,我们必须展平nullspacebasis和particular. generalsolution = t*Flatten[nullspacebasis] + Flatten[particular]{10/3 - 7 t, 1/3 + 5 t, 9 t}作为检验的步骤,我们下面把一般解代回到原来的方程组中.a.x /. x -> {10/3 - 7 t, 1/3 + 5 t, 9 t} // Expand{7,2,9}关于线性方程组求解的算法研究已比较成熟,本文主要利用高斯-约当（Gauss-Jordan）消元法和矩阵LU分解法求解线性方程组.3.2高斯-约当（Gauss-Jordan）消元法求解线性方程组ax=b的高斯-约当方法是基于对增广矩阵[]a b的简化处理,这个过程通过初等变换把矩阵转化为简化行阶梯形的形式.基本初等行变换有下述三种：①交换两行②在一行中乘上一个非零常数③在一行中加上另一行的倍数.容易证明初等行变换对方程组的解没有影响.一个矩阵称为简化行阶梯形形式,是指①每一行中第一个非零元素都是1（称为首1）②在这个1的上面或者下面的元素都是0③对于两个首1的行,下面的行中第一个非零元素比前一个更靠右.在求解线性方程组时,我们首先应该使用初等行变换把增广矩阵转化为简化梯形形式.学过线性代数或者高等代数的读者都知道初等行变换是非常麻烦的,很容易出错.并且只要有一个地方出错,后面的计算就已经无用了.然而,Mathematica数学软件中的RowReduce命令可以自动把任意矩阵转化为简化行梯形形式.RowReduce[A]——把矩阵A 转化为简化行梯形形式.例6： 考虑3⨯4阶矩阵,其元素为ij a i j =-.a=Table[Abs[i-j],{i,1,3},{j,1,4}];a//MatrixForm⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡123012101210 RowReduce[a] //MatrixForm⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎣⎡23021100010001 下面我们演示行简化过程是怎么用于线性方程组的求解.为了便于比较,我们使用的范例采用前面的例3、例4与例5 中考虑的三个例子.例7：（a ）方程组124313427233=++=+-=++z y x z y x z y x （有唯一解） （b ）方程组1143313462=+-=+-=++z y x z y x z y x （无解） （c ）方程组943323472=+-=+-=++z y x z y x z y x (有无穷组解)这三个方程组的增广矩阵分别为 ⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡--=⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡--=⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡-=9433234171123;11433134161122;12431134272331a a aa1={{3,3,2,7},{2,-4,3,1},{1,3,4,12}};a2={{2,1,1,6},{1,-4,3,1},{3,-3,4,11}};a3={{2,1,1,7},{1,-4,3,2},{3,-3,4,9}}; RowReduce[a1]//MatrixForm⎪⎪⎪⎪⎪⎪⎭⎫ ⎝⎛-357610014170107023001 (*当矩阵被简化后,如果把它看作方程组,容易看出x=-23/70,y=17/14,z=76/35*) RowReduce[a2] //MatrixForm795910001-00001⎛⎫ ⎪ ⎪ ⎪⎝⎭(*简化后的矩阵第三行指的是0001x y z ++=,这当然是不可能的.这个矛盾（一行中只有最后一个元素是非零,其余都是零）说明方程组没有解.*) RowReduce[a3] //MatrixForm7109351931001-0000⎛⎫ ⎪ ⎪ ⎪⎝⎭(底行全是0没有任何不对的地方,然而这个方程组没有唯一解,如果令z=t,其中t 为独立参数,那么方程组的解具有如下形式：10715,,3939x t y t z t =-=+=) [注释] 虽然这个解与例5给出的解在形式上稍有些不同,但从它们确定相同的解集的意义上看,它们是等价的.求解线性方程组的算法步骤为：1.对增广矩阵（非齐次）或系数矩阵（齐次）进行初等行变换,将其转化为行最简形矩阵（对非齐次线性方程组判断是否有解);2.确定非自由未知量和自由未知量,将行最简形“翻译”为方程组;3.令自由未知量等于任意常数,得到通解.3.3 矩阵的LU 分解法（Doolittle 法）LU 分解法（杜利特尔Doolittle 方法） 是另外一种经常使用的求解线性方程组的方法,尤其是当具有许多方程组,而且每个方程组的系数相同时,这种方法特别有用. LU 分解法的基本思想是这样的.如果A 为方阵,那么它可以分解为A=LU,其中L 为下三角矩阵,主对角线上元素全是1,U 为上三角矩阵.这样方程组Ax=b 就转化为（LU ）x=b,其可重写为L （Ux ）=b.如果令y=Ux,那么可以从Ly=b 中解出y,一旦有了y,就可以从Ux=y 中解出x.(LU)Ly b Ax b x b Ux y =⎧=⇔=⇔⎨=⎩虽然矩阵LU 分解法把一个方程组的求解转化为两个方程组的求解,但由于现在每个方程组的系数矩阵都是三角矩阵,因此计算起来速度是很快的.这样一来,如果使用LU 分解法求解线性方程组,就包含两个步骤：分解与回代.对应的Mathematica 命令就分别为LUDecomposition 与LUBackSubstitution .LUDecomposition[A]——求出矩阵A 的LU 分解;LUBackSubstitution[数据,b]——利用LUDecomposition[矩阵Ａ]输出的结果求解方程组矩阵.x=bLUDecomposition 的输出为数据,由三部分组成：（1）矩阵L 与U 被压缩到一个矩阵中,（2）一个置换向量,以及（3）矩阵的L ∞条件数.把数据送给LUBackSubstitution 就可以求解出线1112111121421222212221212111n n n n n nn n n nn a a a r r r a a a l r r A a a a l l r ⎛⎫⎛⎫⎛⎫ ⎪ ⎪⎪ ⎪ ⎪⎪== ⎪ ⎪⎪ ⎪ ⎪⎪⎝⎭⎝⎭⎝⎭性方程组.置换向量对行进行重新安排,以保证矩阵具有最大的数值稳定性.本文不再对条件数加以讨论.LUDecomposition 与LUBackSubstitution 不能用于求解具有无穷组解的线性方程组.例8： 利用LU 分解法求解线性方程组2y 743232213x z x y z x y z ++=⎧⎪-+=⎨⎪++=⎩.为了用LU 分解法求解线性方程组,我们首先对系数矩阵进行矩阵分解.2117143;2;32213a b ⎡⎤⎡⎤⎢⎥⎢⎥=-=⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦⎣⎦运行Clear all;a={{2,1,1},{1,-4,3},{3,2,2}};b={{7},{2},{13}};data=LUDecomposition[a]得{{{1, -4, 3}, {2, 9, -5}, {3, 14/9, 7/9}}, {2, 1, 3}, 1}LUBackSubstitution[data,b]{{1},{2},{3}}在例9和例10中演示了LUDecomposition 返回数据的结构.例9 567131232728m ⎛⎫ ⎪= ⎪ ⎪⎝⎭;m={{5,6,7},{1,3,12},{3,27,28}};{lu,m,cond}=LUDecomposition[m]{{{1,3,12},{5,-9,-53},{3,-2,-114}},{2,1,3},1}在本例中,没有进行行置换,因为此时的置换向量p 就是｛1,2,3｝.LUDecomposition[m]的第一部分是以压缩形式给出的两个矩阵,因为我们已知LU 的形式为100100100x x x x x x x x x ⎡⎤⎡⎤⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦⎣⎦; 所以要想指定矩阵L 、U,只需要指定九个元素（用x 表示）.在LUDecomposition[m]的第一部分就是用单个矩阵的形式给出这九个数.lu//MatrixForm234256347⎡⎤⎢⎥⎢⎥⎢⎥⎣⎦这些数虽然是组合在一个矩阵中,但每个数所在的位置却是分别与L 、U 中的位置一致的,因此100234LU 210056341007⎡⎤⎡⎤⎢⎥⎢⎥=⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦⎣⎦例10 211143322m ⎡⎤⎢⎥=-⎢⎥⎢⎥⎣⎦; Clear all;m={{2,1,1},{1,-4,3},{3,2,2}};{lu,p,cond}=LUDecomposition[m]{{{1, -4, 3}, {2, 9, -5}, {3, 14/9, 7/9}}, {2, 1, 3}, 1}lu//MatrixForm714991432953⎛⎫ ⎪ ⎪ ⎪⎝⎭--如果仍然按前面例子同样的步骤处理,就会得到71499100143210;095;3100l u -⎡⎤⎡⎤⎢⎥⎢⎥==-⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦⎣⎦ 但这时l 与u 乘起来并不等于原来的矩阵.l.u//MatrixForm143211322-⎡⎤⎢⎥⎢⎥⎢⎥⎣⎦实际上这时置换向量为p=｛2,1,3｝,表明矩阵中第一行与第二行交换了顺序.如果这时交换l 中相应的行,再乘上U 就会得到原来的矩阵.149210100;31l ⎡⎤⎢⎥=⎢⎥⎢⎥⎣⎦L.U//MatrixForm211143322⎡⎤⎢⎥-⎢⎥⎢⎥⎣⎦重构出L 与U 矩阵的更方便的方法是利用LUMatrices 命令,这条命令包含在软件包LinearAlgebra`MatrixManipulation`中.LUMatrices[lu]——返回一个列表,形式为{l,u},由对应于lu 的矩阵在LU 分解中被置换后的上三角与下三角矩阵组成,其中lu 为运行LUDecomposition[A]所得结果中的第一个元素.例11：自动求解例9中的L 、U运行<<LinearAlgebra`MatrixManipulation`调入软件包. 运行m={{2,1,1},{1,-4,3},{3,2,2}}; {lu,p,cond}=LUDecomposition[m]得数据{{{1, -4, 3}, {2, 9, -5}, {3, 14/9, 7/9}}, {2, 1, 3}, 1} 运行LUMatrices[lu] 得转换后的L 、U 矩阵｛｛1,0,0｝,{2,1,0},{3,14/9,1}｝,{{1,-4,3},{0,9,-5},{0,0,7/9}}} 由此可见可以自动提取转换后的L 、U 矩阵： L= LUMatrices[lu][[1]]; U= LUMatrices[lu][[2]];L[[p]]//MatrixForm (*L[[p]]根据向量p 对矩阵L 的行进行置换*)14921010031⎛⎫ ⎪ ⎪ ⎪⎝⎭U//MatrixForm7914309500-⎛⎫⎪- ⎪ ⎪⎝⎭运行下列命令,从结果可见L 与U 的乘积与原矩阵一致：L[[p]].U//MatrixForm211143322⎛⎫ ⎪- ⎪ ⎪⎝⎭例12：利用LU 分解法求解方程⎪⎪⎪⎭⎫⎝⎛=⎪⎪⎪⎭⎫ ⎝⎛-400320102423xLU 分解法的计算步骤：1.首先将方程组的系数矩阵进行LU 分解,得到A LU =.2.求解LY=b.3.求解UX=Y.(*输入数据*)A = {{3, 2, -4}, {2, 0, 1}, {0, 2, 3}}; b = {0, 0, 4};(*对A 进行LU 分解：A=LU*)Print["待求方程组为", MatrixForm[A], "x=", MatrixForm[b]]; n = Length[b];{lu, p, c} = LUDecomposition[A];Print["LU 分解的紧凑格式为", MatrixForm[lu]]; Print["置换向量为", p];l = lu SparseArray[{i_, j_} /; j < i -> 1, {n, n}] + IdentityMatrix[n]; Print["分离的L 矩阵为", MatrixForm[l]];u = lu SparseArray[{i_, j_} /; j >= i -> 1, {n, n}]; Print["分离的U 矩阵为", MatrixForm[u]]; (*求解Ly=b*)y = LinearSolve[l, b[[p]]]; Print["求解Ly=b 得:y=", y]; (*求解ux=y*)x = LinearSolve[u, y] // N; Print["方程组解为：x=", x];运行以上程序得：4线性方程组求解过程的可读性计算假设AX E =,其中A 为n 阶系数矩阵,E 为单位矩阵,即123=,,,,n E e e e e （）,其中(1,2,)i e i n =为单位列向量,X 为n 个列向量构成的矩阵,即2(,,)i n X x x x =,其中,(1,2,)i x i n =为列向量.于是,可以把等式AX E =看作求解n 个线性方程组,(1,2,)i i Ax e i n ==出了所有i x 之后,也即得到了矩阵X .而由AX E =可知,矩阵X 是A 的逆矩阵,即1X A -=.这样就求出了A 的逆矩阵.于是求逆矩阵的过程被化成了解线性方程组的过程.我们利用Mathematica 软件实现对求逆矩阵的可读性计算,对我们探讨线性方程组的解法有非常重要的意义.4.1 求逆矩阵可读性计算我们编写可读性程序如下： (*求矩阵的逆矩阵*)运行InvMatr[A_] := (If[! MatrixQ[A],Print["A is not a Matrix, We don't calculate the inverse of A."], Print["Matrix A =", MatrixForm[A], ", "]; If[Det[A] == 0,Print["Because of Det[A]=0, There isn't inverse Matrix of A."], n = Length[A];ident = IdentityMatrix[n]; AT = Transpose[A]; Temp = AT;Do[Temp = Append[Temp, ident[[i]]], {i, 1, n}]; AE = Transpose[Temp];Print["(A|E)=", MatrixForm[AE]]; Do[k = Infinity; r = j;Do[ If[Abs[AE[[i, j]]] < k && AE[[i, j]] != 0, k = AE[[i, j]]; r = i], {i, j, n}];If[r != j, TempMatrix = AE[[r]]; AE[[r]] = AE[[j]]; AE[[j]] = TempMatrix];Do[AE[[i]] = AE[[i]] - AE[[j]]*AE[[i, j]]/AE[[j, j]], {i, 1, j - 1}];Do[AE[[i]] = AE[[i]] - AE[[j]]*AE[[i, j]]/AE[[j, j]], {i, j + 1, n}]; AE[[j]] = AE[[j]]/AE[[j, j]];Print["-->", MatrixForm[AE]]; , {j, 1, n} ];B = Transpose[AE]; Temp = {B[[n + 1]]};Do[Temp = Append[Temp, B[[i]]], {i, n + 2, 2 n}];Print["The Inverse Matrix A of is ", Transpose[Temp] // MatrixForm, "."] ] ] )为说明本程序可用于任何对矩阵的逆的各种情形的求解处理,举例计算如下.例13：已知a1={2,3,1,4};a2={3,0,2,5};a3={4,3,8};a4={2,7,4,15};A={a1,a2,a3,a4},求A 的逆阵. 运行Clear; a1={2,3,1,4}; a2={3,0,2,5}; a3={4,3,8}; a4={2,7,4,15}; A={a1,a2,a3,a4}; InvMatr[A]得：A is not a Matrix, We don't calculate the inverse of A. 如果我们输入的A 不是矩阵,我们不可求A 的逆阵.例14：求矩阵A=⎪⎪⎪⎪⎪⎭⎫⎝⎛64401800154012202的逆矩阵. 运行Clear;a1={2,0,2,2}; a2={1,0,4,5}; a3={1,0,0,8}; a4={1,0,4,64}; A={a1,a2,a3,a4}; InvMatr[A]得：Matrix A =⎪⎪⎪⎪⎪⎭⎫ ⎝⎛64401800154012202 Because of Det[A]=0, There isn't inverse Matrix of A.此时,A 的行列式为0,A 没有可逆阵.例15：求矩阵A=⎪⎪⎪⎭⎫ ⎝⎛--242022233的逆矩阵.Clear;a1={3,3,2};a2={2,-2,0}; a3={-2,4,2}; A={a1,a2,a3}; InvMatr[A]运行结果为：例16：求矩阵A= ⎪⎪⎪⎪⎪⎪⎭⎫⎝⎛3000023000023000023000023的逆矩阵Clear;a1={3,2,0,0,0};a2={0,3,2,0,0};a3={0,0,3,2,0};a4={0,0,0,3,2}; a5={0,0,0,0,3};A={a1,a2,a3,a4,a5}; InvMatr[A]通过上面的讨论我们可以看出该程序对任何一个矩阵型数据均能做出相应的处理,对输入的矩阵为可逆时,程序能计算出逆矩阵.进而我们通过Mathemtica 软件对解方程组也可实现可读性计算.4.2 线性方程组的可读性计算虽然Mathematica 软件中有解方程组的函数,但是,用此函数求解方程组存在不足之处：不能看到求解方程组的过程,也不能看到求解方程组的其他信息;输入方程组信息的形式繁琐且输出的内容不符合代数教材中通常的习惯,不能为其他程序提供更多的信息.而通过Mathematica 软件到底能不能对解线性方程组算法的实现将其过程展现出来呢？我们的答案是肯定的.只是在求解过程中,要运用到更多的输出形式函数．其程序相对较复杂.程序中较多的命令都是由简单命令复合而成,现对复合命名的构造过程进行描述.为了使实例17中的增广矩阵2-21-1121-42-2334-103-57812-11-2-1⎛⎫ ⎪ ⎪ ⎪⎪⎝⎭具有分块形式22111|214223|3410357|812112|1⎛⎫⎪⎪⎪ ⎪⎝⎭---------, 就需要在最后一列前插入符号“|”进行分隔,但Mathematica 中的矩阵是由行表构造而成,因此先用Transpose 命令对矩阵AUB 进行转置,并事先生成一个全由“|”为元素的一维表B0,一次性地把B0 插入到AUB 的转置矩阵的最后一行之前,然后再转置,就实现了我们希望的输出效果．完整的命令如下：B0=Table["|",{i,1,m}];MatrixForm[Transpose[Insert[Transpose[AUB],B0,n+1]]]Mathematica对表进行操作时是对表的每个元素进行同一操作,因此在对AUB 中找出的主元AUB(r,j)进行归一化处理时,可用下列命令实现对r行元素各除以AUB(r,j)：AUB[[r]]=AUB[[r]]/AUB[[r,j]];为了把j 列除主元AUB(r,j)外的各元素全变成零,可对1 至r-1 行及r+1 至m 行分别进行第三类初等行变换．完整命令如下：Do[AUB[[i]]=AUB[[i]]-AUB[[r]]*AUB[[i,j]] ,{i,1,r-1}];Do[AUB[[i]]=AUB[[i]]-AUB[[r]]*AUB[[i,j]] ,{i,r+1,m}];处理解唯一时的输出时要解决的问题：AUB 的第n+1 列就是所求的解,但当m>n 时由于方程个数比未知量个数多,直接输出AUB(n+1)将出现多余的0,因此先把这些0 的位置存入表dele 中,再用Delete 命令把AUB(n+1)中多余的0 删除;为了把行向量以列的方式输出,由于Transpose 对一维表无效,因此把删除多余元素后的表放在{}中以增加其维数,最后形成的完整命令如下：dele=Table[{i},{i,n+1,m}];MatrixForm[Transpose[{Delete[Transpose[Simplify[AUB]][[n+1]], dele]}]]]下面是经过完整的求解线性方程组的可读性计算的完整程序：conslist = {C1, C2, C3, C4, C5, C6, C7, C8, C9, C10, C11, C12, C13,C14, C15, C16, C17, C18, C19, C20, C21, C22, C23, C24, C25};varlist = {x1, x2, x3, x4, x5, x6, x7, x8, x9, x10, x11, x12, x13,x14, x15, x16, x17, x18, x19, x20, x21, x22, x23, x24, x25}; MatrixPrint[A_] := (B0 = Table["|", {i, 1, m}];MatrixForm[Transpose[Insert[Transpose[AUB], B0, n + 1]]])PrintSolve[] := (ME = IdentityMatrix[S];basV = Complement[Table[i, {i, 1, n}], fuq];L = Length[basV];(*L 记录基变量的个数*)AUB = Transpose[AUB];eta = Table["", {S}];i = 0; eta[[i]] = AUB[[n + 1]];eta[[i]] = Table[eta[[i]][[j]], {j, 1, L}];Do[eta[[i]] = Insert[eta[[i]], 0, fuq[[k]]], {k, 1, S}];Do[eta[[i]] = -AUB[[fuq[[i]]]];eta[[i]] = Table[eta[[i]][[j]], {j, 1, L}];Do[eta[[i]] = Insert[eta[[i]], ME[[i, k]], fuq[[k]]], {k, 1,S}], {i, 1, S}];Ck = Table[conslist[[i]], {i, 1, n}];solv = MatrixForm[Transpose[{eta[[0]]}]];Do[solv = solv + MatrixForm[Transpose[{eta[[i]]}]]*Ck[[i]], {i, 1,n - r}];Print["X=", solv];)SolutionEquation[A_, B_] := (n = Length[A[[1]]]; m = Length[B];X = Table[varlist[[i]], {i, 1, n}]; AX = Dot[A, X];Print["求解下列线性方程组:"];Do[Print[AX[[i]], "=", B[[i]]], {i, 1, m}];Print["解:先对增广矩阵进行初等行变换:"];AUB = Transpose[Append[Transpose[A], B]];Print[MatrixPrint[Simplify[AUB]]];r = 1; fuq = {};(*fuq 存放自由未知量的下标*)Do[k = Infinity; flag = 0;(*查找列主元素的位置*)Do[If[AUB[[i, j]] != 0, flag = 1], {i, r, m}];If[flag == 0, fuq = Append[fuq, j]; Continue,Do[If[Abs[AUB[[i, j]]] < k && AUB[[i, j]] != 0, k = AUB[[i, j]]; t = i], {i, r, m}];If[t != r, TempMatrix = AUB[[t]];AUB[[t]] = AUB[[r]]; AUB[[r]] = TempMatrix];AUB[[r]] = AUB[[r]]/AUB[[r, j]];Do[AUB[[i]] = AUB[[i]] - AUB[[r]]*AUB[[i, j]], {i, 1, r - 1}];Do[AUB[[i]] = AUB[[i]] - AUB[[r]]*AUB[[i, j]], {i, r + 1, m}];Print["-->", MatrixPrint[Simplify[AUB]]];r = r + 1];, {j, 1, n}];r = r - 1;S = n - r;(*基础解系中解的个数*)Print[" 增广矩阵的秩r=", r, ", 未知量个数 n=", n, ";"];If[r == n, Print["结论：因为增广矩阵的秩等于未知量个数,方程组有唯一解：X=", dele = Table[{i}, {i, n + 1, m}];MatrixForm[Transpose[{Delete[Transpose[Simplify[AUB]][[n + 1]], dele]}]]], If[r == m, Print["结论：因为增广矩阵的秩等于方程个数且小于未知量个数,方程组有无穷多组解："];PrintSolve[],If[AUB[[r + 1, n + 1]] != 0, Print["因为增广矩阵的秩小于未知量个数,但 dr+1!= 0,所给方程组无解"], Print["因为增广矩阵的秩小于未知量个数, 且dr+1=0,方程组有无穷多组解："];PrintSolve[];];];];)说明：本程序可用于变量在25个内的任何线性方程组的求解,如果变量超过25个,只需要增加程序中conslist 与varlist 分量个数即可.例17：求解方程组⎪⎪⎩⎪⎪⎨⎧=-+-+=+-+-=+-+-=+-+-2342310634852232533354321543215432154321x x x x x x x x x x x x x x x x x x x x 在Mathematica 输入并运行：A={{3,-3,1,-1,1},{1,-5,2,-3,2},{5,-8,4,-3,6},{1,3,-2,4,-3}};b={3,2,10,2};SolutionEquation[A,b]即得到如下求解结果,其中所有内容全部由程序自动输出得到. 求解下列线性方程组: 3x1-3x2+x3-x4+x5=3 x1-5x2+2x3-3x4+2x5=2 5x1-8x2+4x3-3x4+6x5=10 x1+3x2−2x3+4x4−3x5=2例18：求解方程组⎪⎪⎩⎪⎪⎨⎧-=+-++=-++--=+-++=+-++1423412231245354321543215432154321x x x x x x x x x x x x x x x x x x x x 在Mathematica 输入并运行：A={{1,3,5,-4,2},{1,3,2,-2,1},{1,-4,1,1,-1},{1,2,1,-4,1}};b={1,-1,3,-1}; SolutionEquation[A,b]即得到如下求解结果,其中所有内容全部由程序自动输出得到.例19：求解线性方程组⎪⎪⎩⎪⎪⎨⎧-=+-=++=--=-+332523233231321321321x x x x x x x x x x x a1={1,2,-1};a2={1,-3,-1};a3={3,1,2};a4={-2,0,3};A={a1,a2,a3,a4};B={3,2,5,-3}; SolutionEquation[A,B]参考文献[1] 北京大学数学系几何与代数教研室前代数小组编.高等代数[M].北京：高等教育出版社,2003[2] 王绍恒,王艺静.利用Mathematica软件实现解线性方程组的可读性计算[J].重庆三峡学院学报,2009(3),25（3）:143-146[3] 姜友谊,应宏,王绍恒.化实对称矩阵为对角矩阵的计算机算法[J]．西南民族学院学报,2002,28(4)：428-432[4] 张禾瑞,郝鈵新,高等代数(第五版)[M].北京：高等教育出版社,2007[5] 刘水强,王绍恒．利用初等行变换解线性矩阵方程[J]．重庆三峡学院学报,2001,17(5)：87-89[6] 徐士良.计算机常用算法[M].北京:清华大学出版社,2000.8．[7] 徐安农,科学计算引论：基于Mathematica的数值分析[M].北京:机械工业出版社,2010.8[8] 冯天祥,数值计算方法理论与实践研究[M].成都：西南交通大学出版社,2005[9] 阳明盛,林建华.Mathematica基础及数学软件[M],大连：大连理工大学出版社,2003[10][美] Rchard J. Gaylord,Samuel N.Kamin ,Paul R. Wellin著,耿勇译.数学软件Mathematica 入门[M].高等教育出版社,2001[11][美]D.尤金.Mathematica使用指南[M].邓建松,彭冉冉,译.北京：科学出版社,2002.。