线性代数

第三章 线性代数

1多项式

多项式是代数学中最基本的对象之一, 它不但与高次方程的讨论有关, 而且是进一步学习代数以及其它数学分支的基础.

1.1 多项式生成及类型测试

在Maple中, 多项式由名称、整数和其他Maple值, 通过+、-、*和^等组合而成. 例如：

> p1:=5*x^5+3*x^3+x+168;

p1 := 5x53x3x168

这是一个整系数单变量多项式. 多元多项式和定义在其他数域上的多项式可以类似构造：

> p2:=3*x*y^2*z^3+2*sqrt(-1)*x^2*y*z+2002;

p2 := 3xy2z32Ix2yz2002

由此可以看出, Maple中多项式的生成与“赋值”命令相似.

另外, 还可以通过函数randpoly生成随机多项式, 生成一个关于vars的随机多项式的格式如下：

randpoly(vars, opts);

其中, vars表示变量或者变量列表或集合, opts为可选项方程或者指定属性的名称. 如： > randpoly(x); ＃随机生成关于x的5次(默认)多项式

42x588x476x365x225x28

> randpoly([x, y], terms=8); ＃随机生成关于[x, y]二元多项式

78xy62x311x2y88x3yxy330y481x4y5x2y3

> randpoly([x, sin(x), cos(x)]);

73sin(x)cos(x)91sin(x)cos(x)2x3cos(x)5x4sin(x)86x2cos(x)343sin(x)4cos(x)

而要随机生成关于[x, y, z]的密集的、均匀的、度为2的多项式的命令为： > randpoly([x,y,z],dense,homogeneous,degree=2);

85x255zx37yx35z297yz50y2

用type命令可以测试多项式的类型：

> type(p1, polynom(integer, x)); #测试p1是否是一个关于x的整系数多项式

true

> type(p2, polynom(complex, {x, y, z})); #测试p2是否是一个关于{x, y, z}的复系数多项式

true

1.2 提取多项式系数

coeff函数用来提取一元多项式的系数, 而多元多项式所有系数的提取用命令coeffs, 指定系数的提取用命令coftayl.

(1) 提取多项式p中x^n的系数使用命令：coeff(p, x^n);或coeff(p, x, n); (2) 提取多项式p中变量x的所有系数并将相应的x幂存于变量t中：coeffs(p, x, ’t’); (3) 返回expr在x=a处的Taylor展式中(x-a)^k的系数: coeftayl(expr, x=a, k); > p:=2*x^2+3*y^3*x-5*x+68;

p := 2x23y3x5x68

> coeff(p, x);

3y35

> coeff(x^4-5*x^2-sin(a)*(x+1)^2, x^2);

5sin(a)

> s:=3*x^2*y^2+5*x*y;

s := 3x2y25xy

> coeffs(s);

5,3

> coeffs(s, x, 't');

5y,3y2

> t;

x,x2

> coeftayl(exp(x), x=0, 10);

1

3628800> p:=3*(x+1)^3+sin(Pi/3)*x^2*y+x*y^3+x-6;

1p := 3(x1)33x2yxy3x6

2> coeftayl(p, x=-1, 1);

13yy3

> coeftayl(p, [x, y]=[0, 0], [1, 0]);

10

返回默认为降序排列的多元多项式的首项和末项系数分别使用命令lcoeff、tcoeff： > lcoeff(p, x);

3

> tcoeff(p, x);

-3

1.3 多项式的约数和根

1.3.1多项式的最大公约因式()/最小公倍因式(lcm)

求多项式的最大公约因式/最小公倍因式的命令与求两个整数最大公约数/最小公倍数命令一样, 都是/lcm. 命令格式分别为：

(p1, p2, 't', 's'); lcm(p1, p2, 't', 's');

其中, 第3个参数t赋值为余因子p1/(p1, p2), 第4个参数s赋值为余因子p2/(p1, p2).

> p1:=x^4+x^3+2*x^2+x+1;

p1 := x4x32x2x1

> p2:=x^2+x+1;

p2 := x2x1

> (p1, p2, 't', 's');

x2x1

> t, s;

x21,1

> lcm(p1, p2);

(x21)(x2x1)

1.3.2多项式的平方根(psqrt)和第n次方根(proot)

求多项式p的平方根, 若不是完全平方, 则返回_NOSQRT：psqrt(p);

求多项式p的n次方根, 若不是完全n次方, 则返回_NOROOT：proot(p, n); > p:=x^4+4*x^3+6*x^2+4*x+1;

p := x44x36x24x1

> psqrt(p);

x22x1

> proot(p, 4);

x1

> proot(p, 8);

_NOROOT

1.3.3 多项式相除的余式(rem)/商式(quo)

计算p1除以p2的余式, 将商式赋值给q的命令格式为：rem(p1, p2, x, 'q'); 计算p1除以p2的商式, 将余式赋值给r的命令格式为：quo(p1, p2, x, 'r');

余式和商式满足：p1=p2*q+r, 其中degree(r, x) rem(x^5+x^3+x, x+1, x, 'q');

-3

> q;

x4x32x22x3

> quo(x^3+x^2+x+1, x-1, x, 'r');

x22x3

> r;

4

1.4 多项式转换及整理

