线性代数 Maple实验 梁燕来,胡源艳.doc

附录 Maple数学实验
Maple是目前世界上最为通用的数学和工程计算软件之一,Maple系统内置高级技术解决建模和仿真中的数学问题,包括世界上最强大的符号计算、无限精度数值计算、创新的互联网连接、强大的4GL语言等,内置超过5000个计算命令,数学和分析功能覆盖几乎所有的数学分支,如微分方程、特殊函数、线性代数、图像声音处理、统计和动力系统等.
Maple不仅仅提供编程工具,更重要的是提供数学知识.Maple都可以帮助您快速、高效地解决类似从简单的数字计算到高度复杂的非线性问题等,用户通过Maple产品可以在单一的环境中完成多领域物理系统建模和仿真、符号计算、数值计算、程序设计、技术文件、报告演示、算法开发、外部程序连接等功能,满足各个层次用户的需要,从高中学生到高级研究人员.
安装Maple软件后,启动Maple进入默认的用户主界面,该界面包括两个主要的窗口:命令窗口（Command directory）,工具窗口（Tool directory）,如图5.1所示.

15

图5.1 用户主界面

图5.2 命令窗口

命令窗口是和其编译程序连接在一起的主要窗口,点击工具栏中“[>”提示符后, Maple就处于准备状态,在提示符后输入正确的的运算表达式后,只需按工具栏的“！”工具(或Enter键),命令窗口就会显示运算结果,如图5.2所示.
Maple语言和其他高级程序设计语言一样，都具有类似的基本程序结构，不管多复杂的计算都是由一系列语句组成,这些语句有的是顺序依次执行的,有的需在特定的地方循环、选择、跳转执行.下面逐一介绍这些基本结构.
1．顺序结构
顺序结构，是程序语句书写的先后顺序，是流程控制语句中最简单的一类。
2．选择结构
所谓选择结构,常称为分支结构,是指在程序执行时,依据不同的条件,从两个或多个程序块中选择其中一个执行.if条件结构的语法为:
if 条件1 then 程序块1；
elif 条件2 then 程序块2；
else 程序块3；
fi;
其中,fi是条件语句的结束标志,也可以写作end if.
此语句是一个三分支结构.若条件1为真,则执行程序块1;若条件1为假,条件2为真,则执行程序块2;若条件1和条件2都为假,则执行程序块3.
例如,判断两个数中较大者:
>bigger:=proc(a,b)
if a>=b then a;
else b;
fi;
end;
>bigger(3,5); #函数调用
运行以上程序后结果为:5
3．for循环
在程序设计中,常常需要把相同或者类似的语句连续执行多次,此时,可以通过for循环结构更便捷地编写程序.
例如,求1至5自然数的和:
>total:=0;
for i from 1 to 5 do
total:=total+i;
od; #循环体的结束标志od(也可以写作end do).
其中,i为循环变量,用于控制循环次数；do与od之间的一条或多条语句,称为循环体,是要反复执行的语句. for循环结构的语法可总结为:
for 循环变量 from初值 to终值 [by步长] do 循环体 od
在for循环结构中,初值和终值一般是整数,若初始值或步长是1则都可以省略.
4．while循环
for循环是在那些已知循环次数,或者循环次数可以用简单表达式计算的情况下比较适用.但有时循环次数并不能简单地给出,要通过计算,判断是否继续循环,这时,可以使用while循环.while循环标准结构为:
while条件表达式do 循环体 od
Maple先判断条件是否成立,若成立就一遍遍执行循环体,直到条件不成立为止.
例如,辗转相除法计算两个自然数的最大公约数.
>GCD:=proc(a::posint,b::posint)
local p,q,r;
p:=max(a,b); q:=min(a,b); r:=irem(p,q);
while r<>0 do
p:=q; q:=r; r:=irem(p,q);
od;
q;
end;
>GCD(6,9);
运行以上程序后结果为: 3
5．递归子程序
正如在一个子程序中可以调用其他的子程序一样(比如系统内部函数),一个子程序也可以在它的内部调用它自己,这样的子程序称为递归子程序.
例如,递归定义Fibonacci数列:
>fib:=proc (a::posint)
local i,x;
if a=1 or a=2 then x:=1;
else x:= fib(a-1)+ fib(a-2); #递归调用
fi
end;
>fib(10); #调用该函数，输出Fibonacci数列的第10项
运行结果为:55
实验一 矩阵的输入与特殊矩阵生成
1．矩阵输入
建立矩阵的方法主要有:
1) matrix(m,n,[a11,a12,…a1n,a21,…a2n,…,am1,…,amn])
例6.1.1 建立一个2行2列矩阵A=.
>A:=matrix(2,2,[1,2,3,4]);
执行以上指令后命令窗口显示为:
2) array(1..m,1..n,[[a11,…,a1n],…,[am1,…,amn]])
例6.1.2 建立矩阵B=.
>B:=array([[sin(alpha),cos(beta)],[sin(beta),cos(alpha)]]);

