maxima之自定义函数库

来源:互联网 发布:2222becom域名升级i 编辑:程序博客网 时间:2024/05/13 03:20

据说maxima是个能够代替maple或mathematica的开源数学软件,使用过后发现解决一般的问题还可以,如果要解决专业一点的问题,恐怕就需要自己来编写函数库了,由于资料严重缺乏,只能做一步看一步。

第一步:让maxima找到自定义函数库

编写自己的函数库之前,首先做的事就是让maxima能够找到你写的函数库,这需要修改变量file_search_maxima,这是个list类型的变量,所以要用append函数来添加新的路径。在终端里打开maxima,输入命令file_search_maxima; 系统输出:

(%o1) [/home/zy/.maxima/###.{mac,mc},
/usr/local/share/maxima/5.20.1/share/###.{mac,mc}, /usr/local/share/maxima/5.2/
0.1/share/{affine,algebra,algebra/charsets,algebra/solver,calculus,colnew,coln/
ew/lisp,combinatorics,contrib,contrib/altsimp,contrib/amatrix,contrib/bitwise,/
contrib/boolsimp,contrib/descriptive,contrib/diffequations,contrib/diffequatio/
ns/tests,contrib/distrib,contrib/ezunits,contrib/finance,contrib/format,contri/
b/fourier_elim,contrib/fractals,contrib/fresnel,contrib/gentran,contrib/gentra/
n/man,contrib/gentran/test,contrib/gf,contrib/graphs,contrib/Grobner,contrib/i/
ntegration,contrib/levin,contrib/lurkmathml,contrib/maximaMathML,contrib/mccli/
m,contrib/namespaces,contrib/noninteractive,contrib/nset,contrib/numericalio,c/
ontrib/pdiff,contrib/prim,contrib/rand,contrib/sarag,contrib/simplex,contrib/s/
implex/Tests,contrib/solve_rec,contrib/state,contrib/stats,contrib/stringproc,/
contrib/unit,contrib/vector3d,contrib/Zeilberger,diffequations,diff_form,draw,/
dynamics,hypergeometric,integequations,integration,lapack,lapack/blas,lapack/l/
apack,lbfgs,linearalgebra,lisp-utils,macro,matrix,minpack,minpack/lisp,misc,my/
lib,numeric,orthopoly,physics,simplification,sym,tensor,trigonometry,utils,vec/
tor}/###.{mac,mc}]

 

猜测###为通配符,试着输入file_search_maxima:append(file_search_maxima,["/home/zy/maxima/###.{mac,mc}"]); 这个意思是把/home/zy/maxima作为新的函数库路径加入,系统输出:

(%o2) [/home/zy/.maxima/###.{mac,mc},
/usr/local/share/maxima/5.20.1/share/###.{mac,mc}, /usr/local/share/maxima/5.2/
0.1/share/{affine,algebra,algebra/charsets,algebra/solver,calculus,colnew,coln/
ew/lisp,combinatorics,contrib,contrib/altsimp,contrib/amatrix,contrib/bitwise,/
contrib/boolsimp,contrib/descriptive,contrib/diffequations,contrib/diffequatio/
ns/tests,contrib/distrib,contrib/ezunits,contrib/finance,contrib/format,contri/
b/fourier_elim,contrib/fractals,contrib/fresnel,contrib/gentran,contrib/gentra/
n/man,contrib/gentran/test,contrib/gf,contrib/graphs,contrib/Grobner,contrib/i/
ntegration,contrib/levin,contrib/lurkmathml,contrib/maximaMathML,contrib/mccli/
m,contrib/namespaces,contrib/noninteractive,contrib/nset,contrib/numericalio,c/
ontrib/pdiff,contrib/prim,contrib/rand,contrib/sarag,contrib/simplex,contrib/s/
implex/Tests,contrib/solve_rec,contrib/state,contrib/stats,contrib/stringproc,/
contrib/unit,contrib/vector3d,contrib/Zeilberger,diffequations,diff_form,draw,/
dynamics,hypergeometric,integequations,integration,lapack,lapack/blas,lapack/l/
apack,lbfgs,linearalgebra,lisp-utils,macro,matrix,minpack,minpack/lisp,misc,my/
lib,numeric,orthopoly,physics,simplification,sym,tensor,trigonometry,utils,vec/
tor}/###.{mac,mc}, /home/zy/maxima/###.{mac,mc}]

果然新的函数库路径被加入了。

但是这个路径没有被保存,如果重启maxima,则file_search_maxima变量恢复原样,为了让系统能够保存,可以在~/.maxima/maxima-init.mac里加入上面的append命令,maxima启动的时候就会自动执行,从而修改file_search_maxima变量。

在终端下输入cd ~/,进入我的帐号主文件夹,发现没有.maxima文件夹,这没关系,自己创建一个就行了。

终端下输入mkdir .maxima,创建了.maxima目录,进去后用emacs创建maxima-init.mac文件,写入file_search_maxima:append(file_search_maxima,["/home/zy/maxima/###.{mac,mc}"]); 保存后启动maxima,发现新的函数库路径被加入。

试着在/home/zy/maxima下创建文件a.mac,并在maxima命令提示符下输入load(a); 系统输出

(%o2)                        /home/zy/maxima/a.mac

说明系统根据新的函数库路径找到了新创建的函数库a.mac

第二步:创建函数库

按照我的理解,当用load命令载入函数库的时候,maxima就像执行脚本一样从头到尾把mac文件里所有的语句都执行了一遍,如果语句可以立即执行,比如赋值语句,那么就当场执行,如果是定义函数,那就解释了放在内存里准备随时调用。所以,函数库真没什么特别的,只要写函数的定义就行了。

在maxima里,定义函数的语法是

fun(x1,x2,...,xn):=(exp1,exp2,...,expm);

其中fun是函数名,x1,x2,...,xn是参数名,exp1,exp2,....,expm是函数体中的m个语句,expm的值作为函数的返回值。

如果返回的语句写在中间,则要用以下形式:

fun(x1,x2,....,xn):=block([a:a0,b,c:c0],exp1,exp2,...,if a>0 then return(exps),...,expm);

注意其中return(exps)不可少,block后的[a:a0,b,c:c0]是局部变量,而且可以用a:a0的方式来赋初值。

要注意的还有maxima的语法和lisp类似,也是函数型语言,所以这个定义中的expx可以看作参数,中间用的是逗号而不是分号,expm尾不需要加逗号,定义末尾也可以用$来代替分号。

最后给个例子来说明一下,这个函数从拉格朗日函数得到拉格郎日方程

 

/*equation: 把拉格朗日函数转换成拉格朗日方程*/
/*其中position velocity都是表*/
/*以下为辅助函数*/
/*getequ: 用递归的方法来生成拉格朗日方程*/
getequ(lagrange,pos,vel):=
block([equlist:[],p,v],
    if not(pos=[] and vel=[]) then
    (
        p:first(pos),
        v:first(vel),
        equlist:[diff(diff(lagrange,v),t)=diff(lagrange,p)],
        equlist:append(equlist,getequ(lagrange,rest(pos),rest(vel)))
       
    ),
    return(equlist)
)$
   
equation(lagrange,position,velocity):=
block(
    if listp(position) and listp(velocity) then
    (
       
        if (length(position)=length(velocity)) then
           (
            return(getequ(lagrange,position,velocity))
           )
        else
           (
            print("the length of position must equal to velocity")
           )
    )
    else
    (
        print("position and velocity must be list")
    )
)$