1.4.1 将多项式转换成Horner形式

将多项式poly转换成关于变量var的Horner形式或者嵌套形式的命令格式如下：

convert(poly, horner, var); > convert(x^5+x^4+x^3+x^2+x+1, horner, x);

1(1(1(1(x1)x)x)x)x

> convert(x^3*y^3+x^2*y^2+x*y+1, horner, [x, y]);

1(y(y2y3x)x)x

1.4.2 将级数转换成多项式形式

将级数(series)转换成多项式(polynom)事实上就是忽略函数的级数展开式中的余项,

其命令格式为：

convert(series, polynom); > s:=series(sin(x), x, 10);

11511s := xx3xx7x9O(x10)

61205040362880> type(s, polynom);

false

> p:=convert(s, polynom);

11511p := xx3xx7x9

61205040362880> type(p, polynom);

true

1.4.3 将级数转换成有理多项式(有理函数)

将级数series(laurent级数或Chebyshev类型级数)转换成有理多项式(有理函数)ratpoly的命令格式为：

convert(series, ratpoly); > series(exp(x^2), x, 15);

11111011211x2x4x6x8xxx14O(x15)

26241207205040> convert(%, ratpoly);

161412xxx1120102

161412xxx11201021.4.4合并多项式系数(合并同类项)

将多项式具有相同次幂的项的系数合并在一起(包括正的、负的或者分数次幂), 即合并同类项(称为多项式的典范形式), 用命令collect:

collect(p, x);

collect(p, x, form, func); collect(p, x, func);

其中x是表示单变量x或多变量x1, x2, …, xn的一个列表或集合. > collect(a*ln(x)-ln(x)*x-x, ln(x));

(ax)ln(x)x > collect(x*(x+1)+y*(x+1), x);;

x2(1y)xy

> collect(x*(x+1)+y*(x+1), y);

x(x1)y(x1)

> p := x*y+a*x*y+y*x^2-a*y*x^2+x+a*x: collect( p, [x, y], recursive );

(1a)yx2((1a)y1a)x

> collect( p, [y, x], recursive );

((1a)x2(1a)x)y(1a)x

> collect( p, {x, y}, distributed );

(1a)x(1a)xy(1a)yx2

其中的参数recureive为递归式的，而distributed为分布式的。 1.4.5 将多项式(或者值的列表)排序

将多项式(或者值的列表)按升(或降)序次方排序时作命令sort. 命令格式： sort(L); sort(L, F); sort(A); sort(A, V);

其中, L—表示要排序的列表; F(可选项)—带两个参数的布尔函数; A—代数表达式; V(可选项)—变量

sort函数将列表L按升序次方, 将代数表达式A中的多项式按降序次方排列.

如果给定了F, 它将用来定义一个排序的列表. 如果F是符号“<”或者数值, 那么L是一个列表类型, 按数值型降序排列; 如果F是符号“>”, L按数值型升序排列; 如果F是一个词典编纂的符号, 那么字符串和符号列表将按词典编纂的次序排列. 另外, F必须是一个有两个参数的布尔函数.

在Maple中, 多项式并不自动按排序次序存储而是按建立的次序存储. 值得注意的是, sort函数对多项式的排序是破坏性的操作, 因为输入的多项式将会按排序后的次序存储.

用来对列表排序的算法是一种带早期排序序列监测的合并排序的递归算法, 而用于对多项式排序的算法是一个原位替换排序. > sort([3, 2, 1, 5, 4]);

[1,2,3,4,5]

> sort(x^3+x^2+1+x^5, x);

x5x3x21

> sort((y+x)^2/(y-x)^3*(y+2*x), x);

(xy)2(2xy) 3(xy)1.4.6 多项式的因式分解

一般情况下, 计算带有整数、有理数、复数或代数数系数的多项式的因式分解使用命令factor, 命令格式为：

factor(expr, K);

其中, expr为多项式, K为因式分解所在的区域 > factor(x^3+y^3);

(xy)(x2xyy2)

> factor(t^6-1, (-3)^(1/2));

1(2t1-3)(2t1-3)(2t1-3)(2t1-3)(t1)(t1) 16> factor(x^3+5, {5^(1/3), (-3)^(1/2)});

1(1/3)(1/3)(1/3)(1/3)(1/3)(2x5-35)(2x5-35)(x5) 4 但是factor不分解整数, 也不分解多项式中的整数系数, 整数或整数系数的分解使

用函数ifactor：

> ifactor(2^99+1);

(3)3 (19) (67) (683) (2420999355987) (20857) (5347)

1.5 多项式运算

1.5.1 多项式的判别式

多项式判别式在多项式求根中具有重要作用, 例如二次多项式的判别式. 求多项式p关于变量x的判别式的命令格式为：

discrim(p, x); > p1:=a*x^2+b*x+c;

p1 := ax2bxc

> discrim(p1, x);

4acb2

> p2:=a*x^3+b*x^2+c*x+d;

p2 := ax3bx2cxd

> discrim(p2, x);

27a2d218adcbb2c24b3d4ac3

1.5.2 多项式的范数

计算关于变元v的多项式p的第n阶范数用命令norm, 格式如下： norm(p, n, v);

