查看“Maxima在线性代数的应用”的源代码

作者:蔡炎龍

原文：[http://math.nccu.edu.tw/~yenlung/mynotes/maximalinear_html html]

简体整理:dbzhang

本文：[http://maxima.sourceforge.net/docs/tutorial/cn/maxima_l_cn.pdf pdf]

状态：整理中...

==简介==
这篇文章,是介绍Maxima 这套数学软件,在学习线性代数的应用。

Maxima 是一个所谓的「电脑代数系统」(Computer Algebra System, CAS),这种系统比较为人熟知的还有Mathematica 和Maple 等等。我们选定Maxima 做为我们使用的程序,主要有三个原因:

;免费 : Maxima 是免费,又是各平台都有的。所有的人可以在自己的电脑上练习。
;功能完整 : Maxima 虽然不要钱,并不代表不好。Maxima 不论计算或图形功能都十分完整。事实上,Maxima 是最早的全功能CAS 系统Macsyma 的后代。
;具代表性 : 许多新的CAS 系统,如Maple, Mathematica 都多少受到Macsyma 的启发。所以学会Maxima,要学会Maple 或Mathematica 等软件都是很容易的事。

这篇文章主要是介绍线性代数相关功能。我们不假设同学已会基本的Maxima 使用方式,所以我们会用到的概念,也许不纯粹是线性代数的,也会一并介绍。专就线性代数而言,我们要会的其实并不多。想要快速进入状况,可以跳过前面的部份,直接看线性代数相关指令,在操作上有问题时,再回头看有问题部份的相关说明即可。

如果同学们比较喜欢使用Mathematica,Maple,或是Matlab 等商业软件也是可以的。我们系上的电脑室有提供这些软件,可以上机试用看看。

== 基本概念 ==

我们先介绍一下Maxima 操作的方式。

=== Maxima 当计算器 ===
我们先来看,如果我们要把Maxima 当计算器用,会是什么情况?

 (%i1) 1+1;                                                                   
 (%o1)                                  2
 (%i2) 3*4*7;
 (%o2)                                 84
 (%i3) 9/3;
 (%o3)                                  3

到目前为止,似乎还没什么特别。除了可以做复杂一点点的运算,和平常的计算机或数值计算软件也没什么不同。以下的例子就不一样了:

 (%i4) 7/3;
                                        7
 (%o4)                                  -
                                        3
 (%i5) 1/2+2/3;
                                        7
 (%o5)                                  -
                                        6

从(%o4)我们看到,7/3这种运算,Maxima 不是告诉我们2.3333...,而是分数的形式!难道Maxima 真的懂分数?不要怀疑,这就是所谓电脑代数系统(CAS) 的特长。我们可以像(%o5)的例子一样,输入个分数的四则运算试试即知。

如果坚持要用浮点数,那只要加个float 指令即可:

 (%i6) float(7/3);
 (%o6)                          2.333333333333334

为了完整,我们顺便再介绍指数,根号,阶乘表示法:

 (%i7) 2^10;
 (%o7)                                1024
 (%i8) sqrt(9);
 (%o8)                                 3
 (%i9) 5!;
 (%o9)                                120

我们可以看出,这些运算不是自然的数学符号,就是和我们平常电脑程序语言的写法。

=== 指令结尾===

在上面的例子中,我们发现,在Maxima 下指令,结束时一定要打上分号「;」,让Maxima 知道我们下的指令已结束。为什么要多这一个动作,主要是为了有时打比较长的指令可以换行之故。另一个结束方式是打入「$」的符号。不同於分号的地方是「运算结果不会显示出来」:

 (%i10) 2+3$
 
 (%i11) 2+3;
 (%o11)                                 5

有一些CAS 程序,如Matehmatica 是用分号表示不显示运算结果。不过Maxima 中分号已用上,必需用其他字符。

===离开Maxima===
离开Maxima 打入“quit();” 即可。

当然,很多人可能会觉得奇怪,为什么不是打入“quit” 就好了呢?原来像这种程序导向的语言,什么动作其实都是执行一个函数。所以我们事实上是执行一个叫「离开」的函数。这函数没有参数,所以就成了quit() 的形式。

=== 结果的引用===
我们时常会需要引用前面的结果,这时就用百分比符号“%” 。比方说:

 (%i12) 7/3;                                                                  
                                        7
 (%o12)                                 -
                                        3
 (%i13) float(%);
 (%o13)                         2.333333333333334

Maxima 也可以指定使用第几个输出的结果,不过自己定一个标签可能是最好的方式。比方说,我们可以这样用:

 (%i14) myresult:34+(65*72)/119;
                                      8726
 (%o14)                               ----
                                      119
 (%i15) float(myresult);
 (%o15)                         73.32773109243698

===重要常数===
Maxima 当然有内建e 或是π 常常用到的数,只是表示法奇怪一点。e 是%e 而π 是%pi 。

===定义变量===
Maxima定义变量的想法有点特别,在定义一个变数时,其时是给某个数字、矩阵,或想要定义的任何式子等等一个标签。让我们来看几个例子:

 (%i16) a: 37;
 (%o16)                                37
 (%i17) a;
 (%o17)                                37
 (%i18) b: 22+100*(375-128);
 (%o18)                               24722
 (%i19) a+b;
 (%o19)                               24759

===函数===
Maxima 函数的定义和使用非常直觉,我们看几个例子就知道:

 (%i20) f(x) := 3*x^2+5;
                                           2
 (%o20)                         f(x) := 3 x  + 5
 (%i21) f(2);
 (%o21)                                17
 (%i22) g(x, y) := sin(x)*cos(y);
 (%o22)                     g(x, y) := sin(x) cos(y)
 (%i23) g(2*%pi, 4);
 (%o23)                                 0

重点就是,在定义函数时要用“:=” 去定义。比较一下和变数定义的不同,想想为什么要有两种不一样的定义方式。

==进阶使用==
===列式而不运算===

我们先计算一个瑕积分,用到无穷大的部份Maxima 是以inf 表示:

 (%i1) integrate(%e^(-x^2),x,0,inf);                                          
                                    sqrt(%pi)
 (%o1)                              ---------
                                        2
还记得这在微积分是怎么积出来的吗?Maxima 居然会积!不过,今天这不是我们的重点。今天重点是,有时你不是要秀答案,只是要列出式子。我们要怎么样让Maxima 不要太自动就算出来呢?答案是加个“”号在前面,例如:

 (%i2) 'integrate(%e^(-x^2),x,0,inf);
                                  inf
                                 /         2
                                 [      - x
 (%o2)                           I    %e     dx
                                 ]
                                 /
                                  0

===kill 指令===
有时我们设定了一堆变数,函数,后来又不想再用下去,可以用kill 指令。而kill(all) 更是把我们定义过的变数,函数全部删除。看些例子就更加清楚:

 (%i3) f(x) := 3*x^2+5;
                                           2
 (%o3)                          f(x) := 3 x  + 5
 (%i4) f(x);
                                       2
 (%o4)                              3 x  + 5
 (%i5) kill(f);
 (%o5)                                done
 (%i6) f(x);
 (%o6)                                f(x)

===ev的使用===
我们可以把Maxima 的ev 指令想成一个独立的环境。有点像在写程序时的函式一样, 并不会影响到其他的运作。第一种ev 的应用是把我们设成不要执行的指令执行:

 (%i7) f: 'integrate(x^2, x);
                                     /
                                     [  2
 (%o7)                               I x  dx
                                     ]
                                     /
 (%i8) ev(f, integrate);
                                        3
                                       x
 (%o8)                                 --
                                       3

另一个很有用的使用方式是, 我们有个式子, 比方说:

 (%i9) f: a*x^2+b*2+c;
                                    2
 (%o9)                           a x  + c + 2 b

假设我们想令一个式子是a = 1, b = −2, c = −8 的情况, 我们当然可以先令各个变数是这样,们问题是这么一来, f 也永远是x2 − 2x − 8, a, b, c 这三个变数也不再是「符号」, 而是有值的。为了避免这个问题, 我们可以用ev 指令, 在下了这个指令后, 我们可以发现, 并没有变动到原来a, b, c或是f :

 (%i10) g: ev(f, a=1, b=-2, c=-8);
                                      2
 (%o10)                              x  - 12
 (%i11) a;
 (%o11)                                 a

==线性代数相关指令==
这节我们正式介绍线性代数相关,也就是矩阵相关的指令。
===矩阵及向量===
我们先来看矩阵和向量的定义方式。前面说过,在Maxima 里,所谓设定一个变数的值,只不过是给某个数字或矩阵等等一个名称。我们这里就举应用在矩阵和向量时的情况:

 (%i1) A:matrix([1,2,3],[-2,8,3],[1,4,9]);
                                  [  1   2  3 ]
                                  [           ]
 (%o1)                            [ - 2  8  3 ]
                                  [           ]
                                  [  1   4  9 ]
 (%i2) v: [2,3,5];
 (%o2)                              [2, 3, 5]

我们可以看出,要定义一个矩阵,就是把矩阵一列列的输入;定义一个向量,其实和我们用手写向量出来也差不多。不过,问题是我们在线性代数常常要把向量写成「行向量」,而非如上的「列向量」表示方式。我们可以用下面两种不同的方式达成:

 (%i3) v:transpose([2,3,4]);
                                      [ 2 ]
                                      [   ]
 (%o3)                                [ 3 ]
                                      [   ]
                                      [ 4 ]
 (%i4) v:matrix([2],[3],[5]);
                                      [ 2 ]
                                      [   ]
 (%o4)                                [ 3 ]
                                      [   ]
                                      [ 5 ]

其实向量应该是一个一列或一行的矩阵, 但是Maxima 提供了简单定义列向量的方法。这里要强调一点, 一般来说因为矩阵乘法的关系, 我们写成列向量和行向量差别很大。不过Maxima 其实不太在意这点: 它可以聪明地发现你要做的事, 并且正确得计算出来!简单的说, 一般而言, 我们不需要麻烦得定义行向量, 用列向量即可。

=== 矩阵的表示和截取===
这节我们讨论矩阵的抽象表示和取出一个矩阵行,列,甚至entry 的方法。这在很多理论和计算的尝试会用到。Maxima 是一个CAS 系统,所以我们可以完全用符号去定义一个矩阵,比方说:

 (%i5) A: matrix([a[1,1],a[1,2]],[a[2,1],a[2,2]]);
                                [ a      a     ]
                                [  1, 1   1, 2 ]
 (%o5)                          [              ]
                                [ a      a     ]
                                [  2, 1   2, 2 ]

你也可以做完全抽象的代数计算:

 (%i6) c*A;
                              [ a     c  a     c ]
                              [  1, 1     1, 2   ]
 (%o6)                        [                  ]
                              [ a     c  a     c ]
                              [  2, 1     2, 2   ]
如此一来,我们要试著导出一些定理就非常方便!

现在,我们重新把A 定义成一个实数矩阵,再看看怎么样找出A 的某一列,某一行,或某个entry。

 (%i7) A: matrix([1,2,3],[-2,8,3],[1,4,9]);
                                  [  1   2  3 ]
                                  [           ]
 (%o7)                            [ - 2  8  3 ]
                                  [           ]
                                  [  1   4  9 ]
 (%i8) row(A,1);
 (%o8)                            [  1   2  3 ]
 (%i9) col(A,2);
                                    [   2   ]
                                    [       ]
 (%o9)                              [   8   ]
                                    [       ]
                                    [   4   ]
 (%i10) A[2,3];
 (%o10)                                 3 

=== 矩阵向量之四则运算===
我们要做矩阵加法、减法、乘法非常直觉而容易。乘法用的运算元是“.”。我们假设有了前面矩阵A 和向量v 的定义,来看以下的例子:

 (%i11) A.v;
                                     [ 23 ]
                                     [    ]
 (%o11)                              [ 35 ]
                                     [    ]
                                     [ 59 ]

也可以定义非向量的矩阵试试矩阵的乘法。比方说,两个矩阵A , B 的乘积是A.B,要注意A*B 并不会得到矩阵相乘的结果!到底A*B 是什么意思,大家不妨自己试试,看可不可以找出其中的意义。

量内积的做法和你想的一样:

 (%i12) w: [2,3,5];
 (%o12)                             [2, 3, 5]
 (%i13) w.w;
 (%o13)                                38

你可能发现了一个问题,那就是我们上面内积的例子是用列向量。那行向量可以吗?可以的!Maxima会聪明的知道你想做什么,不信可以试试看。

矩阵和向量的纯量乘法是用平常的“*” 号:

 (%i14) 2*A;
                                 [  2   4   6  ]
                                 [             ]
 (%o14)                          [ - 4  16  6  ]
                                 [             ]
                                 [  2   8   18 ]

现在我们来看一下有可能会产生误会的地方。假设我们现在要算A · A ,你可能会想是A^2,结果并不正确!其实A^2 是把A 的每一个entry 都平方。正确计算A · A 要用A^^2

=== 矩阵相关函数===

我们要计算矩阵的行列式值,求转置矩阵, 矩阵的秩等等的基本运算,Maxima 当然也都有(A还是我们之前定义的矩阵):

 (%i15) transpose(A);
                                  [ 1  - 2  1 ]
                                  [           ]
 (%o15)                           [ 2   8   4 ]
                                  [           ]
                                  [ 3   3   9 ]
 (%i16) determinant(A);
 (%o16)                                54
 (%i17) rank(A);
 (%o17)                                 3

我们当然也可以手动计算行列式值。但这时需要知道矩阵第i, j 这个位置的子式(minor), 也就是A 矩阵去掉第i 列, 第j 行所成的矩阵, 这指令叫minor:

 (%i18) minor(A,1,1);
                                    [ 8  3 ]
 (%o18)                             [      ]
                                    [ 4  9 ]

矩阵的余因子(cofactor) 在Maxima 中并没有定义, 好在我们自己可以很容易定一个cofactor函数:

 (%i19) cofactor(M,i,j):=(-1)^(i+j)*determinant(minor(M,i,j));
                                     i + j
 (%o19)    cofactor(M, i, j) := (- 1)      determinant(minor(M, i, j))
 (%i20) cofactor(A,1,1);
 (%o20)                                60

我们在计算反矩阵等会用到的古典伴随矩阵(classical adjoint matrix) 也很容易算出来:

 (%i21) adjoint(A);
                               [  60   - 6  - 18 ]
                               [                 ]
 (%o21)                        [  21    6   - 9  ]
                               [                 ]
                               [ - 16  - 2   12  ]

说到反矩阵,要用Maxima 求出来也是易如反掌:

 (%i22) invert(A);
                               [  10     1     1 ]
                               [  --   - -   - - ]
                               [  9      9     3 ]
                               [                 ]
                               [  7     1      1 ]
 (%o22)                        [  --    -    - - ]
                               [  18    9      6 ]
                               [                 ]
                               [   8     1    2  ]
                               [ - --  - --   -  ]
                               [   27    27   9  ]

或是你也可以用前面的方式求反矩阵:

 (%i23) A^^(-1);
                               [  10     1     1 ]
                               [  --   - -   - - ]
                               [  9      9     3 ]
                               [                 ]
                               [  7     1      1 ]
 (%o23)                        [  --    -    - - ]
                               [  18    9      6 ]
                               [                 ]
                               [   8     1    2  ]
                               [ - --  - --   -  ]
                               [   27    27   9  ]
 
在解线性方程组常用到的梯形矩阵也是容易得很:

 (%i24) echelon(A);
                                   [ 1  2  3 ]
                                   [         ]
                                   [       3 ]
 (%o24)                            [ 0  1  - ]
                                   [       4 ]
                                   [         ]
                                   [ 0  0  1 ]

=== 使用模组===
用了Maxima 一阵子,你可能会预期它该会的都会。比方说求一个矩阵的trace,这应该够容易了吧?

事情并不是那么简单。Maxima 本身是「不会」算trace 的!当然我们可以自己写个小程序,不过先别急。我们可以使用使用适当的模组来做这件事。

所谓模组就是一段小程序,通常是增加一些指令,供你使用。你也许会觉得奇怪,那为什么Maxima 不一开始就把这些模组都加进来?那是因为如此一来太占用内存,也许很多对某些人重要的指令你永远也不用去用!

我们要算一个矩阵的trace,要使用ncharl 这个模组,这个模组提供了mattrace 指令去计算trace。

使用的方法如下,先以

 (%i25) load("ncharl");
 
读入ncharl 模组,接著就可以使用这个模组提供的指令:

 (%i26) A:matrix([1,2,3],[2,2,1],[3,3,1]);
                                   [ 1  2  3 ]
                                   [         ]
 (%o26)                            [ 2  2  1 ]
                                   [         ]
                                   [ 3  3  1 ]
 (%i27) mattrace(A);
 (%o27)                                 4

==线性代数应用实例==
    
===特征值和特征向量===
我们这里讨论线性代数很重要的特征值相关的计算。我们定义一个矩阵A , 计算特征值和特征向量时我们都以这个矩阵为主要讨论对象:

 (%i1) A: matrix([4,0,1],[2,3,2],[1,0,4]);
                                   [ 4  0  1 ]
                                   [         ]
 (%o1)                             [ 2  3  2 ]
                                   [         ]
                                   [ 1  0  4 ]

我们计算一下特征值:

 (%i2) eigenvalues(A);
 (%o2)                          [[5, 3], [1, 2]]

 么样,很方便吧...等等,特征值怎么会出来两个向量呢!?原来,真正的特征值是放在结果的第一个list 当中,也就是5 和3。那第二个list 代表什么呢?代表的就是每个特征值的几何重数, 也就是每个特征值对应的特征向量空间之维度。换言之,这是比较完整的特征值资讯!

我们也可以用eigenvectors 计算特征向量。事实上,eigenvectors 也会把特征值列出来,所以是包含前面eigenvalues 功能的指令。不过如果我们一开始就介绍eigenvectors,看到那有点复杂的结果大家可能会昏倒。现在已经会了eigenvalues,大概就没问题了:

 (%i3) eigenvectors(A);
 (%o3)        [[[5, 3], [1, 2]], [1, 2, 1], [1, 0, - 1], [0, 1, 0]]

第一部份和eigenvectors 输出一样,就是说我们有5 有和3 两个特征值,其mutiplicities 分别是1 和2。因此,对於5 应该要有一个对应的特征向量,即[1, 2, 1] ,对於3会有两个,分别是接下来的[1, 0, −1] 和[0, 1, 0] 。这些向量会生成相对应特征值的向量空间。

===手动特征值的计算===
上一节介绍Maxima 内建特征值计算,并不一定每个人都喜欢。比方说显示的方式比较特别,另外就是不是一步一步算的,心里有时也有不踏实的感觉。因此,我们这里介绍一下如何用Maxima 一步一步的把特征值求出来。

我们再用一次上一节的例子:

 (%i4) A: matrix([4,0,1],[2,3,2],[1,0,4]);
                                   [ 4  0  1 ]
                                   [         ]
 (%o4)                             [ 2  3  2 ]
                                   [         ]
                                   [ 1  0  4 ]

我们先求特征多项式,也就是A − tI 的行列式值:

 (%i5) f: charpoly(A,t);
                                               2
 (%o5)                      t + (3 - t) (4 - t)  - 3

如果想要看到比较漂亮的式子, 可以将f 展开:

 (%i6) expand(f);
                               3       2
 (%o6)                      - t  + 11 t  - 39 t + 45

我们还可以将f 做因式分解, 这样就可以清楚看到A 有几个特征值, 和各特征值的代数重数:

 (%i7) factor(f);
                                                2
 (%o7)                         - (t - 5) (t - 3)

这样我们就求得A 的特征值是5 和3 。

另一个解法是,我们可以求f = 0 的零根。做法是使用solve 指令:

 (%i8) solve(f=0, t);
 (%o8)                           [t = 5, t = 3]

当然,我们算法正确,应该是得到和前面一样的结果。

===解线性方程组===
线性代数的核心问题,就是解线性方程组。解线性方程组一样可以用上一节介绍的solve 指令来解。我们来看一个简单的例子,并且用Maxima 来解。

我们考虑下面的线性方程组:

===手动求特征向量空间的基底===
我们在前面介绍过, 使用eigenvectors 指令就可以求出特征向量空间的一组基底。我们再用一次前面的矩阵:

==Maxima 的绘图功能==
===二维绘图===
Maxima 二维绘图的指令是用plot2d。比方说,我们要画4x3 − 2x − 2 这个函数,设定x 轴范围是从-5 到5,就下这个指令:
 (%i1) plot2d([4 * x^3-2 * x-2],[x,-5,5]);

===三维绘图===
三维绘图也一样容易,只要改用plot3d 的指令即可:
 (%i2) plot3d(cos(-x^2+y^3/4),[x,-4,4], [y,-4,4]);

Geomview 是一个UNIX 的软件,Maxima 可以运用Geomview 做出非常漂亮的3D 图形。我们来看上个例子以Geomview 输出的结果。
 (%i3) plot3d(cos(-x^2+y^3/4),[x,-4,4], [y,-4,4], [plot_format,geomview]);

Geomview 不但可以画出漂亮3D 图形,更重要的是它可以弥补Maxima 的一些缺点。比方说,Maxima 本身的3D 绘图不可以同时显示两个或两个以上函数图形(2D可以),但利用Geomview,这样的绘图变成可能。

===点绘图===
有很多绘图的应用,就只需要画出点,或是用一些点来描述一些函数。这事实上比画函数还简单,但是Maxima 直到5.9.2 版才有这样的功能。详情请参考(5.9.2 之后的) 使用手册。

===多个函数的绘图===
如果要比较几个函数,要如何下指令呢?我们来看个例子就明白了:
 (%i4) plot2d([cos(x), sin(x), tan(x)], [x, -2*%pi, 2*%pi], [y,-2,2])$

这个例子会同时画出cos(x), sin(x) 和tan(x) 的图形。

===参数式绘图===
我们仅简单举一参数式绘图之例子, 详情请参考Maxima 使用手册。
 (%i5) plot2d([parametric, cos(t), sin(t), [t,-2*%pi,2*%pi], [nicks,80]]);

==Maxima 的安装==
   
在Maxima 的官方网站有不同版本的Maxima 供各平台使用:
*http://maxima.sourceforge.net/

不过,不同平台可能有一些不同的选择。我概略说明一下我建议的安装方式。不管用Windows, Mac, 或是Linux,我都推荐使用TeXmacs 这个文书处理软件当界面,因为这样可以显示最漂亮的数学符号。

==关於这份文件......==
===原繁体中文版===
Maxima 在线性代数的应用


本文件使用LaTeX2HTML , Version 2002-2-1 (1.71) 转换。

Copyright c 1993, 1994, 1995, 1996, Nikos Drakos, Computer Based Learning Unit,Universityof Leeds.

Copyright c 1997, 1998, 1999, Ross Moore, Mathematics Department, Macquarie University,Sydney.

中文转换方法详见李果正《我的CJK》中「CJK 和LaTeX2HTML 的配合」之说明。本文使用李果正taiwan.perl 套件。

===本文件===
Maxima 在线性代数的应用

由dbzhang转换,并重新排版。