-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathoperator_algebra.mac
245 lines (191 loc) · 8.9 KB
/
operator_algebra.mac
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
/* Maxima code for working with operators, especially quantum mechanical operators.
Barton Willis
Professor Emeritus of Mathematics
University of Nebraska at Kearney
Copyright © Barton Willis, 2021, 2022, 2024
GNU GENERAL PUBLIC LICENSE Version 3, 29 June 2007 */
/* There is no advantage to translating this code. But should a user translate this
file, the following eliminates some warnings from the translator. */
eval_when(translate,
declare_translated(operator_formula,dot_form,operator_apply,operator_simp,operator_express))$
load("simplifying");
load("operator_algebra.lisp");
/* Tell Maxima that 'operator' is a feature.*/
declare(operator, feature);
/* Utilities */
/* Return true iff the input is a declared operator or if the input has the form
U[a,b,...], where U has been declared an operator. */
operatorp(e) := (mapatom(e) and featurep(e,'operator)) or (subvarp(e) and featurep(inpart(e,0) ,'operator));
/* Predicates for detecting expressions whose operator is +,*,.,and ^^.
The function inpart looks at the internal form of an expression, so these
predicates work the same regardless of the setting of inflag.*/
plusp(e) := (not mapatom(e)) and inpart(e,0) = "+";
timesp(e) := (not mapatom(e)) and inpart(e,0) = "*";
exptp(e) := (not mapatom(e)) and inpart(e,0) = "^";
nctimesp(e) := (not mapatom(e)) and inpart(e,0) = ".";
ncexptp(e) := (not mapatom(e)) and inpart(e,0) = "^^";
operator_adjointp(e) := (not mapatom(e)) and inpart(e,0) = 'operator_adjoint;
/* Return the formula (a lambda form) for an operator fn. If no formula
has been declared or if fn isn't an operator, return false. */
get_operator_formula(fn) :=
if symbolp(fn) and operatorp(fn) then (
get(fn, 'formula))
else false;
/* Convert an expression e from a functional form to a "dot" form; for example,
if F & G are operators, then
dot_form(F(G(x)) --> F . G . x and dot_form(107*F(G(x)) --> 107*F^^2 . x */
dot_form(e) := block([inflag : true],
if plusp(e) or listp(e) then map('dot_form, e)
elseif timesp(e) and constantp(first(e)) then first(e)*dot_form(rest(e))
elseif not mapatom(e) and operatorp(inpart(e,0)) then inpart(e,0). dot_form(first(e))
else e);
operator_apply_error(op, arg) :=
if mapatom(arg) and get(arg,'internal) then (
error("I don't know how to apply ", op))
else (
error("I don't know how to apply ", op," to ", arg));
/* Convert from "dot" form to a functional form; for example, if Jx & Jy are
operators, operator_apply(Jx.Jy, ψ) --> Jx(Jy(ψ)). */
operator_apply(e,ψ) := block([opsubst : true,n, inflag : true,ee],
if listp(e) then (
map(lambda([s], operator_apply(s,ψ)),e))
elseif subvarp(e) and not mapatom(e) and operatorp(inpart(e,0)) then (
funmake(e,[ψ]))
elseif constantp(e) then (
e*ψ)
elseif operatorp(e) then (
funmake(e,[ψ]))
elseif plusp(e) then (
map(lambda([q], operator_apply(q,ψ)), e))
/* constant * operator is OK, otherwise, error. */
elseif timesp(e) then (
if constantp(first(e)) then (
operator_apply(first(e), operator_apply(rest(e),ψ)))
else (
operator_apply_error(e,ψ)))
elseif ncexptp(e) and plusp(first(e)) then (
[e,n] : args(e),
if integerp(n) and n > 0 then (
if evenp(n) then (
ee : operator_simp(e^^(n/2)),
ee : operator_apply(ee, operator_apply(ee, ψ)))
else (
ee : operator_simp(e^^((n+1)/2)),
operator_apply(operator_simp(e^^((n-1)/2)),operator_apply(ee, ψ))))
else operator_apply_error(e^^n,ψ))
elseif ncexptp(e) then (
[e,n] : args(e),
if integerp(n) and n > 0 then (
operator_apply(e, operator_apply(e^^(n-1), ψ)))
else ((e^^n)(ψ)))
elseif nctimesp(e) then (
operator_apply(first(e), operator_apply(rest(e), ψ)))
/* Catches, for example, F^n, F/G, ... */
else operator_apply_error(e,ψ));
/* For an expression in dot form, apply the formulae for the operators to evaluate
the expression e. Example
(%i3) ham : p.p/2 + q^^2/2$
(%i4) operator_apply(ham ,exp(- x^2/(2*ħ)))$
(%i5) operator_express(%),factor;
(%o5) (ħ*%e^(-x^2/(2*ħ)))/2 */
operator_express(e) := block([opsubst : true,fn],
if mapatom(e) then e
else (
fn : get_operator_formula(inpart(e,0)),
if fn # false then (
apply(fn,[operator_express(first(e))]))
else
map('operator_express, e)));
simp_op_adjoint(e) := block([estar],
/* e is an operator--look to see if it has a declared operator adjoint */
if operatorp(e) then (
estar : get(e,'operator_adjoint),
if estar = false then simpfuncall('operator_adjoint,e) else estar)
/* operator_adjoint(maptom) = conjugate(mapatom)*/
elseif mapatom(e) or constantp(e) then (
conjugate(e))
/* for a matrix, map operator_adjoint onto the matrix & transpose */
elseif matrixp(e) then (
block([matrix_element_transpose : 'operator_adjoint], transpose(e)))
/* operator_adjoint is an involution */
elseif operator_adjointp(e) then (
first(e))
/* operator_adjointp is additive and multiplicative */
elseif plusp(e) or timesp(e) then (
map('operator_adjoint, e))
/* operator_adjointp(a.b) = operator_adjointp(b).operator_adjointp(a) */
elseif nctimesp(e) then (
apply(".", map('operator_adjoint, reverse(args(e)))))
/* operator_adjointp(a^b) = operator_adjointp(a)^b, where 'b' is a declared integer */
elseif exptp(e) and featurep(second(e), 'integer) then (
operator_adjoint(first(e))^second(e))
/* operator_adjoint(a^^b) = operator_adjoint(a)^^b, where 'b' is a declared integer */
elseif ncexptp(e) and featurep(second(e), 'integer) then (
operator_adjoint(first(e))^^second(e))
/* what did I forget? Not sure--return a operator_adjointp nounform for all others */
else (simpfuncall('operator_adjoint, e)))$
simplifying('operator_adjoint,'simp_op_adjoint)$
/* For an operator in dot form, return a simplified version of e. The operator
e is applied to a gensym, that is converted to a functional form and sent
through the simplifier, and finally, the result is converted by to a dot form.
To supress error messages involving the gensym variables, give 'f' the internal property. */
operator_simp(e) := block([f : gensym(), dotdistrib : true, ans],
put(f,true,'internal),
ans : subst(f=1, dot_form(operator_apply(e,f))),
rem(f,'internal),
ans);
/* Return the commutator of the operators f & g; that is, return f.g - g.f. */
commutator(f,g) := operator_simp(f.g - g.f);
extended_operatorp(e) := block([inflag : true],
operatorp(e) or
constantp(e) or
(plusp(e) or nctimesp(e)) and every('extended_operatorp, args(e)) or
ncexptp(e) and operatorp(first(e)) and featurep(second(e), 'integer) or
timesp(e) and ((constantp(first(e)) and extended_operatorp(rest(e))) or
(operatorp(first(e)) and constantp(rest(args(e))))));
/* Return the n-th order (n < 5) Baker-Campbell-Hausdorf expansion; thus, find the first few
terms of a series expansion for 'z', where 'z' is a solution for exp(x)exp(y) = exp(z).
I haven't attempted to extend this to n >= 5, but it's possible.
The first and second arguments to 'bch' must be operators. */
bch(x,y,n) := block([ans : 0],
if not extended_operatorp(x) then (
error("The first argument to 'bch' must be an operator; found ",x)),
if not extended_operatorp(y) then (
error("The second argument to 'bch' must be an operator; found ",y)),
if n <= 0 or n >= 5 then (
error("Only able to compute the order four or lower BCH formula; requested order ",n)),
if n >= 1 then (
ans : ans + x + y),
if n >= 2 then (
ans : ans + commutator(x,y)/2),
if n >= 3 then (
ans : ans + (commutator(x,commutator(x,y)) + commutator(y, commutator(y,x)))/12),
if n >= 4 then (
ans : ans - commutator(y, commutator(x, commutator(x,y)))/24),
ans);
bch_extended(l, n) := block([ans : 0],
if n <= 0 or n >= 5 then (
error("Only able to compute the order four or lower BCH formula; requested order ",n)),
if not listp(l) then (
error("The first argument to bch_extended must be a Maxima list; found ",l)),
if not every('extended_operatorp, l) then (
error("The first agument to bch_extended must be a list of operators; found ",l)),
lreduce(lambda([f,g], bch(f,g,n)),l,0));
/* Terribly slow!
bch_alternative(x,y,n) := block([t, expx : 0, expy : 0, expz : 0, zpoly : 0, z : gensym(),eqs,sol : []],
n : n+1,
for k : 0 thru n do (
zpoly : zpoly + concat(z, k) * t^k),
declare(t, constant),
for k : 0 thru n do (
expx : expx + t^k * x^^k / k!,
expy : expy + t^k * y^^k / k!,
expz : expz + t^k * zpoly^^k / k!),
eqs : expand(expx . expy - expz),
for k : 1 thru n-1 do (
eq : coeff(eqs,t,k),
sk : solve(eq, concat(z,k-1)),
eqs : expand(subst(sk,eqs)),
push(first(sk), sol)),
rem(t,constant),
xreduce("+", map('rhs, sol))); */