其中p是展开的多项式, n为实常数或无穷大. 对于n≥1, 范数norm定义为norm(p, n, v)=sum(abs(c)^n), c=[coeffs(p, v)]^(1/n). > norm(x^3+x^2+x+1, 1);

4

> norm(x-3*y, 2);

10

> norm(x-3*y, infinity);

3

1.5.3 bernoulli数及多项式

bernoulli函数计算第n阶bernoulli数和关于表达式x的第n阶bernoulli多项式bernoulli(n, x), bernoulli(n, x)通过指数生成函数定义为：

textbernoulli(n,x)nt et1n0n!第n阶bernoulli数定义为：bernoulli(n)= bernoulli(n, 0) > bernoulli(6);

1 42> bernoulli(10, x);

53215x5x47x6x8x105x9 6622> bernoulli(6, 1/2);

-31 13441.5.4 用Bernstein多项式近似一个函数

bernstein多项式起因于stone-weierstrass近似分析理论, 该理论认为任何连续函数都可以在一个闭区间用多项式序列来近似, bernstein多项式就是该理论的一个实例. 若给定p:=(n, i, x)->binomial(n, i)*x^i*(1-x)^(n-i)时, Bernstein多项式定义为：

iBerstein(n,f,x)p(n,i,x)f()

ni0其中, binomial(n, r) = GAMMA(n+1)/GAMMA(r+1)/GAMMA(n-r+1)

在Maple中, 在[0, 1]区间上近似等于函数f(x)的关于x的n次方的Bernstein多项的操作命令是：

nbernstein(n, f, x);

其中, f必须是用一个程序或者操作符指定的单变量函数. > bernstein(3, x->1/(x+1), z);

3311zz2z3

41020> f:= proc(t) if t < 1/2 then 4*t^2 else 2 - 4*t^2 end if end proc:

bernstein(8, f, x);

17x63x5x2119x620x882x7 221.5.5 获取多项式的最高/最低次方

获取多项式的最高/最低次方的Maple命令分别为degree和ldegree, 格式如下： degree(p, x); ldegree(p, x); > p:=3/x^3+x^3-(x-1)^4;

p := 3> degree(p, x); ldegree(p, x);

134 x(x1)x34 -3

> degree(x*sin(x), x);

FAIL

> degree(x*sin(x), sin(x));

1

1.5.6 Euler数及多项式

第n阶Euler多项式E(n, x)由指数生成函数定义为：

2extE(n,x)nt teHn0n!据此, 第n阶Euler数E(n)由指数生成函数定义为：

2E(n)nt tteen0n!由此可见, 第n阶Euler多项式和Euler数的关系为：

1E(n)2nE(n,)

2在Maple中, 获取第n阶Euler多项式和Euler数的命令分别为：

euler(n, x); euler(n);

其中, n为非负整数, x为表达式. > euler(6);

-61

> euler(6, x);

x63x55x33x

> euler(6, 1/2);

-61 1.5.7 多项式插值

函数interp计算次方小于或等于n的过插值点(x[1], y[1]), (x[2], y[2]), …, (x[n+1], y[n+1])关于变量v的多项式, 命令格式如下：

interp(x, y, v);

值得注意的是, 如果相同的x值输入了2次, 不论是否输入了相同的y值都会导致错误, 也就是说所有插值必须不一样. 如过点(1, 2), (2, 4), (3, 6), (4, 8), (5, 10), (6, 12), (7, 14), (8, 16), (9, 18), (10, 22)的关于Z的插值多项式的计算如下： > interp([1,2,3,4,5,6,7,8,9,10], [2,4,6,8,10,12,14,16,18,22], z);

11829756301359452336515299z9zzzzzzzz2 18144040326046801134100812601.5.8 计算自然样条函数

计算一个关于变量z的次数是n(默认为3)的分段多项式来近似X, Y数据值. X值必须唯一且按升序排列, 而对Y值没有特别的约定. 命令格式为：

spline(X, Y, z, n);

> spline([1, 2, 3, 4], [2, 3, 6, 5], x, linear);#生成线性样条插值

1x33x9xx2x3 otherwise> spline([0, 1, 2, 3], [0, 1, 4, 3], x, cubic);＃生成三次样条插值

14xx35514414223xx2x555114151263xxx5555x1x2otherwise

1.6 有理式

所谓有理式, 是指可以表示成f/g形式的式子, 其中f与g均为多项式, 且g0.

1.6.1获取表达式的分子/分母

对于一个有理式x, 可以用numer(x)和denom(x)来获得它的分子(numerator)和分母(denumirator):

> f:=x^2+3*x+2: g:=x^2+5*x+6: h:=f/g;

x23x2h := 2

x5x6> numer(h);denom(h);

x23x2 x25x6

1.6.2 有理表达式的标准化

和数字除法不同, Maple并不自动化简分式, 使分子和分母互不可约, 除非分子和分母都表示成乘积的形式并且有公共部分. 如果要将有理式化简成最简形式即标准化，需要调用函数normal( ) (也称正则化). 事实上, normal函数提供了化简的一种基本形式, 它首先识别在有理数主域上等于0的表达式包括求和、乘积、整数幂次和变量构成的任何表达式. Normal的结果是使有理式转化为约化的标准形式, 即：分子分母是带整数系数的素数多项式. > normal(f/g);

