计算机代数系统有很多，功能也侧重于不同的方面。最有名的当属Wolfram Research出品的Mathematica，其他比较有名的有Maple, OpenAxiom 等，Maxima也是其中比较有名的。

Maxima下载：

https://sourceforge.net/projects/maxima/files/

Android 下也有port ：

https://github.com/YasuakiHonda/Maxima-on-Android

Maxima 是用lisp实现的，性能不太好，大规模计算时要小心优化。下面给出 Maxima 几个简单的例子：

1. 计算Dirac γ 矩阵的迹

Tr.mac

Tr(u):=
block([ele_num:length(u), result:0, i, sublist1:delete(u[1], u), sublist2],if oddp(ele_num) then return(0)else (if ele_num = 2 then return(MT(u))else (for i:2 thru ele_num do (sublist2:delete(u[i], sublist1),result:result+(-1)^(i-1)*MT(u[1], u[i])*Tr(sublist2)),return(result)))
);

上面是定义如何计算trace的函数，要调用它，可以在同一个目录下运行maxima，然后load("Tr.mac")，就可以直接调用了。示例代码如下：

praise@aliyun:~/Template/Maxima$ maximaMaxima 5.24.0 http://maxima.sourceforge.net
using Lisp GNU Common Lisp (GCL) GCL 2.6.7 (a.k.a. GCL)
Distributed under the GNU Public License. See the file COPYING.
Dedicated to the memory of William Schelter.
The function bug_report() provides bug reporting information.
(%i1) load("Tr.mac");
(%o1)                               Tr.mac
(%i2) Tr([a,b,c,d]);
(%o2)  - MT(a, b) MT([c, d]) + MT(a, c) MT([b, d]) - MT(a, d) MT([b, c])
(%i3)

2. 解线性方程

Maxima也有图形界面wxMaxima，下面是在wxMaxima中求解线性方程的代码，将其复制并保存为 .wxm 文件，即可在 wxMaxima 中打开。

/* [wxMaxima: input   start ] */
/* -------------------------------------Solve the linear equation systemA.x = b
--------------------------------------- */
/* [wxMaxima: input   end   ] */
/* [wxMaxima: input   start ] */
/* --------- define A ------------ */A:matrix([1,(T1-delta),(T1-delta)^2,(T1-delta)^3],[1,T1,T1^2,T1^3],[0,1,2*(T1-delta),3*(T1-delta)^2],[0,1,2*T1,3*T1^2]);
/* [wxMaxima: input   end   ] */
/* [wxMaxima: input   start ] */
/* ------------- get inverse A -------------- */Ainv:invert(A);
/* [wxMaxima: input   end   ] */
/* [wxMaxima: input   start ] */
/* ------------- define b -------------- */b:matrix([t0/T1*P+D],[D],[P/T1],[P/T2]);
/* [wxMaxima: input   end   ] */
/* [wxMaxima: input   start ] */
/* ----------------- solve -------------- */result:Ainv.b;
/* [wxMaxima: input   end   ] */
/* [wxMaxima: input   start ] */
/* ------------------ simplify ------------------ */eq:T1-delta=t0;
x:scsimp(result,eq);
/* [wxMaxima: input   end   ] */
/* [wxMaxima: input   start ] */
/* ------------------- check ----------------- */subst([T1=5,T2=10,t0=4.9,delta=0.1,P=0.2], A.x);
subst([T1=5,T2=10,t0=4.9,delta=0.1,P=0.2], b);
/* [wxMaxima: input   end   ] */

可以看到，wxMaxima中的注释采用和c/c++相同的风格。

3. 自然边界条件下的三次样条拟合