例6.1.3 输入3行3列矩阵C,使其元素为.
>f:=(i,j)->i^2+j^3;
>C:=matrix(3,3,f);

2．矩阵元素操作
完成矩阵元素输入后,对矩阵的操作主要包括矩阵的合并、矩阵部分元素删除、矩阵元素提取、矩阵元素的扩充等,下面作简单介绍.
1) 取子阵
submatrix(A,m..n,r..s):取矩阵A的m到n行、r到s列组成的子阵.
例6.1.4 在上例的矩阵C中取1到2行、1到2列得到子阵E.
>with(linalg);
>E:=submatrix(C,1..2,1..2);

2) 取行列
row(A,i):取矩阵A的第i行.
row(A,i..k):取矩阵A的第i到k行.
col(A,i):取矩阵A的第i列.
col(A,i..k):取矩阵A的第i到k列.
3) 矩阵扩展
例6.1.5 对例6.1.1中的矩阵A扩展2行2列并用0填充.
> extend(A,2,2,0);

4) 矩阵合并
例6.1.6 将矩阵A和矩阵B合并成一个矩阵.
>concat(A,B); #或执行<A|B>,水平方向上合并.

>stackmatrix(A,B); #竖直方向上合并.

5) 删除行列
delrows(A,i,k):删除矩阵A中i到k行剩下的子矩阵.
delcols(A,i,k):删除矩阵A中i到k列剩下的子矩阵.
3．矩阵初等变换
矩阵的初等变换是各种消去法的基础,是求解线性方程组的基础.此时就需调用线性代数工具包linalg中的初等变换函数.
1) 行(列)交换
swaprow(A,i,j):互换矩阵A的第i行和第j行.
swapcol(A,i,j):互换矩阵A的第i列和第j列.
>with(linalg);
>A:=matrix(3,3,[1,2,3,4,5,6,7,8,9]);

>swaprow(A,2,3);

>swapcol(A,2,3);

2) 行(列)数乘
mulrow(A,r,expr):用标量expr乘以矩阵A的第r行.
mulcol(A,c,expr):用标量expr乘以矩阵A的第c列.
>A:=matrix(3,3,[1,2,3,4,5,6,7,8,9]);
>mulrow(A,2,5);

>mulcol(A,2,6);

3) 行(列)倍加
addrow(A,r1,r2,m):将矩阵A的第r1行的m倍加到第r2行上.
addcol(A,c1,c2,m):将矩阵A的第c1列的m倍加到第c2列上.
>addrow(A,1,2,-2);

4) 化阶梯形
gausselim(A):用高斯消元法化矩阵A为行阶梯形,结果为上三角阵.
gausselim(A, 'r'):化A为阶梯形后要求返回A的秩r.
gausselim(A,'r','d'):要求返回A的秩r,同时返回A的行列式d.
rref(A):化矩阵A为最简行阶梯形矩阵.
例6.1.7 将矩阵A=化为最简行阶梯形矩阵.
>matrix(3,4,[4,5,8,7,2,6,5,4,5,2,1,8]);
>B:=rref(A);