x1 x3> normal(x^2-(x+1)*(x-1)-1);

0

> normal((f(x)^2-1)/(f(x)-1));

f(x)1

> normal(1/x+x/(x+1), expanded);

1xx2

xx2> normal(2/x+y/3);

16xy

3x1.6.3 将有理式转换为多项式与真分式和的形式

convert/parfrac可将有理函数 f 分解为多项式与真分式和的形式，这是有理函数积分的基础. 命令格式如下：

convert(f, parfrac, x, K);

其中, f为x的有理函数, x为主变量名称, K为可选参数(实数—real, 复数—complex, 扩展域—a field extended, 真—true, 假—false, 无平方项—sqrfree)等. > p:=(x^5+1)/(x^3+1);

x51p := 3

x1> convert(p, parfrac, x);

x2> convert(x/(x-a)^2, parfrac, x);

1x 2xx1a1 2xa(xa)1.6.4 将浮点数转换为接近的有理数

convert/rational将一个浮点数转换成一个近似的有理数, 转换的精度取决于环境变量Digits的值. 命令格式为：

convert(float, rational, Digits); > convert(1.333333333, rational);

4 3> convert(evalf(Pi), rational,6);

355 113> convert(evalf(Pi), rational, 3);

22 72 矩阵基本运算

2.1 数组和向量

在Maple中, 数组(array)由命令array产生, 其下标变量(index)可以自由指定. 下标由1开始的一维数组称为向量(vector), 二维以上的数组称为矩阵(matrix), 因此, 可以说向量与矩阵是数组的特殊形式. 生成数组与向量的命令分别列示如下：

array (m, n); ＃产生一个下标变量为m到n的一维数组

array (m, n, [a1, a2 .. ak]); ＃产生一维数组, 并依次填入初值

array ([a1, a2 ..ak]); ＃产生下标变量为1到k的一维数组, 并依次填入初值 vector(n, [v1, v2, …, vn]); ＃产生n维向量(n可省略) > A:=array(-1..2);

A := array(-1 .. 2,[])

> A[0]:=12;

A0 := 12

> A[1]:=2*x^2+x+4;

A1 := 2x2x4

> A;

A

> lprint(eval(A));

array(-1 .. 2, [(0)=12, (1)=2*x^2+x+4])

> B:=array(0..2, [x, x^2, x^3]);

B := array(0 .. 2,[(0)x(1)x2(2)x3])

> C:=array(1..3, [x, x^2, x^3]);

C := [x,x2,x3]

> print(C);

[x,x2,x3]

> type(B, vector);

false

> type(B, array);

true

> type(C, vector);

true

> type(C, array);

true

> vector(4, [1, x, x^2, x^3]);

[1,x,x2,x3]

> Vector(1..3, 5);

555

注意vector命令使用时大小写的区别：矢量vector(小写)指linalg包的基于数组的矢量, 而Vector(大写)指LinearAlgebra包的基于rtable的矢量.

当我们把一个向量赋给某个变量后, Maple在显示这个向量时, 只会显示出变量名称, 而不会显示出向量的内容, 这个设计可大大地简化向量的显示方式. 如果想要显示向量的内容, 可用evalm指令来强迫对向量求值. 事实上, evalm不仅可以显示向量, 也可以显示矩阵. > A:=array([x, y, z]);

A := [x,y,z]

> B:=array([1, 2, 3]);

B := [1,2,3]

> 2*A+B;

2AB

> evalm(2*A+B);

[2x1,2y2,2z3]

另一个有用的命令是将列表、数组或Vector矢量A转换成vector的函数convert： convert(A, vector);

> A:=array(1..2, 1..2, [[1, 2], [3, 4]]): convert(A, vector);

[1,2,3,4]

> V:=Vector([1, 6, 8]);

1V := 6 8> type(V, vector);

false

> convert(V, vector);

[1,6,8]

> type(%, vector);

true

2.2 矩阵的建立

在使用Maple进行线性代数运算时, 需要先载入线性代数工具包—linalg, 绝大部分线性代数运算函数都存于该工具包中.

在Maple中, 建立矩阵用matrix命令. 如建立一个m行n列的矩阵的命令格式为： matrix(m,n); matrix(m,n,init);

其中，init可以有很多种选择，如程序、数组、列表、方程组、代数表达式或者矩阵的初值等。

> with(linalg):

> A:=array([[1,2,3],[4,5,6],[7,8,9]]);

1A := 47> f:=(i,j)->i^2+j^3;

25836 9f := (i,j)i2j3

> B:=matrix(4,5,f);

25B := 1017912172428313365687380126129 134141> C:=matrix(2,2,[1,2,3,4]);

1C := 32 4 而函数diag可以生成块对角阵, band函数可以建立以常数对角元素的带状矩阵：

> diag(A,C);

1470025800369000001300 024> M:=Matrix([[1,1,1],[2,2,2,2],[3,3,3,3]], shape=band[1,3], scan=band[1,3]);

21M := 003210032100320003

内建于linalg的另一个有用的命令是从线性方程组中提取系数矩阵或增广矩阵, 实

