近来试用Maple 11的Physics Package，主要是关于量子力学计算的部份。不知道是我的水平低，还是Physics Package的功能局限，我觉得算一个么正变换，感觉用手算比用Maple 11算要方便……Maple 11那个Physics Package好像有许多古古怪怪的问题，我也搞不清楚。以下计算在压缩态|s> = SD|0>下算符Y1 = a exp(-i * phi)) + a^+ exp(i * phi)的不确定度sqrt( - ^2)，其中S是压缩态算符，D为位移算符，|0>为真空态，a和a^+分别为湮灭和产生算符，phi为实数。我的计算过程非常繁琐，希望高人多多指点，给我一些简单的解决方法。相关计算过程请参考D.F.Walls, G.J.Milburn, Quantum Optics, (Springer-Verlag, Berlin, 1994).
> restart;
> with(Physics):
设S和D为量子算符：
> Setup(quantumoperator = {S, D});
定义湮灭和产生算符：
> am := Annihilation(n, 1); ap := Creation(n, 1);
                         [quantumoperator = {D, S}]
                                     a-
                                     a+
定义几个么正变换的结果：
> DamD := Dagger(D) . am . D = am + alpha;
> DapD := Dagger(D) . ap . D = ap + conjugate(alpha);
> SamS := Dagger(S) . am . S = am * cosh(r) - ap * exp(2 * I * phi) * sinh(r);
> SapS := Dagger(S) . ap . S = ap * cosh(r) - am * exp(-2 * I * phi) * sinh(r);
还要设定S的么正性：
> SUI := Dagger(S) . S = 1;
                     ((Dagger(D)) . a-) . D = a- + alpha
                     ((Dagger(D)) . a+) . D = a+ + alpha
        ((Dagger(S)) . a-) . S = a- cosh(r) - a+ exp(2 I phi) sinh(r)
       ((Dagger(S)) . a+) . S = a+ cosh(r) - a- exp(-2 I phi) sinh(r)
                             (Dagger(S)) . S = 1
定义真空态：
> VS := Ket(n, 0);
                                  Ket(n, 0)
定义Y1：
> Y1 := am * exp(-I * phi) + ap * exp(I * phi);
                       a- exp(-I phi) + a+ exp(I phi)
准备计算S^+ D^+ Y1 D S：
> tmp := Dagger(S) . (Dagger(D) . Y1 . D) . S;
        exp(-I phi) ((Dagger(S)) . (((Dagger(D)) . a-) . D)) . S + exp(I
phi) ((Dagger(S)) . (((Dagger(D)) . a+) . D)) . S
把相关么正变换的结果代入去：
> tmp := subs({DamD, DapD}, tmp);
              exp(-I phi) ((Dagger(S)) . (a- + alpha)) . S
                        /              /          \\   
           + exp(I phi) \(Dagger(S)) . \a+ + alpha// . S
> tmp := subs({SamS, SapS, SUI}, tmp);
       exp(-I phi) (a- cosh(r) - a+ exp(2 I phi) sinh(r) + alpha)
                       /                                             \
          + exp(I phi) \a+ cosh(r) - a- exp(-2 I phi) sinh(r) + alpha/
便计算完了：
> ry1 := Dagger(VS) . tmp . VS;
                    exp(-I phi) alpha + exp(I phi) alpha
的计算过程也是类似的，也要定义一大堆么正变换的结果，如要令D^+ a a D = D^+ a D D^+ a D，非常麻烦，计算结果太长的我就不作显示了。
> tmp := Dagger(S) . (Dagger(D) . (Y1 . Y1) . D) . S;
> DamapD := Dagger(D) . (am @ ap) . D = (Dagger(D) . am . D) . (Dagger(D) . ap . D);
> DapamD := Dagger(D) . (ap @ am) . D = (Dagger(D) . ap . D) . (Dagger(D) . am . D);
> DamamD := Dagger(D) . (am . am) . D = (Dagger(D) . am . D) . (Dagger(D) . am . D);
> DapapD := Dagger(D) . (ap . ap) . D = (Dagger(D) . ap . D) . (Dagger(D) . ap . D);
> tmp := subs({DamapD, DapamD, DamamD, DapapD}, tmp);
> tmp := subs({DamD, DapD}, tmp);
> SamapS := Dagger(S) . (am @ ap) . S = (Dagger(S) . am . S) . (Dagger(S) . ap . S);
> SapamS := Dagger(S) . (ap @ am) . S = (Dagger(S) . ap . S) . (Dagger(S) . am . S);
> SamamS := Dagger(S) . (am . am) . S = (Dagger(S) . am . S) . (Dagger(S) . am . S);
> SapapS := Dagger(S) . (ap . ap) . S = (Dagger(S) . ap . S) . (Dagger(S) . ap . S);
> tmp := subs({SamapS, SapamS, SamamS, SapapS}, tmp);
> tmp := subs({SamS, SapS, SUI}, tmp);
把算出来了：
> ry2 := Dagger(VS) . tmp . VS;
不确定度sqrt( - ^2)的计算就是好几个代数化简：
> expand(ry2 - ry1^2);
                                           2          2
                  -2 sinh(r) cosh(r) + sinh(r)  + cosh(r)
