maple矩阵求解

。

在载入这个工具包时产生了两个警告，不用为此担心，这是因为Maple系统中已有的两个函数——norm和trace被linalg中的同名函数所覆盖了。原来系统中这两个函数分别是用于计算多项式的范数（norm）和程序调试跟踪的，一般情况下，在解决线性代数问题时用不到。如果确实要用系统中的库函数，可以再用readlib将其重新载入进来。

linalg中的主要函数的详细介绍我们将留到本章第三节中完成，在这一节中，我们主要着眼于矩阵的基本运算，比如加法、乘法和乘幂等。我们先用上一章中介绍过的方法来定义一些矩阵：

矩阵的代数运算最直接而且直观的方法莫过于使用函数evalm( )了，只要把矩阵的代数计算表达式作为它的参数，就可以得到结果的矩阵。

矩阵和标量也可以直接作和差。但要注意的是，和Matlab以及Fortran90中的相应的运算定义不同，在Maple中，矩阵和标量相加（相减），被定义成矩阵的对角元和标量分别相加（相减），也就是矩阵加上（减去）和它相同形状的数量阵。请看下面的例子：

还需要说明的是，运算符“*”仅仅用来表示数乘，对于矩阵乘法，绝对不能使用。这是因为，运算符“*”在Maple中被定义成满足交换律的乘法，而我们知道，矩阵乘法是不满足交换律的。我们可以通过简单的例子来说明这个事实。

可以看到，Maple在化简表达式时，丝毫没有顾及它们的类型，只把它们作为一般的符号变量来计算，甚至不考虑这些矩阵是否可以相乘或相加。如果我们用evalm来计算，也就得不到正确的结果了。在Maple V Release 4及以后的版本中，程序设计者考虑到用户初次使用时可能对此不了解，特别给除了错误信息。

但如果你使用的是早期版本，就不但没有错误信息，还会给出错误的结果，使用时一定要特别注意。上面的错误信息告诉我们，矩阵/向量运算必须使用&*作为乘法运算符。让我们再来试一试。

这次，我们得到了期望的结果。

矩阵也可以进行乘幂运算，矩阵的乘幂运算符和单个表达式的相同——“^”，当然，也需要用evalm来计算得到结果矩阵。

Maple不仅可以计算矩阵的自然数次幂，同样也可以计算矩阵的负整数次幂。

在计算矩阵的0次幂时，却不怎么令人满意，它的结果令人啼笑皆非：

这是由于Maple的表达式自动化简机制造成的，B^0在还没有进入evalm函数之前，就已经被系统化简成了1！

矩阵乘法运算符“&*”有着和数乘“*”、除法“/”相同的优先级，在输入表达式是需要注意，必要时要使用括弧。

还好，Maple没有再次给出1，Maple用“&*( )”来表示单位阵。

4.1矩阵求值

我们早在上一章中就已经注意到，数组和单个的表达式不同，具有数组类型的变量不会自动求值，而需要用eval( )等函数的显式调用来强制地求值。作为数组的特例，矩阵和向量也是一样。更令人吃惊的是，当我们用whattype来查看数组变量的表观数据类型时，返回的居然是符号类型（symbol），和被赋值之前一样！

然而，通过进一步的试验你还会发现，即使用eval( )，还是没办法让矩阵的每一个元素都求值。只有通过我们在上一章中用过的函数map，才可以完成最终的求值。

同样的情况也发生在以查找表（table）和函数或子程序（procedure）为类型的对象上：

在处理这些类型的对象时，Maple采取的是按名求值的方法，而不是像对于单个表达式时采取的完全求值的方法。所谓“按名求值”，是指系统在计算是只将变量的名称代入进行计算，而不读取该变量名所指向的数据。如果需要得到变量名所指的值，则必须强制使用eval甚至是map和eval并用。通过一个简单的例子，也许可以更容易地理解这种求值方式。

我们通过对T的不同层次的求值，获得了不同的对象，很明显，在Maple内部，数据结构是用类似于下图的方式组织的：

如果我们用系统的自动求值，则只能到Q的一层为止，就好象从Q指向数组的指针上有一道难以逾越的屏障：

不管我们怎么赋值，T、P、Q这些变量总是直接或间接地指向同一个数据对象，因为在赋值时，用的也是“按名求值”。用面向对象的观点来看，它们都这个矩阵变量的引用，而不是它的一个拷贝。当我们修改了其中任意一个变量的元素时，其他变量的元素也会跟着改变。在这上面，Maple和有一点向一些面向对象的编程语言，比如Java，但和C++是不同的。