现这一功能的函数是genmatrix, 其命令格式为：

genmatrix(eqns, vars, flag);

> eqns:={x+2*z=a,3*x-5*y=6-z};

eqns := {x2za,3x5y6z}

> A:=genmatrix(eqns,[x,y,z],'flag');

102aA := 3-516

> B:=genmatrix(eqns,[x,y,z]);

102B := 3-51 由此例可以看出, 若加参数'flag'则提取增广矩阵, 否则提取系数矩阵, 这一点对于

求解线性方程组很重要.

另外, Maple中还有一些预先编好的矩阵, 可以直接由函数(命令)获得： > hilbert(3);

> hilbert(3,x-1);

1121312131414x15x16x21233212131415 13x14x15x> toeplitz([1,2,3,4]);

123415x1 6x17x43 21> vandermonde([1,x-1,(x-1)^2,(x-1)^3]);

1111> fibonacci(3);

1x1(x1)2(x1)31111(x1)2(x1)4(x1)610111 03(x1) 6(x1)9(x1)1还有一个有用的命令是随机生成矩阵的函数randmatrix：

> randmatrix(3,3);

63459257-843-59-93 -62同样, 随机生成向量的函数为randvector.

2.3 矩阵的基本运算

矩阵的基本运算命令列于下表, 而矩阵的代数运算最直观的方法莫过于使用函数evalm了, 只要把矩阵的代数计算表达式作为它的参数, 就可以得到结果.

运算 加法 数乘 乘法 逆运算 转置 行列式 秩 迹 > A:=randmatrix(2,2); 函数调用 matadd(A, B) scalarmul(A, expr) multiply(A, B, …) inverse(A) transpose(A) det(A) rank(A) trace(A) 等效的运算符运算 evalm(A+B) evalm(A*expr) evalm(A &*B &* …) evalm(1/A)或evalm(A^(-1)) 无 无 无 无 7266A := -29-91

> B:=randmatrix(2,2);

-53-19B := -4768

> matadd(A,B); evalm(A+B);

1947-76-23 > multiply(A,B); evalm(A*B);

-691831205814-5637 > inverse(%); evalm(1/%);

-520-187969523623476181 -323-115311587273476181> evalf(%);

-.0002702678600-.0001494-.00027871845-.0003316858357 > transpose(A);

72-2966-91 > rank(B); trace(B);

2 15

值得注意的是, 矩阵乘法运算符“&*”有着和数乘“＊”、除法“/”相同的优先级, 在输入表达式时需要注意, 必要时要使用括号, 请看下面实验： > evalm(A&*1/A);

Error, (in evalm/amperstar) &* is reserved for matrix multiplication

> evalm(A&*(1/A));

&*()

上面最后的输出结果表示单位阵.

2.4 矩阵的求值

数组和单个的表达式不同, 具有数组类型的变量不会自动求值, 而需要用eval等函数的显式调用来强制地求值. 作为数组的特例, 矩阵和向量也是一样, 但要让矩阵最终完成求值还需使用函数map, map函数的主要功能是, 对于任意一元函数, 不加任何说明, 就可以作用在整个数组上, 得到的结果是每一个元素分别调用函数所得结果所组成的数组.

> with(linalg):

> Q:=matrix([[cos(alpha), sin(alpha)], [-sin(alpha), cos(alpha)]]);

cos()Q := sin()sin() cos()> alpha:=0;

 := 0

> eval(Q);

cos()sin()10sin() cos()01

> map(eval, Q);

另一个对矩阵中元素变量代换的命令是subs函数(但只能进行一层求值), 例如：

> P:=matrix(2, 2, [1, 2, x, x^3]);

1P := x> subs(x=8, eval(P));

2 3x182 5122.5 矩阵分解

矩阵分解在矩阵运算中的作用最明显的一点是可以大幅度提高运算效率. 矩阵分解算法较多，主要有LU分解、QR分解等。 2.5.1 LU分解

LUdecomp(A, P='p', L='l', U='u' , U1='u1', R='r', rank='ran', det='d');

其中, A为矩阵(长方形), P='p'—主元素因子, L='l'—单位下三角因子, U='u'—上三角因子, U1='u1'—修改的U因子, R='r'—行减少因子, rank='ran'—A的秩, det='d'—U1的行列式. 当然, 此命令的简写形式为：LUdecomp(A); > A:=vandermonde([1,3,5,7]);

11A := 1113571929127 125343>LUdecomp(A, P='p', L='l', U='u' , U1='u1', R='r', rank='ran', det='d');

1000> evalm(l);

12001880126 724800 011111> ran;

012300134

> d;

768

对于代数元素矩阵, 只有首项元素为0时才选主元素, 主元素矩阵返回为P. 2.5.2 QR分解

在QR分解中, 矩阵A被看作乘积Q•R, 此处, Q为正交或归一化矩阵, R为上三角阵. 该分解是对一系列线性无关向量作Gram-Schmidt正交化, 因此, Q包含了正交的向量, R的列记录了产生原向量的线性组合. 命令格式为：

QRdecomp(A, Q='q', rank='r', fullspan=value);

其中, A为矩阵(长方形), Q='q'—A的Q因子, rank='ran'—A的秩, fullspan=value—在Q中包含空格(true或false).