4．特殊矩阵的生成
adj(A):生成A的伴随矩阵.
DiagonalMatrix (V,m,n): 生成m行n列对角矩阵,其中V是一个向量.
>DiagonalMatrix([<1,2>,3,<<4,5>|<6,7>>],6,6);

>diag(1,2,3,4,5); #以1,2,3,4,5为对角元的对角矩阵
fibonacci(m):生成m阶fibonacci矩阵.
HilbertMatrix(m,n,x):生成m行n列希尔伯特矩阵.
>HilbertMatrix(2,3,x)

IdentityMatrix(m):生成m单位阵,或array(identity,1..m,1..m).
>with(linalg);
>IdentityMatrix(3);

randmatrix(m,n):生成一个m行n列的随机矩阵.
VandermondeMatrix(V):生成范德蒙矩阵,其中V是一个向量.
ZeroMatrix(m,n):生成m行n列的零矩阵，或array(1..m,1..n,sparse).
实验题目
1. 输入矩阵A=,并提取矩阵的第2行和第3行元素.
2．生成一个的随机矩阵.

实验二 矩阵的运算
1．矩阵的代数运算
如果已经输入矩阵A和B,则可采用下列命令进行计算.
运算
函数
等效的函数
加法
matadd(A,B)
evalm(A+B)
数乘
scalarmul(b,expr)
evalm(b*expr)
乘法
multiply(A,B,…)
evalm(A&*B&*…)
逆运算
inverse(A)
evalm(1/A)或evalm(A^(-1))
转置
transpose(A)
无
例6.2.1 设A=,B=,求A+B、AB、A-1B、A’.
>A:=matrix(3,3,[2,3,-3,2,0,11,-1,3,4]);
>B:=matrix(3,3,[-1,0, 1,0,2,2, 3,5,1]);
>matadd(A,B);

>multiply(A,B);

>multiply(inverse(A),B);

>transpose(A);

2．矩阵的特征参数运算
在科学计算时,常用到矩阵的特征参数,在Maple中计算命令为：rank(A)求矩阵A的秩,trace(A)求矩阵A的迹, det(A)求矩阵A的行列式,cond(A)求矩阵A的条件数.
实验题目
1. 设A=, B=,求A+B、AB、A2、BA-1、A’.
2. 求向量组(2,-3,2,-1)T,(2,3,-1,4)T,(3,1,1,-2)T,(2,-2,0,3)T,(2,1,-2,1)T的秩.

实验3 线性方程组的求解
利用Maple软件可以轻松实现线性方程组的求解.其过程如下:首先利用Maple软件线性代数函数包linalg中的函数genmatrix,从线性方程组中生成系数矩阵或增广矩阵,然后用函数gausselim将矩阵化为行阶梯矩阵,进一步用函数rref将其化为简化行阶梯矩阵,最后用函数backsub执行回代任务,即可求出线性方程组的解向量.
例6.3.1 求线性方程组的全部解.
>with(linalg);
>eq:={x1-x3+x5=1,x2-x4+x6=2,x1-x2+x5-x6=3,x2-x3+x5=4,x1-x4+x5=5};
>A:= genmatrix(eq, [x1, x2, x3, x4, x5, x6],’flag’);

>B:=gausselim(A); C:=rref(B);

>backsub(C,false,’t’);

可以看到, Maple用辅助变量t1、t2给出了方程组的通解:x6=-1+t2-t1, x5 =2+t2-t1,
x4= 4+t2, x3=t2, x2=6+t1, x1=t1 (Maple显示形式为[-1+t2-t1 2+t2-t1 4+t2 t2 6+t1 t1 ]).
实验题目
1．将矩阵A=化为最简行阶梯形矩阵.
2．求线性方程组的全部解.
实验4 特征值、特征向量与特征多项式
方阵的特征值和特征向量在对角化中和微分方程组等问题中有着广泛的应用.Maple中的命令可以直接求得特征值、特征向量和特征多项式.