当然，有时候我们要的是一个矩阵的拷贝，我们将改变这个拷贝的原素值，而同时保持原来的矩阵不变；这时，我们需要用Maple的拷贝函数copy( )，它可以用来拷贝这些“按名求值”的对象。

4.2矩阵和线性方程组

linalg是Maple中的线性代数工具包，有着很全面的线性代数计算工具。那它里面到底有些什么工具呢？你可以试试在载入它时，不用冒号而是用分号结尾，那样，它的所有函数就会被一一打印出来。这可是一个长长的有序表，因为linalg中的函数多达114个。有一定数学基础的读者，也许看看它的函数名称，就可以大概地知道linalg中间都有些什么工具了。自然，矩阵的基本运算是少不了的，还有矩阵的行、列操作，高斯消元，矩阵的行列式、秩、逆、迹，特征多项式、特征矩阵，矩阵的特征值和特征向量……，举不胜举。从这一节起，我们分类将一些常用的函数作一下简单介绍；其他的一些函数，请读者在使用时参考在线帮助。

4.2.1矩阵基本运算

在上一节中，我们已经介绍过矩阵的基本运算，那时，我们用的是函数evalm( )配合运算符。作为函数式的编程语言，Maple同时也提供了一套函数来完成这些矩阵基本运算，以使不同的场合下可以选用合适的方式进行运算。我们不准备详细的介绍它们的使用，仅将它们列成下表，供读者在使用时参考。

表4.1 基本矩阵运算

有时在处理大矩阵时，常常可以把它们视为是由一些小矩阵组成的，这样，大矩阵的运算就转化成小矩阵间的运算，往往可以带来极大的方便。Maple是一个符号演算软件，在这方面应该具有一定的优势。令人遗憾的是，直到Maple V Release 5为止，Maple中还无法表示一个抽象的矩阵（只有矩阵名称，而不知道矩阵的大小）。因此，严格说来，在Maple中还无法作符号表示的抽象分块矩阵运算。

但是，在Maple中对于分块矩阵在一定程度上还是支持的。Maple对于分块矩阵的支持，主要表现在输入和表示某些大而且部分重复的矩阵上。或者，我们可以用已有的矩阵拼合成一个新的矩阵。

对于除了对角元附近，其他大部分元素为0的矩阵，可以用diag函数更为便捷地生成，它可以将一些方阵（并不需要具有相同大小）或者标量排列在矩阵的对角线上，生成块状对角阵。

4.2.3初等变换

矩阵的初等变换是各种消去法的基础，是解线性方程组的基础。对于符号演算，有时候我们不仅要求最后的求解结果，而且要求中间的求解步骤。这时候我们可以调用Maple线性代数工具包中的初等变换函数。对于初学线性代数的理工科低年级学生，还可以用这些函数透彻、直观地理解各种消元法的实质。

表4.2 矩阵的初等变换

表4.3 辅助变换函数

例4.1 利用初等变换判断矩阵A是否可逆，如果A可逆，求A– 1

首先，我们将A和单位阵合成一个大矩阵AI：

接下来，通过初等变换，把AI的左边部分变成为单位阵。先用第一行消去其他两行的第一个元素：

再用第二行消去其他两行的第二个元素：

最后用第三行消去其他两行的第三个元素：

于是，这个过程说明了A是可逆的，同时也求出了A的逆矩阵——上面最后一个矩阵的右边部分。我们可以加以验证：

4.2.4线性方程组的解

线性方程组也是一般代数方程组的一个特例，所以，我们当然也可以用曾经多次接触过的代数方程求解函数solve( )来对它们进行求解。由于线性方程组的特殊性，在Maple的线性代数工具包中也有很多专门用于处理线性方程组的函数。

在线性代数中，我们通常是把线性方程组转化为与其等价的矩阵问题来解决的。考虑到这一点，Maple软件的编制者在linalg工具包中也提供了在线性代数方程组与矩阵之间相互转换的函数：

同样，用函数geneqns( )可以把矩阵转换成指定未知数的线性方程的集合。

解线性方程组的最简单而且最基本方法是高斯消去法（Gaussian elimination），在linalg工具包中，有与此方法对应的函数gausselim( )。作为例子，我们把上面得到的矩阵利用高斯消去法转化成相抵的上三角阵：