> A:=matrix(3,3,[1,2,3,4,5,6,7,8,10]);

1A := 476600> r;

25836 109766665 1111166> QRdecomp(A,Q='q',rank='r');

1366113111103

> evalm(q);

1666623366766663111111111111111661 63166另外, 对实正定矩阵可采用Cholesky分解, 这是一种形式更对称的LU分解.

3 矩阵的初等变换和线性方程组求解

3.1 矩阵的初等变换

矩阵的初等变换是各种消去法的基础, 是求解线性方程组的基础. 对于符号运算, 有时候我们不仅要求最后的求解结果, 而且要求中间的求解步骤, 此时就需调用线性代数工具包中的初等变换函数. 如下表所示：

行变换 列变换 函数调用 mulrow(A, r, expr) addrow(Ar1, r2, m) swaprow(A, r1, r2) mulcol(A, c, expr) addcol(A, c1, c2) swapcol(A, c1, c2) 对应的初等变换 用标量expr乘以矩阵A的第r行 将矩阵A的第r1行的m倍加到第r2行上 互换矩阵A的第r1行和第r2行 用标量expr乘以矩阵A的第c列 将矩阵A的第c1列的m倍加到第c2列上 互换矩阵A的第c1列和第c2列 除了这些初等变换外, Maple在linalg工具包中还提供了一些矩阵的形状变换, 可以在利用初等变换解决问题时起到辅助作用, 我们也将它们一并列出.

函数调用 contat(A, B, …); stackmatrix(A, B, …) row(A, i) col(A, i) row(A, i .. k) col(A, i .. k) delrows(A, i .. k) 所作的变换或操作 将矩阵A, B, …在水平方向上合并成一个矩阵 将矩阵A, B, …在竖直方向上合并成一个矩阵 取矩阵A的第i行 取矩阵A的第i列 取矩阵A的第i到k行 取矩阵A的第i到k列 删除矩阵A中i到k行剩下的子矩阵 delcols(A, i .. k) extend(A, m, n, x) submatrix(A, m .. n, r .. s)

删除矩阵A中i到k列剩下的子矩阵 扩展矩阵m行n列并用x填充 取矩阵A的m到n行、r到s列组成的子矩阵 例：利用矩阵的初等变换判断矩阵A是否可逆, 若可逆, 求A-1

213A011

120> with(linalg):

A:=matrix([[-2, 1, 3], [0, -1, 1], [1, 2, 0]]);

-2A := 01-2AI := 01> mulrow(AI, 1, -1/2);

1-1231 010001000 1> AI:=concat(A, array(identity, 1..3, 1..3));

1-12310101> addrow(%,1, 3, -1);

-12-12-3210-12000100 01100> mulrow(%, 2, -1);

-12-152-32132-1201201000 1100> addrow(%, 2, 1, 1/2);

-12152-32-132-2-12012-10-10-100 101001052> addrow(%, 2, 3, -5/2);

100100> mulrow(%, 3, 1/4);

100100> addrow(%, 3, 1, 2);

100100> addrow(%, 3, 2, 1);

100100> AA:=submatrix(%, 1..3, 4..6);2-103122-2-12-10412-2-12-101180-14-101180-140181182-10 01-120-105 21-120-1051 843142-1051 843142-31845184

AA := -14181834-3858121414 上述判断矩阵可逆并求逆矩阵的方法即高斯—约当消去法, 从例中可以看出, 按高斯—约当消去法步骤逐一在Maple下操作即可实现. 在这里, 数学基础显然是最重要的. 当然, 上述例子可通过下述两句命令即可完成： > det(A);

若det(A)不为0执行下句即可求出A的逆矩阵： > inverse(A);

3.2 线性方程组的解

线性方程组求解常常转化为与其等价的矩阵问题来解决的. Maple中可借助函数genematrix及高斯消去法函数gausselim联合完成：

> eqns:={x+2*y+3*z=a,8*x+9*y+4*z=b,7*x+6*y+5*z=c};

eqns := {x2y3za,8x9y4zb,7x6y5zc}

> A:=genmatrix(eqns,[x,y,z],'flag');

1A := 87> gausselim(%);

296345ab c1002-703-20487ab8a 158cab77可以看出, 在高斯消去法的结果矩阵中出现了许多分数, 这使表达式看起来不美观, 而且随着矩阵的增大, 分数的分母会越来越大, 如果矩阵中含有符号, 结果矩阵中将出现分式. 为了避免出现上述情况, 可采用无分式高斯消去法(fraction-free Gaussian elimination), 相应的命令为ffgausselim： > ffgausselim(%%);

1002-703a -20b8a-487c15a8b 在把系数矩阵转化成相应的上三角阵后, 就可以用回代(back substitution)的方法求

出各个未知数的值—backsub： > backsub(%);

7a1b19c,1b1a5c,7c5a1b 168341248166和回代相对, 也有前代(forward substitution)的概念, 可用前代函数forwardsub法语解下三角矩阵. 有了这一对函数, 我们就可以用LU分解来线性方程组了, 相应的命令函数为LUdecomp. 关于上一问题也可以利用高斯—约当消去法(Gauss-Jordan elimination), 把系数矩阵变换成单位阵直接得到结果： > gaussjord(A);