> sqrt(factor(%));
                                               (1/2)
                         /                  2 \    
                         \(cosh(r) - sinh(r)) /    
> convert(%, exp);
                                          (1/2)
                              /         2\    
                              \(exp(-r)) /    
> simplify(%) assuming r >= 0;
                                   exp(-r)
于是得出了和Walls书中相同的结果。可见，这种计算过程非常麻烦，主要原因有几个：１、不知道怎么定义算符的么正性，使得如S^+ S这种本该等于1的结果也要自己去做代入；２、不知道如何定义么正变换（或其他什么变换）的结果。如果能用Setup函数去定义D^+ a D的结果，并且能自动根据D的么正性得出D^+ a a D = D^+ a D D^+ a D，免得又要自己设定一个变量去做代入。这种人为代入很容易发生错误，而且设定这种代入式子也很麻烦；３、如果直接算，而不是先算出S^+ D^+ Y1 D S再代入<0|..|0>内，会发生代入不了那些么正变换结果的情况。具体把式子代入S^+ D^+ Y1 D S时，也会发生一些奇怪的代入不成功的情况，使得代入过程要分两步去做。

通过别人的“通风报信”，使用众所周知的方法，经过一个晚上的艰苦工作（是计算机，不是我），Maple 11已经入手了。出于爱好学习的目的，我们决定“试用”一下Maple 11的新包，Physics Package。顾名思义，Physics Package就是拿来处理物理问题的。做物理的人肯定会感觉到，Maple并不能方便地处理日常遇到的很多数学问题，如指定某变量的对易性，抽象矢量分析（即不确定具体表象）等等。如今，这些问题都可以用Maple 11的Physics Package来处理。以下简单地演示对易关系的定义和展开。

众所周知，量子力学中有对易关系[x, p] = x * p – p * x = i * hbar，其中x是位置，p是动量，i是虚数，hbar是Planck常数 / (2 * pi)。那么在Maple 11中该如何定义该对易关系呢？首先载入Physics Package：
> with(Physics);
  [*, ., Annihilation, AntiCommutator, Bra, Bracket, Check, Commutator, Coordinates, Creation, Dagger, Define, Dgamma, FeynmanDiagrams, Fundiff, Gtaylor, Intc, Inverse, Ket, KroneckerDelta, LeviCivita, Parameters, Parity, Projector, Psigma, Setup, Simplify, SpaceTimeVector, Trace, Vectors, ^, dAlembertian, d_, diff, g_]
接着使用该包的Setup函数，设定x和p为量子算符（quantumoperator={x,p}），再定义代数规则[x, p] = I * hbar（algebrarules={%Commutator(x,p)=I*hbar}）：
> Setup(quantumoperator = {x, p}, algebrarules = {%Commutator(x, p) = I*hbar});
                          [algebrarules = {[x, p][-] = hbar I}, quantumoperator = {x, p}]
现在可以测试了。使用函数Commutator算对易式[x, p]，明显有以下的结果：
> %Commutator(x, p) = Commutator(x, p);
                                                 [x, p][-] = hbar I
注意加上“%”号在函数Commutator前，是为了使其处于“inert”的状态，即返回式子，不作具体计算。要计算处于“inert”状态的式子，可使用value或expand函数。我们又可以算[p, x]：
> %Commutator(p, x) = Commutator(p, x);
                                                [p, x][-] = -I hbar
结果刚好多了个负号，是正确的。最后，尝试嵌套地计算对易式，如计算[[x, p^2], x]：
> %Commutator(%Commutator(x, p^2), x);
                                                       2
                                                 [[x, p ][-], x][-]
接着使用Physics Package中的Expand函数对非对易项做展开，注意该函数必须写成长形式，即Physics:-Expand。
> Physics:-Expand(%);
                                                            2
                                                      2 hbar

很久没上来这了，忙着考试，也有地震光缆中断的原因。

刚刚去Maple的官网，发现Maple 11快出来了。其中注意到有一个新包，是关于物理的。下面是官网对该包的介绍：

