个人工具

Maxima简介

来自Ubuntu中文

Dbzhang800讨论 | 贡献2007年5月19日 (六) 19:24的版本 矩阵运算

跳转至: 导航, 搜索

英文原文:http://maxima.sourceforge.net/docs/intromax/intromax.html

翻译:dbzhang

本文pdf版可在http://maxima.sf.net下载

简介

在linux中运行Maxma,键入

maxima <回车>

计算机将显示如下的欢迎词:

Distributed under the GNU Public License. See the file COPYING.
Dedicated to the memory of William Schelter.
This is a development version of Maxima. The function bug_report()
provides bug reporting information.
(%i1)


此处(%i1) 是一个“标签”。每一个输入或输出行都贴有一个标签,每一行都可以在随后的会话中通过标签被 调用。标签i代表该行是你输入的命令,标签o代表该行为机器的响应。 永远不要尝试使用形如 %i1 或 %o5 的变量名,那将会和采用该标签的行相混淆。 Maxima 对字符的大小写是敏感的。所有内建函数的函数名都是小写的 (sin, cos, save, load,等)。内建的常数采用小写形式 (%e, %pi, inf, 等)。如果你键入 SIN(x) 或者 Sin(x), Maxima 认为你指代其他的函数 而不是内建的 sin 函数。用户自定义函数和变量可以采用大写或小写的形式。 foo(XY), Foo(Xy), FOO(xy) 是不同的函数。

特殊键和符号

  1. 要结束一个Maxima会话,键入 quit();.
  2. 要终止一次计算而不退出Maxima,键入 ^C。 (这儿 ^ 代表Ctrl键,因此 ^C 意味着先按住Ctrl键,然后再按下C。)了解这一点在有些情况下对你是很重要的,比如,在一次计算需要耗费太长时间的时候。举例如下:
(%i1) sum (1/x^2, x, 1, 10000);

Maxima encountered a Lisp error:


Console interrupt.