1000100017119abc168115

bac3412751cab48166上面介绍的实际上都是线性方程求解的中间步骤, 事实上我们可以一步到位, 不管它到底用的是什么方法, 而只要求得最终的结果. 这方面, Maple提供了一个非常好的函数linsolve, 这个函数不仅可以用来求解具有唯一解的线性方程组, 而且可以求解有无穷多解的方程组并给出通解. 试看下面的实验：

x13x23x32x41例：求解线性方程组：2x16x29x35x44

x3x3x13231> A:=matrix([[1, 3, 3, 2], [2, 6, 9, 5], [-1, -3, 3, 0]]);

1A := 2-1> B:=vector([-1, 4, 13]);

36-339325 0B := [-1,4,13]

> linsolve(A, B);

[133_t23_t1,_t2,_t1,63_t1]

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

3.3 最小二乘法求解线性方程解

在实际问题中, 由于误差或者其他各方面的原因, 很容易出现无解的方程组. 这样问题的解决则变成求出一个“最优”的近似解. 在线性代数中, 通常采用最小二乘解(least-squares solution). linalg中对应的函数是leastsqrs, 它有两种调用格式：

leastsqrs(A, b); leastsqrs(S, v);

其中，A为矩阵, b为向量，S为线性方程组, v为变量名。 > with(linalg):

> S:={c[0]+c[1]+c[2]-3, c[0]+2*c[1]+4*c[2]-10, c[0]-9/10, c[0]-c[1]+c[2]-3};

9S := {c0c1c23,c02c14c210,c0c1c23,c0}

10> leastsqrs(S,{c[0],c[1],c[2]});

715991{c1,c0,c2}

20020040 再如：

> A:=matrix(3,2,[0,1,1,1,1,1]);

0A := 11> b:=vector([0,1,-1]);

11 1b := [0,1,-1]

> leastsqrs(A,b);

[0,0]

对于相容方程组，leastsqrs命令与solve功能一样：

> leastsqrs({x+y=2,y+z=3,x+z=3},{x,y,z});

{x1,y1,z2}

> solve({x+y=2,y+z=3,x+z=3},{x,y,z});

{x1,y1,z2}

3.4 正定矩阵

在Maple中, 判断一个矩阵是否是正(负)定矩阵非常简单, 只需要使用命令definite

即可, 而且还可求出符号矩阵的正(负)定条件： > with(linalg):

> A := matrix(2, 2, [2,1,1,3]);

2A := 1true

1 3> definite(A, 'positive_def'); > B:=array(1..2, 1..2, 'symmetric'):

definite(B, 'negative_semidef');

B1,10 and B1,1B2,2B1,20 and B2,20

由上述实验可以看出, definite的第1个参数是需要判定的矩阵, 第2个参数是矩阵的正定性, 共有4种情况：'positive_def'(正定)、'positive_semidef'(半正定)、'negative_def'(负定)、'negative_semidef'(半负定). 当判定数值矩阵时, 返回布尔值true/ false; 判定符号矩阵时, 它返回一个布尔表达式, 表示正负定的条件.

24 特征值、特征向量和矩阵的相似

4.1 矩阵的相似

在Maple中可以利用linalg中的issimilar判断两个矩阵是否相似. 命令格式如下： issimilar(A, B, 'P');

其中A, B为方阵, P为转换矩阵(可选). 如果A, B相似, 则返回true, 否则返回false. > with(linalg): > A:=matrix(2, 2);

A := array(1 .. 2,1 .. 2,[])

> B:=matrix(2, 2);

B := array(1 .. 2,1 .. 2,[])

> issimilar(A&*B, B&*A);

true

在这个例子中, Maple作了一定的假设, 由于A的行列式是一个符号表达式, Maple在判定时假设它为非零常数. 不过更多情况下, 这个函数是用在数值矩阵的判定上, 同时也还可以求出转换矩阵P, P赋值为满足A=inverse(P)*B*P的转换矩阵P. > with(linalg,matrix,issimilar,eigenvalues,diag): > A:=matrix(3,3,[1,2,3,4,5,6,7,8,9]);

1A := 47> B:=diag(eigenvalues(A));

25836 90 1533322000B := 0> issimilar(A,B,P);

> print(P);

015333220true

1-116365711117 33333312132619812396571111733333312132619812396> map(normal,evalm(P^(-1)&*B&*P));

14725836 94.2 特征值和特征向量

特征值问题可以说是最常见的线性代数问题了. 在Maple中求解特征矩阵(characteristic matrix)和特征多项式(characteristic polynomial)的函数分别是charmat和charpoly, 求解相应的特征值和特征向量的函数分别是eigenvalues(eigenvals)和eigenvectors. 命令格式如下：

charmat(A, lambda); charpoly(A, lambda); eigenvalues(A); eigenvectors(A);

下面通过几个实例学习上述4个函数的用法：

> A:=randmatrix(2,2);

1772A := -99-85

> charmat(A,lambda);

1799> charpoly(A,lambda);

-72 852685683

> eigenvalues(A);

343I503,343I503

> eigenvectors(A);