该函数还可以有两个可选参数，可以用来返回作为高斯消元的副产品的矩阵的秩，以及矩阵的行列式（仅当矩阵为方阵时有效）。

我们发现，在高斯消元的结果矩阵中出现了许多分数，这使表达式看起来不太美观。而且，通过进一步的试验，你会发现随着矩阵的增大，最后的分数的分母将会愈来愈大；不仅如此，如果系数矩阵中包含有符号变量，那么结果矩阵中将出现分式。使用一般高斯消去法时，这些都将是无法避免的。实际上，我们可以在作初等变换时，避免作除法，而是将被消去的行乘以相应的因子。这种方法称为无分式高斯消去法（fraction-free Gaussian elimination），Maple中也有这样的函数——ffgausselim( )。

在把系数矩阵转化成相抵的上三角阵后，我们就可以用回代（back substitution）的方法求出各个未知数的值。

和回代相对，我们还有前代（forward substitution）的概念。我们用前代函数forwardsub( )求解下三角矩阵。有了这一对函数，我们就可以利用LU分解来解线性方程组了。在linalg工具包中，LU分解函数是LUdecomp( )。上面，我们使用的是高斯消去法，得到的是上三角阵，需要通过回代得到最终结果。我们也可以利用高斯-约当消去法（Gauss-Jordan elimination），把系数矩阵变换成单位阵，就可以直接得到结果了。

上面介绍的实际上都是线性方程组求解的中间步骤，我们当然也可以一步到位，不管它到底用的什么方法，而只要求得到最终的解。这时，我们可以使用线性方程组求解函数linsolve( )。更重要的是，这个函数不仅可以用来求解具有唯一解的线性方程组，而且可以在方程组有多解时给出通解，如下例所示。

例4.4 解方程组

可以看到，Maple用辅助变量_t1，_t2给出了方程组的通解。

4.3.5正定矩阵

在线性代数中，我们把相应的二次型为正定二次型的矩阵称作是正定矩阵，相应的，我们有许多判定矩阵正定性的方法。但是在Maple中，判定一个矩阵是否正（负）定十分简单，只需要调用函数definite( )即可。definite可以用来判定数值矩阵的正（负）定性，也可以求出符号矩阵的正（负）定条件。

definite的第一个参数需要判定的矩阵，第二个参数是'positive_def'，'positive_semidef'，'negative_def'，或者'negative_semidef'之一，分别表示正定、半正定、负定和半负定。在判定数值矩阵时，它返回布尔值true / false；判定符号矩阵时，它返回一个布尔表达式，表示正/负定的条件，例如：

4.3.6特殊矩阵

利用linalg工具包中的一些函数，可以直接生成一些特殊的矩阵，而不需要提供指标函数，使用非常方便。我们将它们列成下表，供读者参考。

表4.4 linalg工具包中的特殊矩阵

4.4.1基本子空间

我们从线性代数中知道与矩阵A相关的R n的四个基本子空间——行空间、列空间、化零空间和左零空间。化零空间，实际上就是齐次线性方程组Ax = 0的解空间。我们可以利用上一节中介绍的函数linsolve( )来求得它的通解，在将通解重的辅助变量依次替换成1，就可以获得它的解空间的一组基。对于左零空间，即方程组ATx = 0的解空间，亦可以用同样的方法得到。

另两个基本子空间，行空间和列空间，在linalg工具包中有专门的函数rowspace( )和colspace( )，可以获得它们的基以及维数。

函数rowspace和colspace的第二个参数是可选参数，用来返回行空间或列空间的维数，它必须是一个变量名——可以是一个未赋值的新变量，也可以是已有值的变量，但要加上延迟求值符“'”（参见第1章）。

4.4.2正交基和Schmidt正交化

在欧氏空间中，我们可以定义两个向量的内积（inner product），在此基础上，我们还可以定义两个向量的夹角。如向量、的夹角可以定义为。在Maple的linalg工具包中，有相应的函数可以计算向量的内积和夹角。下面将通过例子来说明它们的用法。

在线性代数中，我们定义内积为零向量相互之间正交；并且，我们利用Schmidt正交化方法（Gram-Schmidt orthogonalization process），由欧氏空间中一组普通基得到两两正交的基。Maple中相应的函数是GramSchmidt( )，它的输入参数是由一组向量组成的集合（或有序表），它将给出Schmidt正交化后的向量集合（或有序表）。输入的向量必须是线性无关的，否则，结果向量间也将线性相关。函数并不对向量进行单位化，如果需要得到一组正交标准化基，还需要用map方法对这些向量使用单位化函数normalize( )。