Automatically continuing.
To reenable the Lisp debugger set *debugger-hook* to nil.
(%i2)
  1. 要告诉Maxima你已经完成了命令的输入,键入分号(;)并回车。注意到单独一个回车并不代表输入的结束。
  2. 另一个可以代替分号 (;) 的终止符是一个美元符号($),而且,它可以取消 MAXIMA 计算结果的显示,当你在进行一次有着很长结果的计算,并且你不想浪费时间显示结果的时候,这会很有用。
  3. 如果你想重复一条你已经给出的命令,比如说在 (%i5) 行,你可以在上述的行号前加两个单引号 (") 的方法来避免再次输入,比如, %i5。(注是简单地输入 %i5 并不能实现这种效果──试试吧。)
  4. 如果你想引用Maxima上一步计算的结果,你可以用它的 o 标签,也可以使用专门的百分号(%)。
  5. 标准量 e (自然对数的底数),i (-1 的平方根) 和 p (3.14159?) 分别表示成 %e,%i, 和 %pi。注意这里% 只是作为一个前缀使用, 与用 % 来查询先前计算结果的用法完全无关。
  6. 为了把一个值赋给一个变量,Maxima使用冒号 (:),而不是等号。等号被用来表示方程或等式。

算术

常见的算术操作符有:

  • + 加法
  • - 减法
  • * 标量乘法
  • / 除法
  • ^或 ** 指数
  • . 矩阵乘法
  • sqrt(x) x的平方根

Maxima输出的特点是严格的算术(有理)运算。例如:

(%i1) 1/100 + 1/101;
                          201
(%o1)                         -----
                         10100

如果计算中涉及无理数,它们将保持符号形式:

(%i2) (1 + sqrt(2))^5;

                                 5
(%o2)                    (sqrt(2) + 1)
(%i3) expand (%);
(%o3)                    29 sqrt(2) + 41

尽管如此,将计算结果用小数显示出来往往是有用的。这可以通过在你想要展开的表达式后面加上“ ,numer” 来实现。

(%i4) %, numer;
(%o4)                   82.01219330881976


注意在这里使用了% 来引用上一步的结果. 在这个版本的 Maxima 里, numer 给出 16 位的有效数字, 但最后 一位往往是不可靠的. 尽管如此, Maxima 可以通过使用 bfloat函数来提供 任意高的的精度:

(%i5) bfloat (%o3);
(%o5)                  8.201219330881976B1

显示的有效数字的位数由Maxima的变量fpprec控制,它的缺省值是16:

(%i6) fpprec;
(%o6)                          16


在这儿我们重置 fpprec 以产生100个有效数字:

(%i7) fpprec: 100;
(%o7)                          100
(%i8) %i5;
(%o8) 8.20121933088197564152489730020812442785204843859314941221#
2371240173124187540110412666123849550160561B1

注意在(%i8)中用两个单引号()来重复命令 (%i5)。Maxima可以处理非常大的数而不用近似值:

(%i9) 100!;
(%o9) 9332621544394415268169923885626670049071596826438162146859#
2963895217599993229915608941463976156518286253697920827223758251#
185210916864000000000000000000000000

代数

当我们看到 MAXIMA 在代数上的便利后,作为一个使符号运算的计算机工具,它的重要性就变得更明显了。 这里有一个分解多项式的例子:

(%i1) (x + 3*y + x^2*y)^3;

                     2             3
(%o1)                   (x  y + 3 y + x)
(%i2) expand (%);
  6  3      4  3       2  3       3      5  2       3  2
(%o2) x  y  + 9 x  y  + 27 x  y  + 27 y  + 3 x  y  + 18 x  y

                                    2      4        2      3
                            + 27 x y  + 3 x  y + 9 x  y + x


现在假定我们想把上面表达式中的 x用 5/z 替代:

(%i3) %o2, x=5/z;
      2        3                 2               3

 135 y    675 y    225 y   2250 y    125   5625 y    1875 y
(%o3) ------ + ------ + ----- + ------- + --- + ------- + ------

   z         2       2        3       3       4         4
            z       z        z       z       z         z
                                        2          3
                                  9375 y    15625 y        3

                                + ------- + -------- + 27 y
                                     5          6
                                    z          z


Maxima 内建函数ratsimp 用来进行通分操作:

(%i4) ratsimp (%);
      3  6        2  5         3           4
(%o4) (27 y  z  + 135 y  z  + (675 y  + 225 y) z

     2         3          3            2         2
+ (2250 y  + 125) z  + (5625 y  + 1875 y) z  + 9375 y  z

     3   6
+ 15625 y )/z

表达式也可以使用factor进行因式分解:

(%i5) factor (%);
                      2              3
                (3 y z  + 5 z + 25 y)
(%o5)                ----------------------

                           6
                          z

Maxima可以求解非线性代数方程系统的精确解。在这个例子中,我们使用函数solve解一个三元方程组, 三个未知数分别为 a, b, c:

(%i6) a + b*c = 1;
(%o6)                      b c + a = 1
(%i7) b - a*c = 0;
(%o7)                      b - a c = 0
(%i8) a + b = 5;
(%o8)                       b + a = 5
(%i9) solve ([%o6, %o7, %o8], [a, b, c]);

       25 sqrt(79) %i + 25      5 sqrt(79) %i + 5
(%o9) [[a = -------------------, b = -----------------,

       6 sqrt(79) %i - 34       sqrt(79) %i + 11
sqrt(79) %i + 1        25 sqrt(79) %i - 25
c = ---------------], [a = -------------------,

     10               6 sqrt(79) %i + 34
5 sqrt(79) %i - 5        sqrt(79) %i - 1
b = -----------------, c = - ---------------]]

sqrt(79) %i - 11               10

注意上面的显示由一个“列表”组成,也就是说,由两个方括号[ ... ]括住的一些表达式, 它本身还包含两个列表,每个列表包含方程组的一组解。Maxima可以轻松处理三角问题,函数trigexpand利用sum-angles公式使得每个三角函数的参数尽可能简单。

(%i10) sin(u + v) * cos(u)^3;

                     3
(%o10)                 cos (u) sin(v + u)
(%i11) trigexpand (%);
           3
(%o11)       cos (u) (cos(u) sin(v) + sin(u) cos(v))

与此相反,函数 trigreduce把一个表达式转换成几项和的形式,每一项只含有1个 sin 或 cos:

(%i12) trigreduce (%o10);
  sin(v + 4 u) + sin(v - 2 u)   3 sin(v + 2 u) + 3 sin(v)
(%o12) --------------------------- + -------------------------

               8                            8

函数 realpart 和 imagpart 返回一个复数表达式的实部和虚部:

(%i13) w: 3 + k*%i;
(%o13)                      %i k + 3
(%i14) w^2 * %e^w;

                          2   %i k + 3
(%o14)               (%i k + 3)  %e
(%i15) realpart (%);
           3       2               3
(%o15)        %e  (9 - k ) cos(k) - 6 %e  k sin(k)

微积分

Maxima可以计算导数和积分,按 Taylor 级数展开,求极限和普通的微分方程的精确解。我们从把 符号 f 定义为如下的 x 的函数开始:

(%i1) f: x^3 * %e^(k*x) * sin(w*x);
                    3   k x
(%o1)                   x  %e    sin(w x)

我们计算 f 对自变量 x 的导数:

(%i2) diff (f, x);
    3   k x               2   k x
(%o2) k x  %e    sin(w x) + 3 x  %e    sin(w x)

                                            3   k x
                                       + w x  %e    cos(w x)

现在我们求 f 对 x 的不定积分:

(%i3) integrate (f, x);
       6      3  4      5  2    7   3
(%o3) (((k w  + 3 k  w  + 3 k  w  + k ) x

  6      2  4      4  2      6   2
+ (3 w  + 3 k  w  - 3 k  w  - 3 k ) x

       4       3  2      5         4       2  2      4
+ (- 18 k w  - 12 k  w  + 6 k ) x - 6 w  + 36 k  w  - 6 k )

k x                 7      2  5      4  3    6     3
%e    sin(w x) + ((- w  - 3 k  w  - 3 k  w  - k  w) x

    5       3  3      5     2
+ (6 k w  + 12 k  w  + 6 k  w) x

  5       2  3       4              3       3      k x
+ (6 w  - 12 k  w  - 18 k  w) x - 24 k w  + 24 k  w) %e

        8      2  6      4  4      6  2    8
cos(w x))/(w  + 4 k  w  + 6 k  w  + 4 k  w  + k )


稍微改变一下语法可以得到定积分:

(%i4) integrate (1/x^2, x, 1, inf);
(%o4)                           1
(%i5) integrate (1/x, x, 0, inf);

Integral is divergent
-- an error.  Quitting.  To debug this try debugmode(true);


下面我们使用符号 f (在 %i1中已定义)和双曲正弦函数来定义符号g,并在 x = 0 点展开为泰勒 级数(高达3阶):

(%i6) g: f / sinh(k*x)^4;

                    3   k x
                   x  %e    sin(w x)
(%o6)                   -----------------
                          4
                      sinh (k x)
(%i7) taylor (g, x, 0, 3);

                   2    3   2         2    3   3
    w    w x   (w k  + w ) x    (3 w k  + w ) x
(%o7)/T/ -- + --- - -------------- - ---------------- + . . .

     4    3            4                 3
    k    k          6 k               6 k

当 x 趋向于0时 g 的极限计算如下:

(%i8) limit (g, x, 0);
                          w
(%o8)                          --
                           4
                          k

Maxima也允许导数已非计算形式表示(注意引号):

(%i9) 'diff (y, x);
                          dy
(%o9)                          --

                          dx

(%i9) 中的引号操作符表示(不求值)。没有它的话,Maxima将得到0:

(%i10) diff (y, x);
(%o10)                          0

使用引号操作符我们可以写微分方程:

(%i11) 'diff (y, x, 2) + 'diff (y, x) + y;

                      2
                     d y   dy
(%o11)                    --- + -- + y
                       2   dx
                     dx

Maxima的 ode2 函数可以求解一阶和二阶的常微分方程:

(%i12) ode2 (%o11, y, x);
        - x/2          sqrt(3) x            sqrt(3) x
(%o12) y = %e      (%k1 sin(---------) + %k2 cos(---------))

                           2                    2

矩阵运算

Maxima 可以计算行列式,以及带符号元素(也就是说,带有代数变量的元素)的矩阵的逆, 特征值和 特征向量。我们从一个元素一个元素地输入一个矩阵 m开始:

(%i1) m: entermatrix (3, 3);

Is the matrix  1. Diagonal  2. Symmetric  3. Antisymmetric  4. General
Answer 1, 2, 3 or 4 :

4;
Row 1 Column 1:
0;
Row 1 Column 2:
1;
Row 1 Column 3:
a;
Row 2 Column 1:

1;
Row 2 Column 2:
0;
Row 2 Column 3:
1;
Row 3 Column 1:
1;
Row 3 Column 2:

1;
Row 3 Column 3:
0;

Matrix entered.
                      [ 0  1  a ]
                      [         ]
(%o1)                      [ 1  0  1 ]

                      [         ]
                      [ 1  1  0 ]

下面我们求它的转置,行列式和逆矩阵:

(%i2) transpose (m);
                      [ 0  1  1 ]

                      [         ]
(%o2)                      [ 1  0  1 ]
                      [         ]
                      [ a  1  0 ]
(%i3) determinant (m);
(%o3)                         a + 1
(%i4) invert (m), detout;

                   [ - 1   a    1  ]
                   [               ]
                   [  1   - a   a  ]
                   [               ]
                   [  1    1   - 1 ]
(%o4)                   -----------------

                         a + 1

在 (%i4)中,修饰符 detout 将行列式的值保持在逆矩阵元素的外边。作为检验,我们用 矩阵 m 乘以它的逆(注意这儿用小数点表示矩阵的乘法):

(%i5) m . %o4;

                          [ - 1   a    1  ]
                          [               ]
                          [  1   - a   a  ]
            [ 0  1  a ]   [               ]

            [         ]   [  1    1   - 1 ]
(%o5)            [ 1  0  1 ] . -----------------
            [         ]         a + 1

            [ 1  1  0 ]
(%i6) expand (%);
    [   a       1                                 ]
    [ ----- + -----        0              0       ]

    [ a + 1   a + 1                               ]
    [                                             ]
    [                  a       1                  ]
(%o6)    [       0        ----- + -----        0       ]

    [                a + 1   a + 1                ]
    [                                             ]
    [                                 a       1   ]
    [       0              0        ----- + ----- ]

    [                               a + 1   a + 1 ]
(%i7) factor (%);
                      [ 1  0  0 ]
                      [         ]
(%o7)                      [ 0  1  0 ]

                      [         ]
                      [ 0  0  1 ]

要求得矩阵 m 的特征值和特征向量,我们使用函数 eigenvectors:

(%i8) eigenvectors (m);
      sqrt(4 a + 5) - 1  sqrt(4 a + 5) + 1
(%o8) [[[- -----------------, -----------------, - 1],

              2                  2
             sqrt(4 a + 5) - 1    sqrt(4 a + 5) - 1
[1, 1, 1]], [1, - -----------------, - -----------------],

                  2 a + 2              2 a + 2
sqrt(4 a + 5) + 1  sqrt(4 a + 5) + 1
[1, -----------------, -----------------], [1, - 1, 0]]

    2 a + 2            2 a + 2

在 %o8中,第一个元组(triple)给出了m的特征值,第二个给出了它们各自的重数(此处都是不重复的)。下面的三个元组给出了 m 相应的特征向量。为了从这些表达式中提取一个特征向量,我们使用 part 函数:

(%i9) part (%, 2);
           sqrt(4 a + 5) - 1    sqrt(4 a + 5) - 1
(%o9)     [1, - -----------------, - -----------------]

                2 a + 2              2 a + 2

Maxima编程

到现在为止,我们已经在交互模式下用过了 Maxima,就象个计算器一样,然而,对于那些扯进了反 复控制次序的计算,还是运行一个程序来得方便。这里我们展示了一个短且简单的程序,用来计算一 个有着两个变量 x 和 y 的函数 f 的临界点。这个程序提示用户输入函数 f,然后它计算 fx 和 fy 的偏导数,然后,它使用 Maxima 命令 solve 去获得 fx=fy=0的解。这个程序是在 Maxima 之外用一个文本编辑器写得,然后用 batch 命令装载进 Maxima。下面是程序列表:

/* ------------------------------------------------------------------------
this is file critpts.max:
as you can see, comments in maxima are like comments in C

Nelson Luis Dias, [email protected]
created 20000707
updated 20000707
------------------------------------------------------------------------- */
critpts():=(
print("program to find critical points"),
/* -------------------------------------------------------------------------

asks for a function
------------------------------------------------------------------------- */
f:read("enter f(x,y)"),
/* -------------------------------------------------------------------------
echoes it, to make sure
------------------------------------------------------------------------- */
print("f = ",f),
/* -------------------------------------------------------------------------

produces a list with the two partial derivatives of f
------------------------------------------------------------------------- */
eqs:[diff(f,x),diff(f,y)],
/* -------------------------------------------------------------------------
produces a list of unknowns

------------------------------------------------------------------------- */
unk:[x,y],
/* -------------------------------------------------------------------------
solves the system
------------------------------------------------------------------------- */
solve(eqs,unk)
)$

这个程序 (实际上是个没有参数的函数) 叫做 critpts。每一行都是一个有效的。可以由键盘上执行 的Maxima 命令。它们由逗号隔开。偏导数被贮存在一个叫做 eqs 的变量中。未知数贮存在 unk 中。这里是它运行的例子:


(%i1) batch ("critpts.max");

batching #p/home/robert/tmp/maxima-clean/maxima/critpts.max
(%i2) critpts() := (print("program to find critical points"),
f : read("enter f(x,y)"), print("f = ", f),

eqs : [diff(f, x), diff(f, y)], unk : [x, y], solve(eqs, unk))
(%i3) critpts ();
program to find critical points

enter f(x,y)
%e^(x^3 + y^2)*(x + y);
           2    3
          y  + x
f =  (y + x) %e   

(%o3) [[x = 0.4588955685487 %i + 0.35897908710869,
y = 0.49420173682751 %i - 0.12257873677837],
[x = 0.35897908710869 - 0.4588955685487 %i,

y = - 0.49420173682751 %i - 0.12257873677837],
[x = 0.41875423272348 %i - 0.69231242044203,
y = 0.4559120701117 - 0.86972626928141 %i],

[x = - 0.41875423272348 %i - 0.69231242044203,
y = 0.86972626928141 %i + 0.4559120701117]]

Maxima函数的不完全列表

看Maxima的参考手册doc/html/maxima_toc.html (在Maxima的安装目录下)。在Maxima运行 时,你可以使用 describe(函数名)。

allroots(a)

   求解多项式方程a所有的(复数)根,并把它们以数值格式(i.e. 采用16位有效数 字)列出
   来。 

append(a,b)

   将列表b追加到列表a,产生一个单一列表。 

batch(a)

   加载并运行一个文件名为a的程序。 

coeff(a,b,c)

   给出表达式a中b的c次方项的系数。 

concat(a,b)

   生成符号ab,比如concat(y,44)的结果为y44。 

cons(a,b)

   将a加入列表b的头部。 

demoivre(a)

   将表达式a中的复指数项变换为等价的三角形式。 

denom(a)

   给出表达式a的分母。 

depends(a,b)

   声明a是自变量b的函数。这在书写微分方程的时候很有用。 

desolve(a,b)

   使用拉普拉斯变换求解线性常微分方程a的未知量b。 

determinant(a)

   给出方阵a的行列式。 

diff(a,b1,c1,b2,c2,...,bn,cn)

   给出a对变量bi的ci阶偏导数。diff(a,b,1)可简写为diff(a,b)。’diff(...) 代表不经过计算(unevaluated)的求导,这在书写微分方程的时候很有用。 

eigenvalues(a)

   返回两个列表,第一个列表是矩阵a的本征值,第二个是本征值对应的重复次 数。 

eigenvectors(a)

   包含eigenvalues所有功能,并且计算矩阵a的本征向量。 

entermatrix(a,b)

   引导用户一个一个元素地输入一个a × b 的矩阵。 

ev(a,b1,b2,...,bn)

   在bi 的条件下计算表达式a的值。bi可以是方程、方程构成的列表(比如solve 返 回 的 结 果 ) 或 者 赋 值 语 句 , 在 这 种 情况 下 ,ev将bi“插 入”到 表 达 式a中 。bi还 可 以 是关键词numer(它让结果以数值格式显示),detout(它使任一矩阵的逆矩阵把行列
   式 的值作为系数放在矩阵外),或者diff(它要求所有的微分都必须计算,即’diff被diff替 代 ) 。 对manual command(即 , 不 在 用 户 自 定 义 函 数 内),ev可 以 省 略 ,
   于 是 可 简 写 为a,b1,b2,...,bn. 

expand(a)

   展开表达式a。 

exponentialize(a)

   将a中的所有三角函数转换为它们对应的复指数形式。 

factor(a)

   对表达式a进行因式分解。 

freeof(a,b)

   如果a不是表达式b的一部分,返回true。 

grind(a)

   Displays a variable or function a in a compact format.When used with writefile and an editor outside of Maxima, it offers a scheme for producing batch files which include Maxima- generated expressions. 

ident(a)

   返回一个a × a的单位矩阵。 

imagpart(a)

   返回a的复数部分。 

integrate(a,b)

   计算a对变量b的不定积分。 

integrate(a,b,c,d)

   计算a在区间b ∈ [c, d]上的定积分。积分限c,d可以分别取minf(负无穷大), inf(正
   无穷大)。 

invert(a)

   计算方阵a的逆矩阵。 

kill(a)

   从当前的Maxima环境中移除变量a以及它的属性。 

limit(a,b,c)

   计算当b趋近于c时a的极限。与积分函数integrate一样,c可以取inf或minf。 

lhs(a)

   给出方程a的等号左边部分。 

loadfile(a)

   从 磁 盘 的 当 前 目 录 中 加 载 文 件 名 为a的 文 件 。 该文件必须具有正
   确 的 格 式 (i.e. 由save命令创建)。 

makelist(a,b,c,d)

   创 建 一 个a( 假 定a以b为 自 变 量 ) 的 列 表 , 从b=c到b=d依 次 将a追 加 到 列 表。 

map(a,b)

   Maps the function a onto the subexpressions of b. 

matrix(a1,a2,...,an)

   创建一个以ai为行向量的矩阵a,每一个行向量是一个包含m个元素的列 表[b1, b2, ..., bm]。 

num(a)

   给出表达式a的分子。 

ode2(a,b,c)

   求解一阶或二阶常微分方程a,其中b是c的函数。 

part(a,b1,b2,...,bn)

   首先取表达式的第b1部分,然后再在该部分中再取b2部分,依次... 

playback(a)

   Displays the last a (an integer) labels and their associated expressions. If a is omitted, all lines are played back. See the Manual for other options. 

ratsimp(a)

   化简并以两个多项式的商的形式显示。 

realpart(a)

   返回a的实部。 

rhs(a)

   给出方程a的等号右边部分。 

save(a,b1,b2,..., bn)

   在磁盘的当前目录下创建包含变量、函数或矩阵bi的文件a。该文件可以 在以后的会话中用loadfile命令重新载入。如果b1取all的话,每一个符号(包括标签)都可 以得以保存。 

solve(a,b)

   求解关于未知数b的代数方程a,将返回一个根的列表。简单起见,如果方程a是c=0的 形式(即
   方程右侧为0),a可以用表达式c替代。 

string(a)

   将表达式a转换为Maxima的线性表示(类似a的Fortran表示)就像是它被输入并放入 一个缓
   冲区以用来进行可能的编辑。这样string以后的表达式不能用于后续的计算。 

stringout(a,b1,b2,...,bn)

   在当前缺省磁盘目录下创建关于变量bi(比如labels)的文件a。该文 件采用文本格式并且不
   能被Maxima再次读入。尽管如此,这种字符串化的输出只需经过稍 许修改就可用于Fortran,Basic或C程序。 

subst(a,b,c)

   将表达式c中的b用a来替换。 

taylor(a,b,c,d)

   将表达式a在b=c处展开为泰勒级数,展开的级次不超过(b − c)d 。Maxima也支 持超过一个自变量的泰勒展开,详细资料查看手册。 

transpose(a)

   给出矩阵a的转置。 

trigexpand(a)

   这是一个三角化简的函数,它采用sum-of-angle公式使得每一个sin和cos函数的变 量尽可
   能简单。例如:trigexpand(sin(x+y))的结果为cos(x) sin(y) + sin(x) cos(y). 

trigreduce(a)

   这是一个三角化简的函数,它采用三角恒等式将乘积或幂函数变换为sin或cos的 和的形式,每
   一项只含有一个sin或cos。例如:trigreduce(sin(x)^2)的结果为(1 - cos(2x))/2. 

trigsimp(a)

   这也是一个三角化简的函数,它将表达式中的tan,sec等函数变换为cos和sin函数。 变换过程中它也使用恒等式sin()2 + cos()2 = 1。