Matlab中的ODE Solvers只能解一阶微分方程（组），如我们遇到高阶的常微分方程时，必须把其重写为一阶微分方程组，才能使用Matlab求解。
如果方程的形式为y'' = f(t, y, y')，我们可以通过y1 = y，y2 = y'代入原方程，使原方程化为：
y1' = y2
y2' = f(t, y1, y2)
从而可以使用Matlab的ODE Solvers求解。
这其中的重写过程我们可以通过调用Maple的函数来完成。我们使用Matlab帮助中的一个例子，即van der Pol Equation来重写（Example: Solving an IVP ODE (van der Pol Equation, Nonstiff)）。其中van der Pol Equation为：
y(t)'' - mu(1 - y(t)^2)y(t)' + y(t) = 0
我们现在把它重写为两条一阶常微分方程组。首先载入Maple的DEtools：
> maple('with(DEtools)');
然后使用Maple的函数convertsys来重写方程：
> maple('convertsys(diff(y(t), t, t) - mu * (1 - y(t)^2) * diff(y(t), t) + y(t) = 0, [], y(t), t)')
ans =
[[YP[1] = Y[2], YP[2] = mu*(1-Y[1]^2)*Y[2]-Y[1]], [Y[1] = y(t), Y[2] = diff(y(t),t)], NaN, []]
在方程形式为yp(n) = f(t, y, y'...yp(n-1))的情况下（yp(n)为y对t的n次导数），一般convertsys函数后两个参数不用变化。其中convertsys函数中的第二个参数为原方程初值条件，设定的形式为y(0) = ...，D(y)(0) = ...。
从ans中的第一项可以看出结果为：
YP[1] = Y[2]
YP[2] = mu*(1-Y[1]^2)*Y[2]-Y[1]
这便是两条一阶常微分方程组了，YP[1]即y1对t的一阶导数，Y[1]即是y1。
而ans中的第二项便是y1和y2的定义，即：
Y[1] = y(t)
Y[2] = diff(y(t),t)
经过重写，便可使用Matlab作数值求解了。
另外，Matlab 7中新增了ode15i函数，对“隐式”的“一阶”常微分方程直接做数值求解。
相关内容可参加Matlab帮助中的"Initial Value Problems for ODEs and DAEs"和Maple帮助中的"DEtools[convertsys]"。

如果要Maple一启动就算动执行某些命令，如你希望Maple每次启动时虚数单位都使用"i"而非"I"时，可以使用初始化档案maple.ini。
在默认安装时，在C:\Program Files\Maple 10\LIB\中建立maple.ini文件，以记事本打开，里面写上相关的语句，那么以后Maple每次运行时，都会首先执行里面的命令。
所以，我们在maple.ini内写上：
#Modify the imaginary unit to "i"
interface(imaginaryunit=i):
然后运行Maple，便自动运行了“interface(imaginaryunit=i):”这句命令。我们通常在命后加上冒号，免得启动时出现执行结果。现在我们可以测试一下：
>(1+i)*(1-i);
                                                   2
明显，"i"已经是虚数单位了，如要回复原状，可令：
>interface(imaginaryunit=I):
>(1+i)*(1-i);
                                             (1+i)(1-i)
现在"i"已不是虚数单位了，展开上式：
>expand(%);
                                                 1-i^2

如果你觉得Maple启动得太慢，占记忆体太多，或要进行很复杂的计算时，你可以考虑使用Command-line Maple。其启动方法为Start->All Programs->Maple 10>Command-line Maple 10或Maple安装目录下的bin.win\cmaple.exe。
Command-line Maple就是在命令行的形式下使用Maple，其指令的输入基本上和Standard Worksheet中一致，但也有少许不同。如Maple 10中的Standard Worksheet可以在命令后不加分号，在Command-line Maple中则不可以。
在Command-line Maple中可以运行大部份在Standard Worksheet中的命令，算式用prettyprint=1的形式表达（上回说到的），而且可以画图。在画图方面，二维图可以较好地表达出来，尽管也能画三维图，但模样较差。另外，Command-line Maple有些命令不能使用，如animate这种动态画图命令。
其中在Command-line Maple中画图时，其画图的device为char。即在Standard Worksheet时，令：
>interface(plotdevice = char):
此时在Standard Worksheet中画的图，即和在Command-line Maple中的一致。如：

>plot(sin(x), x=0..2*Pi);

这时，你就可以在各大讨论区中和别人讨论问题时，方便地加上“图”。要注意的是，要令图更整齐地显示，要把字体设置为Courier New。此时，英文字和相关符号能严格对齐。
如要回复默认画图device，可令：
>interface(plotdevice = default): 
详细设置可参考帮助：
>?interface

今天在google上的maple groups中看到一个小技巧

使用Classic Worksheet,输入公式执行,把输入和输出一起copy,就可以得到如下结果:

可见结果还是蛮漂亮的,只是有一点错位,改过来不难.

 

ref. http://groups.google.com/group/comp.soft-sys.math.maple/browse_thread/thread/112764cb26ba6cdd/c2cdfee874039c62#c2cdfee874039c62