343I503,1,{1,171I503},343I503,1,{1,171I503} 24242424> A:=matrix(3, 3, [1, 2, 2, 2, 1, 2, 2, 2, 1]);

1A := 22> charmat(A, lambda);

21222 1-2-2

11-2-2> charpoly(A, lambda);

-21-233295

> eigenvals(A);

5,-1,-1

> eigenvectors(A);

[-1,2,{[-1,1,0],[-1,0,1]}],[5,1,{[1,1,1]}]

可以看到, 特征向量函数eigenvectors在给出特征向量的时候, 还给出了特征值及特征值重数. 一般情况下, 如果同时需要得到特征值和特征向量用eigenvectors函数即可. 另外, 这个函数还可处理符号矩阵, 此时, 通常以根式的形式给出结果, 如果加上可选参数'implicit', 结果就以RootOf的形式给出.

函数eigenvalues不仅可以求解普通特征值问题, 还可求解广义特征值问题, 即求解方程det(CA)0相应的命令函数为eigenvals(A, C).

由特征值理论, 任何一个实对称阵, 都可以用求特征值的方法将它对角化, 但对于非对称阵或者是复矩阵, 就没有这样的保证了. 不过对于任意复矩阵A, 还是可以把它化成相似的Jordan标准型, 相应的命令函数是jordan, 而且还可给出转换矩阵P. > A:=matrix(3, 3, [2, -1, 1, 2, 2, -1, 1, 2, -1]);

2A := 21> jordan(A, 'p');

-1221-1 -101 1100> print(p);

11003312110 0在Maple中, 可以调用函数JordanBlock(,k)构造约当块J(,k), 利用该命令和diag函数连用, 就可以生成具有约当标准型的矩阵: > JordanBlock(3, 4);

3000> linalg[JordanBlock](x, 5);

1300013000 13x00001x00001x00001x00001x

5 线性空间基本理论

5.1 基本子空间

与矩阵A相关的Rn的四个子空间分别是行空间、列空间、化零空间和左零空间. 化零空间, 实际上就是齐次线性方程组AX=0的解空间, 可利用linsolve求得其通解, 再将通解中的辅助变量依次替换成1, 就可以获得解空间的一组基. 对于左零空间, 即方程组ATX=0的解空间, 亦可用同样的方法得到. 行空间和列空间, 在linalg中也有专门的函数rowspace和colspace, 可以获得它们的基和维数. 例举如下： > with(linalg):

> A:=matrix(3, 2,[2, 0, 3, 4, 0, 5]);

2A := 30> rowspace(A, 'dim');

04 5{[0,1],[1,0]}

> dim;

2

> colspace(A);

5-15{0,1,,1,0,} 48函数rowspace和colspace的第二个参数是可选参数, 用来返回行空间或列空间的维

数, 它必须是一个变量名—可以是一个未赋值的新变量, 也可以是已有值的变量, 但要加上延迟求值符“'”.

而寻找向量空间的基的函数是basis: > v1:=vector([1,0,0]):

v2:=vector([0,1,0]): v3:=vector([0,0,1]): v4:=vector([1,1,1]): basis({v1,v2,v3});

{v2,v1,v3}

> basis({vector([1,1,1]),vector([2,2,2]),vector([1,-1,1]),vector([2,-2,2]),vector([1, 0, 1]), vector([0, 1, 1])});

{[1,1,1],[0,1,1],[2,-2,2]}

5.2 正交基和Schmidt正交化

在欧氏空间中, 我们可以定义两个向量的内积(inner product), 在此基础上, 我们还可以定义两个向量的夹角. 如向量、的夹角可以定义为：

arccos,

在Maple的linalg工具包中, 有相应的函数innerprod、angel可以计算向量的内积和夹角. 如：

> alpha:=vector([1, 2, -1, 11]);

 := [1,2,-1,11]

> beta:=vector([2, 3, 1, -1]);

 := [2,3,1,-1]

> innerprod(alpha, beta);

-4

> angle(alpha, beta);

4 arccos127151905求三维向量的向量积使用命令crossprod: > alpha:=vector([1, 2, 1]);

 := [1,2,1]

> beta:=vector([2, 3, 1]);

 := [2,3,1]

> crossprod(alpha,beta);

[-1,1,-1]

定义内积为零的向量间正交, 此时可以利用Schmit正交化方法, 由欧氏空间中一组普通基得到两两正交的基. Maple中的相应函数是GramSchmidt, 它的输入参数是由一组向量组成的集合(或有序表), 将给出Schmidt正交化后的向量集合(或有序表). 输入的向量必须是线性无关的, 否则, 结果向量之间也将线性相关. 函数并不对向量进行单位化. 如果需要得到一组正交标准基, 还需要用map方法对这些向量使用单位化函数normalize.

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

1,1,1,2T, 5,8,2,3T, 3,9,3,9T

> A:={vector([1, 1,-1, -2]), vector([5, 8, -2, -3]), vector([3, 9, 3, 9])};

A := {[5,8,-2,-3],[3,9,3,9],[1,1,-1,-2]}

> simplify(map(normalize, GramSchmidt(A)));

23113,{291,291,291,2919118218291

17,17,17,27,239,539,139,139}739777393913事实上, 向量的Schmidt正交化过程实际上给出了矩阵的QR分解.