maple矩阵求解.docx

maple矩阵求解.docx

《maple矩阵求解.docx》由会员分享，可在线阅读，更多相关《maple矩阵求解.docx（20页珍藏版）》请在冰豆网上搜索。

maple矩阵求解

。

在载入这个工具包时产生了两个警告，不用为此担心，这是因为Maple系统中已有的两个函数——norm和trace被linalg中的同名函数所覆盖了。

原来系统中这两个函数分别是用于计算多项式的范数（norm）和程序调试跟踪的，一般情况下，在解决线性代数问题时用不到。

如果确实要用系统中的库函数，可以再用readlib将其重新载入进来。

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

我们先用上一章中介绍过的方法来定义一些矩阵：

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

矩阵和标量也可以直接作和差。

但要注意的是，和Matlab以及Fortran90中的相应的运算定义不同，在Maple中，矩阵和标量相加（相减），被定义成矩阵的对角元和标量分别相加（相减），也就是矩阵加上（减去）和它相同形状的数量阵。

请看下面的例子：

还需要说明的是，运算符“*”仅仅用来表示数乘，对于矩阵乘法，绝对不能使用。

这是因为，运算符“*”在Maple中被定义成满足交换律的乘法，而我们知道，矩阵乘法是不满足交换律的。

我们可以通过简单的例子来说明这个事实。

可以看到，Maple在化简表达式时，丝毫没有顾及它们的类型，只把它们作为一般的符号变量来计算，甚至不考虑这些矩阵是否可以相乘或相加。

如果我们用evalm来计算，也就得不到正确的结果了。

在MapleVRelease4及以后的版本中，程序设计者考虑到用户初次使用时可能对此不了解，特别给除了错误信息。

但如果你使用的是早期版本，就不但没有错误信息，还会给出错误的结果，使用时一定要特别注意。

上面的错误信息告诉我们，矩阵/向量运算必须使用&*作为乘法运算符。

让我们再来试一试。

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

矩阵也可以进行乘幂运算，矩阵的乘幂运算符和单个表达式的相同——“^”，当然，也需要用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基本矩阵运算

运算

函数调用

等效的运算符运算

加法

matadd（A,B）

evalm（A+B）

数乘

scalarmul（A,expr）

evalm（A*expr）

乘法

multiply（A,B,...）

evalm（A&*B&*...）

逆矩阵

inverse（A）

evalm（1/A）或evalm（A^（-1））

转置

transpose（A）

无

行列式

det（A）

无

秩

rank（A）

无

迹

trace（A）

无

4.2.2分块矩阵

有时在处理大矩阵时，常常可以把它们视为是由一些小矩阵组成的，这样，大矩阵的运算就转化成小矩阵间的运算，往往可以带来极大的方便。

Maple是一个符号演算软件，在这方面应该具有一定的优势。

令人遗憾的是，直到MapleVRelease5为止，Maple中还无法表示一个抽象的矩阵（只有矩阵名称，而不知道矩阵的大小）。

因此，严格说来，在Maple中还无法作符号表示的抽象分块矩阵运算。

但是，在Maple中对于分块矩阵在一定程度上还是支持的。

Maple对于分块矩阵的支持，主要表现在输入和表示某些大而且部分重复的矩阵上。

或者，我们可以用已有的矩阵拼合成一个新的矩阵。

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

4.2.3初等变换

矩阵的初等变换是各种消去法的基础，是解线性方程组的基础。

对于符号演算，有时候我们不仅要求最后的求解结果，而且要求中间的求解步骤。

这时候我们可以调用Maple线性代数工具包中的初等变换函数。

对于初学线性代数的理工科低年级学生，还可以用这些函数透彻、直观地理解各种消元法的实质。

表4.2矩阵的初等变换

函数调用

对应的初等变换

行变换

mulrow（A,r,expr）

用标量expr乘以矩阵A的第r行

addrow（A,r1,r2,m）

将矩阵A的第r1行的m倍加到第r2行上

swaprow（A,r1,r2）

互换矩阵A的第r1行和r2行

列变换

mulcol（A,c,expr）

用标量expr乘以矩阵A的第c列

addcol（A,c1,c2,m）

将矩阵A的第c1列的m倍加到第c2列上

swapcol（A,c1,c2）

互换矩阵A的第c1列和c2列

除了这些初等变换外，Maple在linalg工具包中还提供了一些矩阵的形状变换函数，可以在利用初等变换解决问题时起到辅助作用，我们也将它们一并列出：

表4.3辅助变换函数

函数调用

所作的变换或操作

concat（A,B,...）

将矩阵A，B，…在水平方向上合并成一个矩阵

stackmatrix（A,B,...）

将矩阵A，B，…在竖直方向上合并成一个矩阵

row（A,i）