例4.5 用Schmidt正交化方法，由下列向量组构造出一组标准正交向量组

(1, 1, –1, –2)T, (5, 8, –2, –3)T, (3, 9, 3, 8)T

向量的Schmidt正交化过程实际上就给出了矩阵的QR分解：A = QR，其中Q是正交矩阵，R是上三角阵，这个分解是存在而且唯一的。在Maple中，有QR分解函数QRdecomp( )，利用它可以完成数值或符号矩阵的QR分解。

QRdecomp的返回值是上三角阵R，正交阵Q和矩阵的秩可以通过函数的可选参数返回，但这里这两个参数必须指明它们的形参名——Q和rank，这和Fortran语言的函数调用有一定的相似之处。有关这方面的知识，我们将在程序设计章节中专门介绍。

那么这里的Q到底是不是一个正交阵呢？可以使用函数orthog( )来判定它的正交性：

4.4.3线性方程组的最小二乘解

在实际问题中，由于误差或者其他各方面的原因，很容易出现无解的方程组。但问题还是要解决的，我们必须给出一个“最优”的近似解。在线性代数中，我们通常采用最小二乘解（least-squares solution）。linalg中对应的函数是leastsqrs( )，它有两种调用格式，可以用来求解矩阵形式的方程组，和linsolve类似；也可以直接求借用解析表达式给出的方程组，这时调用的格式和solve函数一样。

例4.6 求不相容方程组

在最小二乘意义下的最优解。

4.5特征值、特征向量和相似标准型

4.5.1矩阵的相似

我们知道，对于方阵A，B，如果存在一个可逆的方阵P满足B = P–1AP，则称A，B相似。在Maple中，我们可以利用linalg中的函数issimilar( )判断两个矩阵是否相似。作为例子，我们来验证一个简单的定理——在A可逆时，AB和BA相似。

这里，Maple作了一定的假设，由于A的行列式是一个符号表达式，Maple在判定时假设它为非零常数。不过更多的情况下，这个函数是用在数值矩阵的相似判定上，同时，它还可以求出可逆矩阵P。

4.5.2特征值和特征向量

特征值问题可以说是最常见的线性代数问题了，许多数学问题（比如常微分方程组的求解、二阶张量的主分量）和工程问题（比如振动模态辨识、系统辨识）都会涉及到特征值问题。求方阵的特征值，实际上是求与之相似的对角阵。在线性代数中，求解矩阵A的特征值，我们用求解方程det( I – A ) = 0方法。矩阵I – A称为A的特征矩阵（characteristic matrix），而行列式det( I – A )展开得到的多项式称为特征多项式（characteristic polynomial），在linalg工具包中，分别有函数charmat( )和charpoly( )与它们相对应。

得到了特征多项式以后，我们就可以利用solve求出矩阵的特征值了。不过，在linalg中也有可以直接求出矩阵的特征值和特征向量的函数eigenvalues( )和eigenvectors( )。

我们看到，特征向量函数eigenvectors在给出特征向量的同时，还给出了对应的特征值和特征值的重数。一般情况下，如果同时需要得到矩阵的特征值和特征向量，可以直接调用函数eigenvectors。这两个函数不仅可以用来求解数值矩阵的特征值问题，也可以处理符号矩阵；在处理符号矩阵时，它们通常都会以根式的形式给出结果，如果加上可选参数'implicit'，结果就以RootOf的形式给出。

函数eigenvals( )不仅可以求解普通特征值问题，还可以求解广义特征值问题，也就是求解方程det( C – A ) = 0。广义特征值问题的调用格式是eigenvals( A, C )。

由特征值理论，任何一个实对称阵，都可以用求特征值的方法将它对角化；但对于非对称矩阵，或者是复矩阵，就没有这样的保证了。不过对于任意复矩阵A，还是可以把它化成相似的约当标准型J，使得P– 1AP = J。linalg中的约当标准型函数是jordan( )。它在求得矩阵的约当标准型的同时，还可以给出转换矩阵P。

约当标准型的定义是：由若干约当块组成的块状对角阵，即

其中，约当块J(, k)的定义如下：

在Maple中，可以调用函数JordanBlock( , k )可以构造约当块J( , k )，利用它和diag函数连用，就可以生成具有约当标准型的矩阵。