一款开源计算机代数系统（CAS）Maxima

最新推荐文章于 2024-11-05 19:59:11 发布

转载最新推荐文章于 2024-11-05 19:59:11 发布 · 1.1k 阅读

1 ·

CC 4.0 BY-SA版权

原文链接：https://my.oschina.net/propagator/blog/798502

文章标签：

#java #操作系统 #python

本文介绍了使用Maxima进行数学计算的几个示例，包括计算Dirac γ矩阵的迹、求解线性方程组以及三次样条拟合的自然边界条件。通过具体的代码示例展示了如何在Maxima和wxMaxima环境中进行这些计算。

2019独角兽企业重金招聘Python工程师标准>>>

计算机代数系统有很多，功能也侧重于不同的方面。最有名的当属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$ maxima

Maxima 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 system
           A.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 
do
    block(
    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
do
    if i = 0 then
        block(
        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 then
        block(
        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])
        )
    else
        block(
        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
do
    if i = 0 then
        block(
        ci[i] : aij[i,1]/aij[i,0],
        di[i] : bi[i]/aij[i,0]
        )
    else
        block(
        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
do
    if i = point_num then
        kcoe[i] : di[i]
    else
        kcoe[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
do
    block(
    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
do
    if i = 0 then
        block(
        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 then
        block(
        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])
        )
    else
        block(
        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
do
    if i = 0 then
        block(
        ci[i] : aij[i,1]/aij[i,0],
        di[i] : bi[i]/aij[i,0]
        )
    else
        block(
        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
do
    if i = point_num then
        kcoe[i] : di[i]
    else
        kcoe[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
do
    block(
    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 then
    fun1(1,x)
elseif x < 2 then
    fun1(2,x)
elseif x < 3 then
    fun1(3,x)
elseif x < 4 then
    fun1(4,x)
elseif x < 5 then
    fun1(5,x)
elseif x < 6 then
    fun1(6,x)
elseif x < 7 then
    fun1(7,x)
else
    fun1(8,x)
));
/* [wxMaxima: input   end   ] */
/* [wxMaxima: input   start ] */
define(myy(x),
block(
if x < 1 then
    fun2(1,x)
elseif x < 2 then
    fun2(2,x)
elseif x < 3 then
    fun2(3,x)
elseif x < 4 then
    fun2(4,x)
elseif x < 5 then
    fun2(5,x)
elseif x < 6 then
    fun2(6,x)
elseif x < 7 then
    fun2(7,x)
else
    fun2(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
do
    block(
        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   ] */

转载于:https://my.oschina.net/propagator/blog/798502