星期六, 一月 01, 2005

maxima 笔记

Maxima

GNU Maxima 是一种用LISP编 写的自由计算机代数系统, 用于公式推导和符号计算。它由麻省理工学院在 美国能源部的支持下于60年代末创造的 Macsyma 演变而来。 从1982年开始, Bill Schelter 教授默默无闻地维护 Maxima 直到2001年他去世。

1998年,Schelter获得能源部许可,以GNU GPL许可证发行Maxima。现在的 Maxima 由独立的用户和开发人员维护。

上面是在wikipedia上copy下来的简介.

本页导航: tips, 编程, 扩展, 命令, 变量, 问题, 参考资料

我 粗略地估算了一下, 大概有1600多个命令, 估算方法是按个首字母, 然后按tab键自动补齐, 估计每个首字母的命令个数, 然后加起来. ft, 是1806个, 直接按tab键它就会问我是否显示所有1806个可能选项.

tips

1. 在emacs中运行maxima(这个既方便又好看):
M-x imaxima (但是现在运行有问题, 等一会儿, 然后按C-g, C-x b *maxima*, 才行, 而且是字符显示而不是latex公式, 所以还不如用命令行; 另外现在输入输出提示符不是'C1,D1'而是'%i1,%o1', 不知道是否因为升级了, 检查版本号又好象没有升级, 是11.6号, 这个可能对emacs出现不兼容)
前提是装好了emacs的Maxima mode包, 并在.emacs文件中加入:
(autoload 'maxima "maxima" "Maxima interaction" t)
一些快捷键: M-tab: 自动补齐; 鼠标中键可以将前面的输入送到当前, 以便编辑;


2. 常用命令:
退出: quit(); (注意要有括号和分号) 不能用ctrl-c(新版已经可以用Ctrl-c来取消了), 会进入调试模式, 退出用:q

3. 一般用法:
3.1 .max文件的用法
如果foo.max文件(不一定要这个后缀名)是一个一串maxima执行语句, 可在命令行用 "maxima <> foo.out" 将输出结果直接写入foo.out文件中.

3.2 save and load:
比较郁闷的是maxima中没有保存的功能, 保存下来的东西是文本形式的, 再次打开就不认
识. 可以用batch代替, 作法是:先写好一个maxima能执行的文件"foo", 然后在maxima中batch("foo"); 就会批量执行, 如果要修改什么的话, 去修改原来的文件, 再在这边batch.
例:
文件foo:
c : 2$
c^2;
maxima中:batch("路径/foo"); 会输出:
2
4

4. 怎么写下标?
Example: n\\_1即可,就是要两个斜线。 参见14.1"与Latex结合"

5. Maximia直接将单词转化为希腊字母,如输入pi;输出则为希腊字母的“pi”。但是后面跟了东西就不变 了,如 gamma_3。
但是eta和ETA都是小写的,gamma和GAMMA都是大写的,怎么区分大小写呢?

6. 关于保存然后打开:

6.1
(C1) x:'y;
(D1) y
(C2) foo(x):=1+x;
(D2) foo(x) := 1 + x
(C3) save("myworld.lisp",all);
(D3) myworld.lisp
You restore a saved session by loading the corresponding file in the
usual way
(C1) load("myworld.lisp");
(D3) myworld.lisp
(C4) foo(x);
(D4) y + 1
(C5) fundef(foo);
(D5) foo(x) := 1 + x
但是我这样做就不成功, 初步原因是保存后的文件最开始是三个分号";;;", 在lisp里是说明, 但是在maxima里是结束的标志, maxima看到这个后就报错, 因为分号前没有任何动作, 然后就不往下继续了. 不知道怎么解决.

6.2 用writefile("foo"); 和 closefile(); 配对, 可以把这二者之间的内容原封不懂地保存到foo文件中, 用来以后查看. 注意, 只能看.

7. 启动的maxima总是自己开一个页面,使得顶上的tab栏中其它的buffer都看不到了,解决办法:
1. M-x maxima/imaxima 2. quit(); 3. M-x imaxima

8. 加撇号表示不演算
如: diff(x^2,x); 的结果是2x, 'diff(x^2,x) 是 d x^2 / d x

9. 常数的表示法:
前面加%, 如%pi, %e. 好象也就只有这两个常数.

10. 自带例子:/usr/share/maxima/5.9.1/demo/
用法: demo("/usr/share/maxima/5.9.1/demo/demo.dem");
demo和batch功能相同, 区别是demo是每执行一步要敲一下回车, batch一次性执行到结束或错误处.

11. 获取帮助的方式:
(1) help(); (2) describe(); example(); (3) demo(); (4) info maxima 内容与(2)相同 (5) maximabook (6)maillist 几乎就这么多了. info里边有个function and variable index, 如果想详细了解某个命令的的说明, 可以参考, 比较实用.

12. 任意精度
关键词: bfloat, FPPREC, block
例(copy from mailist):
(C1) float(%pi);
(D1) 3.141592653589793
(C2) bfloat(%pi);
(D2) 3.141592653589793B0
(C3) float(%pi*10^5);
(D3) 314159.2653589793
(C4) bfloat(%pi*10^5);
(D4) 3.141592653589793B5
(C5) FPPREC:100;
(D5) 100
(C6) float(%pi);
(D6) 3.141592653589793
(C7) bfloat(%pi);
(D7)
3.14159265358979323846264338327950288419716939937510582097494459230781640#

6286208998628034825342117068B0
(C8) FPPREC:5;
(D8) 5
(C9) float(%pi);
(D9) 3.1416
(C10) bfloat(%pi);
(D10) 3.1416B0
(C11) FPPREC:16;
(D11) 16
(C12) float(%pi);
(D12) 3.141592653589793
(C13) bfloat(%pi);
(D13) 3.141592653589793B0
(C14) FPPREC:40;
(D14) 40
(C15) float(%pi);
(D15) 3.141592653589793
(C16) bfloat(%pi);
(D16) 3.141592653589793238462643383279502884197B0
(C17) block([FPPREC:100],bfloat(%pi));
(D17)
3.1415926535897932384626433832795028841971693993751058209749445923078164#

06286208998628034825342117068B0
(C19) bfloat(%pi);
(D19) 3.141592653589793238462643383279502884197B0

小结: float和bfloat的精度默认都是16, 如果FPPREC>16, 只影响bfloat(float的最大精度为16, 可以小但不能大), 而如果FPPREC<16, 则float和bfloat都会受影响; 但是FPPRINTPREC的值无论多少都只对bfloat有效.

但是用tex来显示比较小的指数的时候并不听话:
(%i2360) FPPREC:6;

(%o2360) 6
(%i2361) tex(float(1/3));

$$0.3333$$
(%o2361) FALSE
(%i2362) tex(float(1/3*10^-24));

$$3.3333333333333335 \times 10^{-25}$$
(%o2362) FALSE
(%i2363) tex(float(1/3*10^-2));

$$0.0033$$
(%o2363) FALSE

13. 关于注释
方法与c语言相同, 是 "/*...*/". 在emacs里只要/* 后面的就都变颜色了, 但是语法上还是要两个对应, 否则maxima会报错.

14. 与latex结合

14.1 如何输入latex中的A_{2,1/2}?
A_\{2\,1\/2\}
注意那个逗号,除号前面都要加斜杠(花我好长时间找帮助啊, 最后还是试出来的)=_=.
例2: A_{*,-1}要表示成A_\{\*\,\-1\}
如果是特殊字符, 在latex中是一斜杠, 在maxima中要用两斜杠, 如
\alpha 要写为 \\alpha

14.2 更方便地输入一个比较复杂的符号
(%i2) load("mactex-utilities")$
(%i3) texput(st,"\\Sigma_t");
(%o3) \Sigma_t
(%i4) tex(st^2);
$$\Sigma_t^2$$

14.3 对于下标只有一个数字的情况
(%i9) tex(k_1);

$$k__1$$
(%o9) FALSE
(%i10) tex(k_\{1\});

$$k_{1}$$
(%o10) FALSE

不知道是不是一个bug.

14.4 tex()是一个将maxima的表达式变成latex表示的工具
一个需要小注意一下的:
system("echo 'The flux density for the case $\\nu < \\nu_c < \\nu_m$
$F_\{\\nu\,j\}$: ' >> jet.tex")$
可以!
system("echo 'The flux density for the case $\\nu < \\nu_c < \\nu_m$ $F_\{\\nu\,j\}$: '
>> jet.tex")$
不行!
想一想, 在命令行中下面的也是不行的, 不准许还没写完, 也没有续行符就直接按enter的. (与tex关系不大)

14.5
EXPTDISPFLAG:FALSE;
使得tex(x^(-2))的输出是$x^{-2}$而不是$1\over x^2$.


15. 几个图形界面
xmaxima, wxmaxima, symaxx
xmaxima: 独立于maxima, 但界面并不比maxima好看, 一样不能顺眼地显示公式, 是用文本凑的.
wxmaxima: 是maxima的X界面, 以maxima为核心, 显示的公式好看点.
symaxx: 也是maxima的X界面, 多年没有更新了, 不过好象也还不错.
其中wxmaxima可能最好, 不过还是没有结合emacs和latex的imaxima好.

16. 怎样是变量以1/2为指数而不是开根号出现?
有一个专门的变量 SQRTDISPFLAG, 设成 false即可, 另外
DISPLAY_FORMAT_INTERNAL : true;
这个变量还有如下功能:
(%i108) a-b;

(%o108) a + (- 1) b
(%i109) a/b;

- 1
(%o109) a b
(%i110) sqrt(a);

1/2
(%o110) a

maxima中的flag是个很有用的东东, 设置flag的值可以用来控制各种各样的设置, 因为有各种flag, 象:number, exponential, expansion, fatoring. flag就是variable了. 主要控制显示的, 关键词: display.

17. emaxima
在latex中写mamixa代码并执行, 然后都显示在latex文件中. 参考maxima book 3.2.4节. 这里还有一个介绍EMaximaIntro.ps.
例子:
\begin{document}
\usepackage{emaxima}
\begin{document}
\beginmaxima
x^2;
\endmaxima
\end{documnent}
注意要在emacs环境中做, 这样才能用M-x emacs-update-all 来执行其中的maxima代码.
maxima book本身就是用emaxima写的, 所以是最好的example. 现在还有一个例子是自己写的orphan.tex.

\beginmaxima 和 \beginmaixmasession 的区别在于前者把编译前\maximaoutput之前的内容也显示出来了, 而后者没有.




18. 物理常数
load(physconst)$ 就可以直接引用物理常数了.

19. forget和kill
assume(x>1); forget(x);
assume(x>y); forget(x>y);
y:3; kill(y);
y(x):**; kill(y);

小注意: 用kill(all); 来把前面的变量都kill掉, 如果想重新做什么东东的话. 一个好的习惯是在用某个量之前先初始化, 要么赋值(当然可以是表达式), 要么kill它. 也可以kill(a,b);这样kill多个变量的.

20. 随时加个圆括号()好象是没有问题的, 就象latex里没事加个大括号{}样.


21. 赋值, 适合于解出方程后赋值
lhs(x=3)::rhs(x=3); 就相当于x:3;了, 但是不能用lhs(x=3):rhs(x=3);
例如:
(%i1) aa:solve(x-3=0);

(%o1) [x = 3]
(%i2) x;

(%o2) x
(%i3) lhs(aa[1])::rhs(aa[1]);

(%o3) 3
(%i4) x;

(%o4) 3

22. 关于大小写
虽然命令不分大小写, 但是变量是分的, 例如
(%i193) I:2;

(%o193) 2
(%i194) i;

(%o194) i

Maxima开发者好象在做关于大小写的事, 也许是做得更符合习惯吧.

23. noun和verb
maxima中分名词形式和动词形式, noun就是原样照抄, verb就是执行. 不用太在意, 只不过是里边的一个说法而已, 便于表达. 如'sin()就是noun, sin(x)就是verb, 而'diff用的比较普遍.

24. 对前文的引用(quote): %,%%,%i*,%o*,%th(i)
%: 表示前一个输出的结果;
%%: 表示前前一个输出的结果;
%th(i): 表示前i个输出的结果, 注意是相对的, 例%th(1)与%相同;
%i*: 第*个输入, 这是系统自带的仅有的一个引用输入结果的;
%o*: 第*个输出, 和%i是两个绝对量;
另外, 自定义一个, 在mail list上问到的, Viktor回答的,
thi(i):=block([tmp],tmp:concat('%i,linenum-i),ev(tmp,nouns));
然后thi(i)就表示前i个输入的结果了, 当然也可以定义为%thi(), 不过以%为前缀的一般是系统自有的, 最好不要强人家生意:)
这个定义的解释(也是他解释的): concat('%i,linenum-i)将%i和(linenum-i)连接起来, 其中linenum是当前的行数, linenum-i就是所指目标的行数, '%i表示%i不执行(即名词形式), 将这个连接起来的量暂时存到tmp中, ev(tmp,nouns), 将tmp动词化(nouns的作用), 就是所想要的那一行. 实际上就是借用了%i*而已.
发挥一下thi(i):=block([tmp],tmp:concat('%i,i),ev(tmp,nouns)); 就是引用标号为i的那一行的输入. 简化一下就是
thi(i):=(tmp:concat('%i,i),ev(tmp,nouns)); 好象block可有可无, [tmp]也不知道有什么用--申明为本地变量, 而不是全局变量的.

25. 画图
最简单的例子:
plot2d(sin(x),[x,-%pi,%pi]);

一个比较完整的例子:
beta_1: 1.0$
c : 3.D10$
R : 1.D18$
theta_p : 0.50$
gamma: 10$
beta: sqrt(1-1.0/gamma^2)$
/*ksi: c*t/R$*/
ksi0: 0.00001$
assume(ksi>0+ksi0, ksi f_nu(ksi):= (1.0-beta*(1.0-(ksi-ksi0)))^(-(2.0+beta_1)) *sqrt(1.0-(ksi-ksi0)/2.0) /(sqrt(1.0-(1.0-(ksi-ksi0))^2)*sqrt(1.0-(ksi-ksi0)/2.0/sin(theta_p)^2));
plot2d(f_nu,[ksi,1.D-6+ksi0,ksi0+2.0*sin(theta_p)^2-1.D-8], [GNUPLOT_PREAMBLE, "set logscale xy 10"], [GNUPLOT_TERM, PS], [gnuplot_out_file, "tail.eps"]);
_____________
其中[GNUPLOT_PREAMBLE, "set logscale xy 10"], 可改写为[GNUPLOT_PREAMBLE, "set logscale x 10; set logscale y 10"], 注意是分号.

26. 浮点数
设置浮点数的精度:FPPREC, 默认值:16;
设置浮点数的显示长度:FPPRINTPREC, 默认值:0, 即浮点数的事迹精度, 仅对任意精度的bfloat有效;

27. 从表达式中抽值
(%i2) FIRST (3*x*y);
(%o2) 3
还有second...tenth, 比如这里就可以把一个单项式的系数抽出来.

28. SQRTDISPFLAG: FALSE $
使得开根号也以指数形式表示.

29: 局域变量
(%i19) b:1;
(%o19) 1
(%i20) block([b],b:2);
(%o20) 2
(%i21) b;
(%o21) 1
可见block([b]是把b当作局部变量的, 注意block(b,...)不可以, ([b],b:2);也不行.

30: RESET ();
将变量重新设定. 可以通过这个命令看有那些环境变量.

31: 一个编程的好习惯
每一行的内容都可以直接在行首和行尾用/*...*/说明掉. 也就是每一行都是独立的, 这样, 有一些括号, 特别是block的回括号, 就要单独成行(当然, 这个会括号不能随便说明掉).


用maxima 编程
方法一: 写好一个文件, 这个文件是一些命令的堆积, 然后用batch命令执行, 例:
print(The first program);
diff(x^2,x);
然后保存为first.max, 其实扩展名是什么无所谓, 在命令行中输入maxima,
(%i1)batch("first.max");
这样即可.

方法二: 用个括号括起来, 里边的命令之间用逗号隔开, 然后batch它. 好象前面加不加block无所谓.

需要用到的命令: block, ev, map
map: 定义一个过程, 并作用到右面的表达式上,
例:
(%i3) MAP(F,x-y);
(%o3) F(x)-F(y)
(%i4) MAP("=",[A,B],[-0.5,3]);
(%o4) [A = - 0.5, B = 3]

maxima 包
(在/usr/share/maxima/5.9.1/share中找)
用法: load("***.mac"); 也许要加上绝对路径名.
说明: 一般会有一个跟mac文件同名的usg文件, 是用法的说明, 有例子, 也有可能包文件是个lisp文件.

physics
dimension.mac 用来做量纲分析的, 还有一个pdf的说明文件;
physconst.mac 存了很多物理常数在里边;
units.mac 用来做单位转换的;
dimen.mac I don't know.






命 令的用法:
block
block([x,y],arg1,arg2...);
用来将一堆命令堆在一起, 顺序执行, 中间用逗号隔开, 返回最后一个值, 其中方括号的内容表示是本地变量。
和(。。。)的区别就在于block中[x]可以被认出是本地量, 如果不加block这几个字, 也能顺序执行。
多数是用block来写一个子程序或函数, 如
f(x):block([x],x+3);
ev
在ev环境中一步一步地往下做, 输出的是最后的结果, 例
(%i22) ev(x:1/3,float);
(%o22) 0.33333333333333

radcan
化简得比较彻底一点. 基本上就是把小数化成分数, 把指数, 对数中的东西化简.
小心:
(%i832) RADCAN(1.E-38);
RAT replaced 0.0 by 0//1 = 0.0
(%o832) 0
会把比较小的量变成0.

变 量的含义:
numer
将含有特殊数字的字符转换成数字, 如果numer为false(默认), 就不转换了.
(%i113) NUMER;

(%o113) FALSE
(%i114) %pi;

(%o114) %PI
(%i115) NUMER:true;

(%o115) TRUE
(%i116) %pi;

(%o116) 3.141592653589793

问 题:
& nbsp;lisp的文件怎么加载? 肯定不能batch, 因为lisp的说明是用; 表示的, 而;在maxima里是表示语句结束, 我试过, 失败了.
貌似是用load("***.lisp");
解出方程后, 变量怎么在其他地方用? 如结果是[x=4], 但用x的时候, 它还是它自己:( --已解决, 用"::".
在maxima环境中tex然后tab, 看到了tex, texput, texinit, texend, 但是关于后面三个的说明在任何地方都没有找到, 到哪儿去找呢? 这是一个更一般的问题, 就是在info里还找不到的东东, 到哪儿去找?

参考资料:
作 为参考, maxima info是不可少的, 特点是全, 在emacs里边看也很方便(C-h i); 而maximabook是很实用的, 特点是实用*.*, 一直在更新; mecref.pdf后面的索引比较有查阅价值, 既全面, 全文安排的又有逻辑解构; 下面还有不少其他的参考资料, 作为补充.

这里有一份参考DOE-Maxima Reference Manual, 一个叫Mike的人针对maxima5.9写的参考手册, 是比较全面而又有逻辑体系的.
而一个叫refman16.pdf的文件(本地改名为Macsyma-refman16-1996.pdf, 下载自Richard Fateman 的主页:
http://www.cs.berkeley.edu/~fateman/macsyma/docs), 是Macsyma的官方文档, 内容很丰富, 包括了programming, 目前其他参考文献都没有的.

这 里是一份macsyma的介绍, file:///usr/share/maxima/5.9.1/doc/html/intromax.html 叫introduction to MACSYMA, 是2000年的.
而这个才是maxima的介绍, file:///usr/share/maxima/5.9.1/doc/html/maxima_toc.html, 比上面intromax.html要详细很多, 也很新, 2004年10月更新的, 但好象就是maxima info的html版.
Minimal Maxima, 目前的维护者Robert Dodier写的20页的简短而又全面的介绍.
这里是一个" Micro introduction": http://www.math.harvard.edu/computing/maxima/index.html. 一页, 所以弄到本地了, 有一些例子.

有一个叫MaximavsMupad的pdf文档, 介绍Maxima和Mupad的异同, 不过里边举了很多例子, 可以用来学习Maxima的.

TeXmacs与maxima结合,显示的效果就要好不少, 这里(本地)有一份别人的介绍, 原网址是: http://www.inp.nsk.su/~grozin/max-tut/

maxima的mailing list: http://dir.gmane.org/gmane.comp.mathematics.maxima.general, 可以搜索的. 注意, http://gmane.org/,搜集了各种free software的mailing list.

一些个人写的介绍也值得看的, 如王垠的: http://learn.tsinghua.edu.cn/homepage/2001315450/maxima.html, 好象现在还只有他的, 呵呵.

这上面有一些用maxima写成的程序, http://www.mines.edu/fs_home/whereman/software.html, 不过主要是mathematica的程序; 还有 http://www.risc.uni-linz.ac.at/research/combinat/risc/software/Zeilberger/, 是解关于超几何函数的, 要发e-mail申请密码;
这里有一些maxima的功能扩展, 三个画图的, 一个用龙格库塔法数值解微分方程组的, 是葡萄牙人写的 http://fisica.fe.up.pt/maxima/index.html;

别人的笔记.
http://www.vttoth.com/maxima.html maxima一个开发者的网站, 他关心的主要是和相对论相关的矩阵方面.

3 条评论:

Y. C. Zou 说...

新的版本的maxima里貌似GNUPLOT_PREAMBLE 这些东西都要小写

Y. C. Zou 说...

如果要画图的时候写些特殊字符什么的, 可参考这里, 虽然是gnuplot的
http://gnuplot.sourceforge.net/demo/enhancedtext.html

一个复杂点的例子是:

c : 3.0D10$
q_e : 4.803242D-10$
m_e : 9.109534D-28$
m_p : 1.6726485D-24$
pi : 3.14159265358979D0$
sigma['T] : 6.652448D-25$

assume(Y>0);

f_nu_gamma_p:2.7e-25; /*Jy->cgs*/
f_nu_opt:2.87e-22;
nu_opt:5.45e14;
nu_gamma_p:1.6e20;
D_l:1.9e28;

z:0.937;

g_e:1000$
R_BB(Gamma):=(f_nu_opt*Gamma*D_l^2/(2*pi*nu_opt^2*g_e*m_e))^(1/2)$
dt:1$
R_dt(Gamma):=2*Gamma^2*c*dt$
n:1$
A:3e33$
E_tot:1e55$
R_dec_ISM(Gamma):=(3*E_tot/(4*pi*n*Gamma^2*m_p*c^2))^(1/3)$
R_dec_wind(Gamma):=E_tot/(4*pi*A*Gamma^2*m_p*c^2)$


plot2d([R_BB(Gamma),R_dt(Gamma),R_dec_ISM(Gamma),R_dec_wind(Gamma)],[Gamma,100,10000],[logx],[logy], [style,[lines,2],[lines,1],[points,2],[points,2]], [xlabel, "{/Symbol G}"], [ylabel, "R"], [legend,"{/=16F_{sa}=F_{opt}}"," {/=16 {/Symbol d} t=R/{/Symbol G}^2 c}","{/=16 R=R_{dec} (ISM)}","{/=16 R=R_{dec} (wind)}"], [gnuplot_term, ps], [gnuplot_out_file, "R-Gamma-ISM.eps"]
);

Y. C. Zou 说...

难怪上面画个图那么困难的, 原来用错命令了! 要和GNUPLOT接合, 应该是draw, 而不是plot!
这里有详细的介绍:
http://www.telefonica.net/web2/biomates/maxima/gpdraw/