/* [wxMaxima: input   start ] */
/* -------------------------------------------------------------this is a script for cubic spline with nature boundary condition------------------------------------------------------------ */
/* [wxMaxima: input   end   ] */
/* [wxMaxima: input   start ] */
/* --------------- Definition ----------------- */point_num : 8;array(knots,point_num,2);
array(kcoe,point_num);
array(aij,point_num,2);
array(bi,point_num);
/* [wxMaxima: input   end   ] */
/* [wxMaxima: input   start ] */
/* --------------- Test data set --------------- */for i:0 thru point_num step 1
doblock(knots[i,0]: i,knots[i,1]: -8+i*16/point_num,knots[i,2]: sin(knots[i,1]/4*%pi));
/* [wxMaxima: input   end   ] */
/* [wxMaxima: input   start ] */
/* =============== The first variable spline ================ */
/* [wxMaxima: input   end   ] */
/* [wxMaxima: input   start ] */
/* -------------- Initiate the linear system --------------- */for i:0 thru point_num step 1
doif i = 0 thenblock(aij[i,0] : 2,aij[i,1] : 1,aij[i,2] : 0,bi[i] : 3*(knots[1,1]-knots[0,1])/(knots[1,0]-knots[0,0]))elseif i = point_num thenblock(n : point_num,aij[i,0] : 1,aij[i,1] : 2,aij[i,2] : 0,bi[i] : 3*(knots[n,1]-knots[n-1,1])/(knots[n,0]-knots[n-1,0]))elseblock(aij[i,0] : 1/(knots[i,0]-knots[i-1,0]),aij[i,1] : 2*(1/(knots[i,0]-knots[i-1,0])+1/(knots[i+1,0]-knots[i,0])),aij[i,2] : 1/(knots[i+1,0]-knots[i,0]),bi[i] : 3*((knots[i,1]-knots[i-1,1])/(knots[i,0]-knots[i-1,0])^2+(knots[i+1,1]-knots[i,1])/(knots[i+1,0]-knots[i,0])^2));
/* [wxMaxima: input   end   ] */
/* [wxMaxima: input   start ] */
/* ---- Solve the linear system by tridiagonal algorithm ---- */array(ci, point_num);
array(di, point_num);for i:0 thru point_num step 1
doif i = 0 thenblock(ci[i] : aij[i,1]/aij[i,0],di[i] : bi[i]/aij[i,0])elseblock(tmp : aij[i,1]-aij[i,0]*ci[i-1],ci[i] : aij[i,2]/tmp,di[i] : (bi[i]-aij[i,0]*di[i-1])/tmp);
/* [wxMaxima: input   end   ] */
/* [wxMaxima: input   start ] */
/* -------------- Continue -------------- */for i:point_num thru 0 step -1
doif i = point_num thenkcoe[i] : di[i]elsekcoe[i] : di[i]-ci[i]*kcoe[i+1];
/* [wxMaxima: input   end   ] */
/* [wxMaxima: input   start ] */
/* -------------- Define the middle coefficients ------------- */array(a, point_num);
array(b, point_num);for i:1 thru point_num step 1
doblock(a[i] : kcoe[i-1]*(knots[i,0]-knots[i-1,0])-(knots[i,1]-knots[i-1,1]),b[i] : -kcoe[i]*(knots[i,0]-knots[i-1,0])+(knots[i,1]-knots[i-1,1]));
/* [wxMaxima: input   end   ] */
/* [wxMaxima: input   start ] */
/* ---------------- Get the interpolation functions ------------- */t : (x-knots[i-1,0])/(knots[i,0]-knots[i-1,0]);
define(fun1(i, x),(1-t)*knots[i-1,1]+t*knots[i,1]+t*(1-t)*(a[i]*(1-t)+b[i]*t)
);
/* [wxMaxima: input   end   ] */
/* [wxMaxima: input   start ] */
/* =============== The second variable spline ================ */
/* [wxMaxima: input   end   ] */
/* [wxMaxima: input   start ] */
/* -------------- Initiate the linear system --------------- */for i:0 thru point_num step 1
doif i = 0 thenblock(aij[i,0] : 2,aij[i,1] : 1,aij[i,2] : 0,bi[i] : 3*(knots[1,2]-knots[0,2])/(knots[1,0]-knots[0,0]))elseif i = point_num thenblock(n : point_num,aij[i,0] : 1,aij[i,1] : 2,aij[i,2] : 0,bi[i] : 3*(knots[n,2]-knots[n-1,2])/(knots[n,0]-knots[n-1,0]))elseblock(aij[i,0] : 1/(knots[i,0]-knots[i-1,0]),aij[i,1] : 2*(1/(knots[i,0]-knots[i-1,0])+1/(knots[i+1,0]-knots[i,0])),aij[i,2] : 1/(knots[i+1,0]-knots[i,0]),bi[i] : 3*((knots[i,2]-knots[i-1,2])/(knots[i,0]-knots[i-1,0])^2+(knots[i+1,2]-knots[i,2])/(knots[i+1,0]-knots[i,0])^2));
/* [wxMaxima: input   end   ] */
/* [wxMaxima: input   start ] */
/* ---- Solve the linear system by tridiagonal algorithm ---- */array(ci, point_num);
array(di, point_num);
for i:0 thru point_num step 1
doif i = 0 thenblock(ci[i] : aij[i,1]/aij[i,0],di[i] : bi[i]/aij[i,0])elseblock(tmp : aij[i,1]-aij[i,0]*ci[i-1],ci[i] : aij[i,2]/tmp,di[i] : (bi[i]-aij[i,0]*di[i-1])/tmp);
/* [wxMaxima: input   end   ] */
/* [wxMaxima: input   start ] */
/* -------------- Continue -------------- */for i:point_num thru 0 step -1
doif i = point_num thenkcoe[i] : di[i]elsekcoe[i] : di[i]-ci[i]*kcoe[i+1];
/* [wxMaxima: input   end   ] */
/* [wxMaxima: input   start ] */
/* -------------- Define the middle coefficients ------------- */array(a, point_num);
array(b, point_num);
for i:1 thru point_num step 1
doblock(a[i] : kcoe[i-1]*(knots[i,0]-knots[i-1,0])-(knots[i,2]-knots[i-1,2]),b[i] : -kcoe[i]*(knots[i,0]-knots[i-1,0])+(knots[i,2]-knots[i-1,2]));
/* [wxMaxima: input   end   ] */
/* [wxMaxima: input   start ] */
/* ---------------- Get the interpolation functions ------------- */s : (x-knots[i-1,0])/(knots[i,0]-knots[i-1,0]);
define(fun2(i, x),(1-s)*knots[i-1,2]+s*knots[i,2]+s*(1-s)*(a[i]*(1-s)+b[i]*s)
);
/* [wxMaxima: input   end   ] */
/* [wxMaxima: input   start ] */
/* ===================== Finalize and Plot ==================== */
/* [wxMaxima: input   end   ] */
/* [wxMaxima: input   start ] */
define(myx(x),
block(
if x < 1 thenfun1(1,x)
elseif x < 2 thenfun1(2,x)
elseif x < 3 thenfun1(3,x)
elseif x < 4 thenfun1(4,x)
elseif x < 5 thenfun1(5,x)
elseif x < 6 thenfun1(6,x)
elseif x < 7 thenfun1(7,x)
elsefun1(8,x)
));
/* [wxMaxima: input   end   ] */
/* [wxMaxima: input   start ] */
define(myy(x),
block(
if x < 1 thenfun2(1,x)
elseif x < 2 thenfun2(2,x)
elseif x < 3 thenfun2(3,x)
elseif x < 4 thenfun2(4,x)
elseif x < 5 thenfun2(5,x)
elseif x < 6 thenfun2(6,x)
elseif x < 7 thenfun2(7,x)
elsefun2(8,x)
));
/* [wxMaxima: input   end   ] */
/* [wxMaxima: input   start ] */
load(draw);
draw2d(parametric(myx(x),myy(x),x,knots[0,0],knots[point_num,0]));
/* [wxMaxima: input   end   ] */
/* [wxMaxima: input   start ] */
/* =================== Test and Check ================= */
/* [wxMaxima: input   end   ] */
/* [wxMaxima: input   start ] */
for i:0 thru point_num step 1
do display(knots[i,1], knots[i,2]);
/* [wxMaxima: input   end   ] */
/* [wxMaxima: input   start ] */
for i:0 thru 8 step 0.1
do display(myy(i));
/* [wxMaxima: input   end   ] */
/* [wxMaxima: input   start ] */
for i:0 thru 8 step 0.1
do display(myx(i));
/* [wxMaxima: input   end   ] */
/* [wxMaxima: input   start ] */
plot2d(myx(x),[x,0,8]);
/* [wxMaxima: input   end   ] */
/* [wxMaxima: input   start ] */
for i:0 thru point_num step 1
doblock(display(kcoe[i]),display(bi[i]),display(aij[i,0],aij[i,1],aij[i,2]),display(a[i], b[i]),print("-----"));
/* [wxMaxima: input   end   ] */