取矩阵A的第i行

col（A,i）

取矩阵A的第i列

row（A,i..k）

取矩阵A由第i行到第k行

col（A,i..k）

取矩阵A由第i列到第k列

delrows（A,r..s）

删除矩阵A中第r行到第s行后剩下的子矩阵

delcols（A,r..s）

删除矩阵A中第r列到第s列后剩下的子矩阵

extend（A,m,n,x）

扩展矩阵m行n列，并用x填充

submatrix（A,m..n,r..s）

取矩阵A的m-n行，r-s列组成的子矩阵

我们通过一个简单的例题来演示这些函数的用法。

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

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

接下来，通过初等变换，把AI的左边部分变成为单位阵。

先用第一行消去其他两行的第一个元素：

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

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

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

我们可以加以验证：

4.2.4线性方程组的解

线性方程组也是一般代数方程组的一个特例，所以，我们当然也可以用曾经多次接触过的代数方程求解函数solve（）来对它们进行求解。

由于线性方程组的特殊性，在Maple的线性代数工具包中也有很多专门用于处理线性方程组的函数。

在线性代数中，我们通常是把线性方程组转化为与其等价的矩阵问题来解决的。

考虑到这一点，Maple软件的编制者在linalg工具包中也提供了在线性代数方程组与矩阵之间相互转换的函数：

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

解线性方程组的最简单而且最基本方法是高斯消去法（Gaussianelimination），在linalg工具包中，有与此方法对应的函数gausselim（）。

作为例子，我们把上面得到的矩阵利用高斯消去法转化成相抵的上三角阵：

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

我们发现，在高斯消元的结果矩阵中出现了许多分数，这使表达式看起来不太美观。

而且，通过进一步的试验，你会发现随着矩阵的增大，最后的分数的分母将会愈来愈大；不仅如此，如果系数矩阵中包含有符号变量，那么结果矩阵中将出现分式。

使用一般高斯消去法时，这些都将是无法避免的。

实际上，我们可以在作初等变换时，避免作除法，而是将被消去的行乘以相应的因子。

这种方法称为无分式高斯消去法（fraction-freeGaussianelimination），Maple中也有这样的函数——ffgausselim（）。

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

和回代相对，我们还有前代（forwardsubstitution）的概念。

我们用前代函数forwardsub（）求解下三角矩阵。

有了这一对函数，我们就可以利用LU分解来解线性方程组了。

在linalg工具包中，LU分解函数是LUdecomp（）。

上面，我们使用的是高斯消去法，得到的是上三角阵，需要通过回代得到最终结果。

我们也可以利用高斯-约当消去法（Gauss-Jordanelimination），把系数矩阵变换成单位阵，就可以直接得到结果了。

上面介绍的实际上都是线性方程组求解的中间步骤，我们当然也可以一步到位，不管它到底用的什么方法，而只要求得到最终的解。

这时，我们可以使用线性方程组求解函数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.4linalg工具包中的特殊矩阵

函数调用

所生成的矩阵

wronskian（f,v）

由变量v的向量函数f生成的朗斯基矩阵

bezout（p,q,x）

由变量x的多项式p、q生成的Bezout阵

hessian（expr,vars）

由变量向量vars的表达式expr生成的海赛阵

hilbert（n）

生成n阶希尔伯特阵

sylvester（p,q,x）

由变量x的多项式p、q生成的西尔维斯特阵

toeplitz（L）

由表达式有序表L生成的Toeplitz阵

4.3线性空间基本理论

4.4.1基本子空间

我们从线性代数中知道与矩阵A相关的Rn的四个基本子空间——行空间、列空间、化零空间和左零空间。

化零空间，实际上就是齐次线性方程组Ax=0的解空间。

我们可以利用上一节中介绍的函数linsolve（）来求得它的通解，在将通解重的辅助变量依次替换成1，就可以获得它的解空间的一组基。

对于左零空间，即方程组ATx=0的解空间，亦可以用同样的方法得到。

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

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

4.4.2正交基和Schmidt正交化

在欧氏空间中，我们可以定义两个向量的内积（innerproduct），在此基础上，我们还可以定义两个向量的夹角。

如向量α、β的夹角θ可以定义为

。

在Maple的linalg工具包中，有相应的函数可以计算向量的内积和夹角。

下面将通过例子来说明它们的用法。

在线性代数中，我们定义内积为零向量相互之间正交；并且，我们利用Schmidt正交化方法（Gram-Schmidtorthogonalizationprocess），由欧氏空间中一组普通基得到两两正交的基。

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-squaressolution）。

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的特征矩阵（characteristicmatrix），而行列式det（λI–A）展开得到的多项式称为特征多项式（characteristicpolynomial），在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函数连用，就可以生成具有约当标准型的矩阵。