函数
说明
特征值
eigenvalues(A)
或eigenvals(A)
求矩阵A的特征值,加radical显根号形 式,implicit显复数形式.
特征向量
eigenvectors(A)
返回了特征向量、特征值及其重数.
特征矩阵
charmat(A,lambda);

特征多项式
charpoly(A,lambda)

正定性
definite(A,option)
判定矩阵A的正定性,参数option有4种情况:正定(positive_def),半正定(positive_semidef),负定(negative_def)和半负定(negative_semidef).判定数值矩阵时,返回布尔值true/false;判定符号矩阵时,它返回一个布尔表达式,表示正负定的条件.
相似性
issimilar(A,B,'P')
判断两个方阵A,B是否相似,P为变换矩阵(可选),如 A,B相似,则返回true,否则返回false.
例6.4.1产生随机矩阵,计算其特征矩阵、特征多项式、特征值和特征向量.
>A:=randmatrix(2,2);

>charmat(A,lambda);

>charpoly(A,lambda);

>eigenvalues(A);

>eigenvectors(A);

例6.4.2 将方阵B=对角化.
>B:=matrix(3,3,[1,2,2,2,1,2,2,2,1]);

>C:=eigenvals(B);

>diag(C);

实验题目
1、对于方阵A=,计算其特征矩阵、特征多项式、特征值和特征向量.
2、将方阵A=对角化.
实验5 综合练习
例6.5.1 求线性方程组的全部解.
>with(linalg);
>eq:={x1+x3+x4=1,2*x1-3*x3+x4=2,x1+x2-x3+x4=1};
>A:= genmatrix(eq,[x1,x2,x3,x4],’flag’);

>B:=gausselim(A);

>C:=rref(B);

>backsub(C,false,’t’);

可以看到,Maple用辅助变量t1给出了方程组的通解:x4=1-4t1/5, x3=-2t1/5, x2=-t1/5, x1= t1.
学习中发现,手工将一个实对称高阶矩阵正交化是较为麻烦的.利用Maple软件可以辅助解决这个问题.首先利用Maple函数charpoly(A)求A的特征多项式,再通过函数eigenvector(A)求出矩阵的特征向量.进一步用函数gramschmidt将特征向量Schmit正交化,得到A的n个两两正交的特征向量.下例说明了矩阵正交化过程.
例6.5.2 将矩阵A=正交化.
>with(linalg); with(LinearAlgebra);
>A:=array([[0,-1,1],[-1,0,1], [1,1,0]]);
>eigenvectors(A);

可以看到，函数eigenvectors()返回了矩阵A的特征向量<-1,1,0>,<1,0,1>,<-1,-1,1>、特征值及其重数、.
>Q:=GramSchmidt([<-1,1,0>,<1,0,1>,<-1,-1,1>],normalized);

>Q:=array([[-sqrt(2)/2,sqrt(6)/6,-sqrt(3)/3],[sqrt(2)/2,sqrt(6)/6,-sqrt(3)/3],[0,sqrt(6)/3,sqrt(3)/3]]);

> multiply(inverse(Q),A,Q);

即Q-1AQ= ..
=
例6.5.3利用顺序主子式判定二次型的正定性.
>with(linalg); with(LinearAlgebra); with(MTM);
v:=vector([0.5,0.5,0.5,0.5,0.5,0.5]);
A:=diag(v,1);
B:= diag([1,1,1,1,1,1,1]);
C:= diag(v,-1);
E:=matadd(C, matadd(A, B));

>flag:=0;
for i from 1 to 7 do
M:=submatrix(E,1..i,1..i);
if det(M)<=0 then
flag:=1;break;
fi;
od;
if flag=0 then printf(“正定二次型”);
else printf(“非正定二次型”);
fi;
运行结果：正定二次型
实验题目
1. 求非齐次线性方程组的通解.
2. 判定二次型的正定性.