Items tagged with mechanics


This presentation is on an undergrad intermediate Quantum Mechanics topic. Tackling the problem within a computer algebra worksheet in the way shown below is actually the novelty, using the Physics package to formulate the problem with quantum operators and related algebra rules in tensor notation.


Quantization of the Lorentz Force


Pascal Szriftgiser1 and Edgardo S. Cheb-Terrab2 

(1) Laboratoire PhLAM, UMR CNRS 8523, Université Lille 1, F-59655, France

(2) Maplesoft


We consider the case of a quantum, non-relativistic, particle with mass m and charge q evolving under the action of an arbitrary time-independent magnetic field "B=Curl(A(x,y,z)), "where `#mover(mi("A",mathcolor = "olive"),mo("→"))` is the vector potential. The Hamiltonian for this system is

H = (`#mover(mi("p",mathcolor = "olive"),mo("→"))`-q*`#mover(mi("A",mathcolor = "olive"),mo("→"))`(X))^2/(2*m)

where `#mover(mi("p",mathcolor = "olive"),mo("→"))` is the momentum of the particle, and the force acting in this particle, also called the Lorentz force, is given by


`#mover(mi("F",mathcolor = "olive"),mo("→"))` = m*(diff(v(t), t))


where `#mover(mi("v",mathcolor = "olive"),mo("→"))` is the quantized velocity of the particle, and all of  H, `#mover(mi("p",mathcolor = "olive"),mo("→"))`, `#mover(mi("v",mathcolor = "olive"),mo("→"))`, `#mover(mi("B",mathcolor = "olive"),mo("→"))`, `#mover(mi("A",mathcolor = "olive"),mo("→"))` and `#mover(mi("F",mathcolor = "olive"),mo("→"))` are Hermitian quantum operators representing observable quantities.


In the classic (non-quantum) case, `#mover(mi("F"),mo("→"))` for such a particle in the absence of electrical field is given by


`#mover(mi("F"),mo("→"))` = `&x`(q*`#mover(mi("v"),mo("→"))`, `#mover(mi("B"),mo("→"))`) ,


Problem: Departing from the Hamiltonian, show that in the quantum case the Lorentz force is given by [1]


`#mover(mi("F",mathcolor = "olive"),mo("→"))` = (1/2)*q*(`&x`(`#mover(mi("v",mathcolor = "olive"),mo("→"))`, `#mover(mi("B",mathcolor = "olive"),mo("→"))`)-`&x`(`#mover(mi("B",mathcolor = "olive"),mo("→"))`, `#mover(mi("v",mathcolor = "olive"),mo("→"))`))


[1] Photons et atomes, Introduction à l'électrodynamique quantique, p. 179, Claude Cohen-Tannoudji, Jacques Dupont-Roc et Gilbert Grynberg - EDP Sciences janvier 1987.




We choose to tackle the problem in Heisenberg's picture of quantum mechanices, where the state of a system is static and only the quantum operators evolve in time according to


diff(O(t), t) = I*Physics:-Commutator(H, O(t))/`ℏ`


Also, the algebraic manipulations are simpler using tensor abstract notation instead of the standard 3D vector notation. We then start setting the framework for the problem, a system of coordinates X, indicating the dimension of the tensor space to be 3 and the metric Euclidean, and that we will use lowercaselatin letters to represent tensor indices. In addition, not necessary but for convenience, we set the lowercase latin i to represent the imaginary unit and we request automaticsimplification so that the output of everything comes automatically simplified in size.


restart; with(Physics); interface(imaginaryunit = i)

Setup(mathematicalnotation = true, automaticsimplification = true, coordinates = X, dimension = 3, metric = Euclidean, spacetimeindices = lowercaselatin, quiet)

[automaticsimplification = true, coordinatesystems = {X}, dimension = 3, mathematicalnotation = true, metric = {(1, 1) = 1, (2, 2) = 1, (3, 3) = 1}, spacetimeindices = lowercaselatin]



Next we indicate the letters we will use to represent the quantum operators with which we will work, and also the standard commutation rules between position and momentum, always the starting point when dealing with quantum mechanics problems


Setup(quantumoperators = {F}, hermitianoperators = {A, B, H, p, r, v, x}, realobjects = {`ℏ`, m, q}, algebrarules = {%Commutator(p[k], p[n]) = 0, %Commutator(x[k], p[l]) = I*`ℏ`*KroneckerDelta[k, l], %Commutator(x[k], x[l]) = 0})

[algebrarules = {%Commutator(p[k], p[n]) = 0, %Commutator(x[k], p[l]) = I*`ℏ`*Physics:-KroneckerDelta[k, l], %Commutator(x[k], x[l]) = 0}, hermitianoperators = {A, B, H, p, r, v, x}, quantumoperators = {A, B, F, H, p, r, v, x}, realobjects = {`ℏ`, m, q, x1, x2, x3, %dAlembertian, Physics:-dAlembertian}]



Note that we start not indicating F as Hermitian, in order to arrive at that result. The quantum operators A, B, and F are explicit functions of X, so to avoid redundant display of this functionality on the screen we use


CompactDisplay((A, B, F)(X))

A(x1, x2, x3)*`will now be displayed as`*A


B(x1, x2, x3)*`will now be displayed as`*B


F(x1, x2, x3)*`will now be displayed as`*F


Define now as tensors the quantum operators that we will use with tensorial notation (recalling: for these, Einstein's sum rule for repeated indices will be automatically applied when simplifying)


Define(x, p, v, A, B, F, quiet)

{A, B, F, p, v, x, Physics:-Dgamma[a], Physics:-Psigma[a], Physics:-d_[a], Physics:-g_[a, b], Physics:-KroneckerDelta[a, b], Physics:-LeviCivita[a, b, c], Physics:-SpaceTimeVector[a](X)}


The Hamiltonian,

H = (`#mover(mi("p",mathcolor = "olive"),mo("→"))`-q*`#mover(mi("A",mathcolor = "olive"),mo("→"))`(X))^2/(2*m)

in tensorial notation, is given by

H = (p[n]-q*A[n](X))^2/(2*m)

H = (1/2)*Physics:-`^`(p[n]-q*A[n](X), 2)/m


Generally speaking to arrive at  ```#mover(mi("F",mathcolor = "olive"),mo("→"))` = (1/2)*q*(`&x`(`#mover(mi("v",mathcolor = "olive"),mo("→"))`, `#mover(mi("B",mathcolor = "olive"),mo("→"))`)-`&x`(`#mover(mi("B",mathcolor = "olive"),mo("→"))`, `#mover(mi("v",mathcolor = "olive"),mo("→"))`)) what we now need to do is

1) Express this Hamiltonian (5) in terms of the velocity


And, recalling that, in Heisenberg's picture, quantum operators evolve in time according to

diff(O(t), t) = I*Physics:-Commutator(H, O(t))/`ℏ`


2) Take the commutator of H with the velocity itself to obtain its time derivative and, from `#mover(mi("F",mathcolor = "olive"),mo("→"))` = m*(diff(v(t), t)) , that commutator is already the force up to some constant factors.


To get in contact with the basic commutation rules between position and momentum behind quantum phenomena, the quantized velocity itself can be computed as the time derivative of the position operator, i.e as the commutator of x[k] with H

I*Commutator(H = (1/2)*Physics[`^`](p[n]-q*A[n](X), 2)/m, x[k])/`ℏ`

I*Physics:-Commutator(H, x[k])/`ℏ` = (1/2)*(I*q^2*Physics:-AntiCommutator(A[n](X), Physics:-Commutator(A[n](X), x[k]))-I*q*Physics:-AntiCommutator(p[n], Physics:-Commutator(A[n](X), x[k]))-2*(q*A[n](X)-p[n])*Physics:-KroneckerDelta[k, n]*`ℏ`)/(`ℏ`*m)


This expression for the velocity, that involves commutators between the potential A[n](X), the position x[k] and the momentum p[n], can be simplified taking into account the basic quantum algebra rules between position and momentum. We assume that A[n](X)(X) can be decomposed into a formal power series (possibly infinite) of the x[k], hence all the A[n](X) commute between themselves as well as with all the x[k]


{%Commutator(A[k](X), x[l]) = 0, %Commutator(A[k](X), A[l](X)) = 0}

{%Commutator(A[k](X), x[l]) = 0, %Commutator(A[k](X), A[l](X)) = 0}


(Note: in some cases, this is not true, but those cases are beyond the scope of this worksheet.)


Add these rules to the algebra rules already set so that they are all taken into account when simplifying things


Setup(algebrarules = {%Commutator(A[k](X), x[l]) = 0, %Commutator(A[k](X), A[l](X)) = 0})

[algebrarules = {%Commutator(p[k], p[n]) = 0, %Commutator(x[k], p[l]) = I*`ℏ`*Physics:-KroneckerDelta[k, l], %Commutator(x[k], x[l]) = 0, %Commutator(A[k](X), x[l]) = 0, %Commutator(A[k](X), A[l](X)) = 0}]


Simplify(I*Physics[Commutator](H, x[k])/`ℏ` = (1/2)*(I*q^2*Physics[AntiCommutator](A[n](X), Physics[Commutator](A[n](X), x[k]))-I*q*Physics[AntiCommutator](p[n], Physics[Commutator](A[n](X), x[k]))-2*(q*A[n](X)-p[n])*Physics[KroneckerDelta][k, n]*`ℏ`)/(`ℏ`*m))

I*Physics:-Commutator(H, x[k])/`ℏ` = (-A[k](X)*q+p[k])/m


The right-hand side of (9) is then the kth component of the velocity tensor quantum operator, the relationship is the same as in the classical case

v[k] = rhs(I*Physics[Commutator](H, x[k])/`ℏ` = (-A[k](X)*q+p[k])/m)

v[k] = (-A[k](X)*q+p[k])/m


and with this the Hamiltonian (5) can now be rewritten in term of the velocity completing step 1)

simplify(H = (1/2)*Physics[`^`](p[n]-q*A[n](X), 2)/m, {SubstituteTensorIndices(k = n, (rhs = lhs)(v[k] = (-A[k](X)*q+p[k])/m))})

H = (1/2)*m*Physics:-`^`(v[n], 2)


For step 2), to compute

 `#mover(mi("F",mathcolor = "olive"),mo("→"))` = m*(diff(v(t), t)) and m*(diff(v(t), t)) = I*m*Physics:-Commutator(H, v(t)[k])/`ℏ` 


we need the commutator between the different components of the quantized velocity which, contrary to what happens in the classical case, do not commute. For this purpose, take the commutator between (10) with itself after replacing the free index

Commutator(v[k] = (-A[k](X)*q+p[k])/m, SubstituteTensorIndices(k = n, v[k] = (-A[k](X)*q+p[k])/m))

Physics:-Commutator(v[k], v[n]) = -q*(Physics:-Commutator(A[k](X), p[n])+Physics:-Commutator(p[k], A[n](X)))/m^2


To simplify (12), we use the fact that if f  is a commutative mapping that can be decomposed into a formal power series in all the complex plan (which is assumed to be the case for all A[n](X)(X)), then

Physics:-Commutator(p[k], f(x, y, z)) = -I*`ℏ`*`∂`[k](f(x, y, z))

where p[k]"=-i `ℏ` `∂`[k] " is the momentum operator along the x[k] axis. This relation reads in tensor notation:

Commutator(p[k], A[n](X)) = -I*`ℏ`*d_[k](A[n](X))

Physics:-Commutator(p[k], A[n](X)) = -I*`ℏ`*Physics:-d_[k](A[n](X), [X])


Add this rule to the rules previously set in order to automatically take it into account in (12)

Setup(Physics[Commutator](p[k], A[n](X)) = -I*`ℏ`*Physics[d_][k](A[n](X), [X]))

[algebrarules = {%Commutator(p[k], p[n]) = 0, %Commutator(p[k], A[n](X)) = -I*`ℏ`*Physics:-d_[k](A[n](X), [X]), %Commutator(x[k], p[l]) = I*`ℏ`*Physics:-KroneckerDelta[k, l], %Commutator(x[k], x[l]) = 0, %Commutator(A[k](X), x[l]) = 0, %Commutator(A[k](X), A[l](X)) = 0}]


Physics[Commutator](v[k], v[n]) = -q*(Physics[Commutator](A[k](X), p[n])+Physics[Commutator](p[k], A[n](X)))/m^2

Physics:-Commutator(v[k], v[n]) = -I*q*`ℏ`*(Physics:-d_[n](A[k](X), [X])-Physics:-d_[k](A[n](X), [X]))/m^2


Also add this other rule so that it is taken into account automatically

Setup(Physics[Commutator](v[k], v[n]) = -I*q*`ℏ`*(Physics[d_][n](A[k](X), [X])-Physics[d_][k](A[n](X), [X]))/m^2)

[algebrarules = {%Commutator(p[k], p[n]) = 0, %Commutator(p[k], A[n](X)) = -I*`ℏ`*Physics:-d_[k](A[n](X), [X]), %Commutator(v[k], v[n]) = -I*q*`ℏ`*(Physics:-d_[n](A[k](X), [X])-Physics:-d_[k](A[n](X), [X]))/m^2, %Commutator(x[k], p[l]) = I*`ℏ`*Physics:-KroneckerDelta[k, l], %Commutator(x[k], x[l]) = 0, %Commutator(A[k](X), x[l]) = 0, %Commutator(A[k](X), A[l](X)) = 0}]


Recalling now the expression of the Hamiltonian (11) as a function of the velocity, one can compute the components of the force operator  "()Component(v*B,k)=m (v[k])=(i m [H,v[k]][-])/`ℏ`"

F[k](X) = I*m*%Commutator(rhs(H = (1/2)*m*Physics[`^`](v[n], 2)), v[k])/`ℏ`

F[k](X) = I*m*%Commutator((1/2)*m*Physics:-`^`(v[n], 2), v[k])/`ℏ`


Simplify this expression for the quantized force taking the quantum algebra rules (16) into account

Simplify(F[k](X) = I*m*%Commutator((1/2)*m*Physics[`^`](v[n], 2), v[k])/`ℏ`)

F[k](X) = (1/2)*q*(-Physics:-`*`(Physics:-d_[n](A[k](X), [X]), v[n])+Physics:-`*`(Physics:-d_[k](A[n](X), [X]), v[n])-Physics:-`*`(v[n], Physics:-d_[n](A[k](X), [X]))+Physics:-`*`(v[n], Physics:-d_[k](A[n](X), [X])))


It is not difficult to verify that this is the antisymmetrized vector product `&x`(`#mover(mi("v",mathcolor = "olive"),mo("→"))`, `#mover(mi("B",mathcolor = "olive"),mo("→"))`). Departing from `#mover(mi("B",mathcolor = "olive"),mo("→"))` = `&x`(VectorCalculus[Nabla], `#mover(mi("A",mathcolor = "olive"),mo("→"))`) expressed using tensor notation,

B[c](X) = LeviCivita[c, n, m]*d_[n](A[m](X))

B[c](X) = -Physics:-LeviCivita[c, m, n]*Physics:-d_[n](A[m](X), [X])


and taking into acount that

 Component(`&x`(`#mover(mi("v",mathcolor = "olive"),mo("→"))`, `#mover(mi("B",mathcolor = "olive"),mo("→"))`), k) = `ε`[b, c, k]*v[b]*B[c](X) 

multiply both sides of (19) by `ε`[b, c, k]*v[b], getting

LeviCivita[k, b, c]*v[b]*(B[c](X) = -Physics[LeviCivita][c, m, n]*Physics[d_][n](A[m](X), [X]))

Physics:-LeviCivita[b, c, k]*Physics:-`*`(v[b], B[c](X)) = -Physics:-LeviCivita[b, c, k]*Physics:-LeviCivita[c, m, n]*Physics:-`*`(v[b], Physics:-d_[n](A[m](X), [X]))


Simplify(Physics[LeviCivita][b, c, k]*Physics[`*`](v[b], B[c](X)) = -Physics[LeviCivita][b, c, k]*Physics[LeviCivita][c, m, n]*Physics[`*`](v[b], Physics[d_][n](A[m](X), [X])))

Physics:-LeviCivita[b, c, k]*Physics:-`*`(v[b], B[c](X)) = Physics:-`*`(v[m], Physics:-d_[k](A[m](X), [X]))-Physics:-`*`(v[n], Physics:-d_[n](A[k](X), [X]))


Finally, replacing the repeated index m by n 

SubstituteTensorIndices(m = n, Physics[LeviCivita][b, c, k]*Physics[`*`](v[b], B[c](X)) = Physics[`*`](v[m], Physics[d_][k](A[m](X), [X]))-Physics[`*`](v[n], Physics[d_][n](A[k](X), [X])))

Physics:-LeviCivita[b, c, k]*Physics:-`*`(v[b], B[c](X)) = Physics:-`*`(v[n], Physics:-d_[k](A[n](X), [X]))-Physics:-`*`(v[n], Physics:-d_[n](A[k](X), [X]))


Likewise, for

 Component(`&x`(`#mover(mi("v",mathcolor = "olive"),mo("→"))`, `#mover(mi("B",mathcolor = "olive"),mo("→"))`), k) = `ε`[b, c, k]*B[b]*B[c](X) 

multiplying (19), this time from the right instead of from the left, we get

Simplify(((B[c](X) = -Physics[LeviCivita][c, m, n]*Physics[d_][n](A[m](X), [X]))*LeviCivita[k, b, c])*v[b])

Physics:-LeviCivita[b, c, k]*Physics:-`*`(B[c](X), v[b]) = Physics:-`*`(Physics:-d_[k](A[m](X), [X]), v[m])-Physics:-`*`(Physics:-d_[n](A[k](X), [X]), v[n])


SubstituteTensorIndices(m = n, Physics[LeviCivita][b, c, k]*Physics[`*`](B[c](X), v[b]) = Physics[`*`](Physics[d_][k](A[m](X), [X]), v[m])-Physics[`*`](Physics[d_][n](A[k](X), [X]), v[n]))

Physics:-LeviCivita[b, c, k]*Physics:-`*`(B[c](X), v[b]) = Physics:-`*`(Physics:-d_[k](A[n](X), [X]), v[n])-Physics:-`*`(Physics:-d_[n](A[k](X), [X]), v[n])


Simplifying now the expression (18) for the quantized force taking into account (22) and (24) we get

simplify(F[k](X) = (1/2)*q*(-Physics[`*`](Physics[d_][n](A[k](X), [X]), v[n])+Physics[`*`](Physics[d_][k](A[n](X), [X]), v[n])-Physics[`*`](v[n], Physics[d_][n](A[k](X), [X]))+Physics[`*`](v[n], Physics[d_][k](A[n](X), [X]))), {(rhs = lhs)(Physics[LeviCivita][b, c, k]*Physics[`*`](v[b], B[c](X)) = Physics[`*`](v[n], Physics[d_][k](A[n](X), [X]))-Physics[`*`](v[n], Physics[d_][n](A[k](X), [X]))), (rhs = lhs)(Physics[LeviCivita][b, c, k]*Physics[`*`](B[c](X), v[b]) = Physics[`*`](Physics[d_][k](A[n](X), [X]), v[n])-Physics[`*`](Physics[d_][n](A[k](X), [X]), v[n]))})

F[k](X) = (1/2)*q*Physics:-LeviCivita[b, c, k]*(Physics:-`*`(v[b], B[c](X))+Physics:-`*`(B[c](X), v[b]))



`#mover(mi("F",mathcolor = "olive"),mo("→"))` = (1/2)*q*(`&x`(`#mover(mi("v",mathcolor = "olive"),mo("→"))`, `#mover(mi("B",mathcolor = "olive"),mo("→"))`)-`&x`(`#mover(mi("B",mathcolor = "olive"),mo("→"))`, `#mover(mi("v",mathcolor = "olive"),mo("→"))`))

in tensor notation. Finally, we note that this operator is Hermitian as expected

(F[k](X) = (1/2)*q*Physics[LeviCivita][b, c, k]*(Physics[`*`](v[b], B[c](X))+Physics[`*`](B[c](X), v[b])))-Dagger(F[k](X) = (1/2)*q*Physics[LeviCivita][b, c, k]*(Physics[`*`](v[b], B[c](X))+Physics[`*`](B[c](X), v[b])))

F[k](X)-Physics:-Dagger(F[k](X)) = 0


Download:,   Quantization_of_the_Lorentz_force.pdf

Edgardo S. Cheb-Terrab
Physics, Differential Equations and Mathematical Functions, Maplesoft
Editor, Computer Physics Communications

The presentation below is on undergrad Quantum Mechanics. Tackling this topic within a computer algebra worksheet in the way it's done below, however, is an exciting novelty and illustrates well the level of abstraction that is now possible using the Physics package.


Quantum Mechanics: Schrödinger vs Heisenberg picture

Pascal Szriftgiser1 and Edgardo S. Cheb-Terrab2 

(1) Laboratoire PhLAM, UMR CNRS 8523, Université Lille 1, F-59655, France

(2) Maplesoft


Within the Schrödinger picture of Quantum Mechanics, the time evolution of the state of a system, represented by a Ket "| psi(t) >", is determined by Schrödinger's equation:

I*`ℏ`*(diff(Ket(psi, t), t)) = H*Ket(psi, t)

where H, the Hamiltonian, as well as the quantum operators O__S representing observable quantities, are all time-independent.


Within the Heisenberg picture, a Ket Ket(psi, 0) representing the state of the system does not evolve with time, but the operators O__H(t)representing observable quantities, and through them the Hamiltonian H, do.


Problem: Departing from Schrödinger's equation,


a) Show that the expected value of a physical observable in Schrödinger's and Heisenberg's representations is the same, i.e. that

Bra(psi, t)*O__S*Ket(psi, t) = Bra(psi, 0)*O__H(t)*Ket(psi, 0)


b) Show that the evolution equation of an observable O__H in Heisenberg's picture, equivalent to Schrödinger's equation,  is given by:

diff(O__H(t), t) = (-I*Physics:-Commutator(O__H(t), H))*(1/`ℏ`)

where in the right-hand-side we see the commutator of O__H with the Hamiltonian of the system.

Solution: Let O__S and O__H respectively be operators representing one and the same observable quantity in Schrödinger's and Heisenberg's pictures, and H be the operator representing the Hamiltonian of a physical system. All of these operators are Hermitian. So we start by setting up the framework for this problem accordingly, including that the time t and Planck's constant are real. To automatically combine powers of the same base (happening frequently in what follows) we also set combinepowersofsamebase = true. The following input/output was obtained using the latest Physics update (Aug/31/2016) distributed on the Maplesoft R&D Physics webpage.


Physics:-Setup(hermitianoperators = {H, O__H, O__S}, realobjects = {`ℏ`, t}, combinepowersofsamebase = true, mathematicalnotation = true)

[combinepowersofsamebase = true, hermitianoperators = {H, O__H, O__S}, mathematicalnotation = true, realobjects = {`ℏ`, t}]


Let's consider Schrödinger's equation

I*`ℏ`*(diff(Ket(psi, t), t)) = H*Ket(psi, t)

I*`ℏ`*(diff(Physics:-Ket(psi, t), t)) = Physics:-`*`(H, Physics:-Ket(psi, t))


Now, H is time-independent, so (2) can be formally solved: psi(t) is obtained from the solution psi(0) at time t = 0, as follows:

T := exp(-I*H*t/`ℏ`)



Ket(psi, t) = T*Ket(psi, 0)

Physics:-Ket(psi, t) = Physics:-`*`(exp(-I*t*H/`ℏ`), Physics:-Ket(psi, 0))


To check that (4) is a solution of (2), substitute it in (2):

eval(I*`ℏ`*(diff(Physics[Ket](psi, t), t)) = Physics[`*`](H, Physics[Ket](psi, t)), Physics[Ket](psi, t) = Physics[`*`](exp(-I*H*t/`ℏ`), Physics[Ket](psi, 0)))

Physics:-`*`(H, exp(-I*t*H/`ℏ`), Physics:-Ket(psi, 0)) = Physics:-`*`(H, exp(-I*t*H/`ℏ`), Physics:-Ket(psi, 0))


Next, to relate the Schrödinger and Heisenberg representations of an Hermitian operator O representing an observable physical quantity, recall that the value expected for this quantity at time t during a measurement is given by the mean value of the corresponding operator (i.e., bracketing it with the state of the system Ket(psi, t)).

So let O__S be an observable in the Schrödinger picture: its mean value is obtained by bracketing the operator with equation (4):

Dagger(Ket(psi, t) = Physics[`*`](exp(-I*H*t/`ℏ`), Ket(psi, 0)))*O__S*(Ket(psi, t) = Physics[`*`](exp(-I*H*t/`ℏ`), Ket(psi, 0)))

Physics:-`*`(Physics:-Bra(psi, t), O__S, Physics:-Ket(psi, t)) = Physics:-`*`(Physics:-Bra(psi, 0), exp(I*t*H/`ℏ`), O__S, exp(-I*t*H/`ℏ`), Physics:-Ket(psi, 0))


The composed operator within the bracket on the right-hand-side is the operator O in Heisenberg's picture, O__H(t)

Dagger(T)*O__S*T = O__H(t)

Physics:-`*`(exp(I*t*H/`ℏ`), O__S, exp(-I*t*H/`ℏ`)) = O__H(t)


Analogously, inverting this equation,

(T*(Physics[`*`](exp(I*H*t/`ℏ`), O__S, exp(-I*H*t/`ℏ`)) = O__H(t)))*Dagger(T)

O__S = Physics:-`*`(exp(-I*t*H/`ℏ`), O__H(t), exp(I*t*H/`ℏ`))


As an aside to the problem, we note from these two equations, and since the operator T = exp((-I*H*t)*(1/`ℏ`)) is unitary (because H is Hermitian), that the switch between Schrödinger's and Heisenberg's pictures is accomplished through a unitary transformation.


Inserting now this value of O__S from (8) in the right-hand-side of (6), we get the answer to item a)

lhs(Physics[`*`](Bra(psi, t), O__S, Ket(psi, t)) = Physics[`*`](Bra(psi, 0), exp(I*H*t/`ℏ`), O__S, exp(-I*H*t/`ℏ`), Ket(psi, 0))) = eval(rhs(Physics[`*`](Bra(psi, t), O__S, Ket(psi, t)) = Physics[`*`](Bra(psi, 0), exp(I*H*t/`ℏ`), O__S, exp(-I*H*t/`ℏ`), Ket(psi, 0))), O__S = Physics[`*`](exp(-I*H*t/`ℏ`), O__H(t), exp(I*H*t/`ℏ`)))

Physics:-`*`(Physics:-Bra(psi, t), O__S, Physics:-Ket(psi, t)) = Physics:-`*`(Physics:-Bra(psi, 0), O__H(t), Physics:-Ket(psi, 0))


where, on the left-hand-side, the Ket representing the state of the system is evolving with time (Schrödinger's picture), while on the the right-hand-side the Ket `ψ__0`is constant and it is O__H(t), the operator representing an observable physical quantity, that evolves with time (Heisenberg picture). As expected, both pictures result in the same expected value for the physical quantity represented by O.


To complete item b), the derivation of the evolution equation for O__H(t), we take the time derivative of the equation (7):

diff((rhs = lhs)(Physics[`*`](exp(I*H*t/`ℏ`), O__S, exp(-I*H*t/`ℏ`)) = O__H(t)), t)

diff(O__H(t), t) = I*Physics:-`*`(H, exp(I*t*H/`ℏ`), O__S, exp(-I*t*H/`ℏ`))/`ℏ`-I*Physics:-`*`(exp(I*t*H/`ℏ`), O__S, H, exp(-I*t*H/`ℏ`))/`ℏ`


To rewrite this equation in terms of the commutator  Physics:-Commutator(O__S, H), it suffices to re-order the product  H  exp(I*H*t/`ℏ`) placing the exponential first:

Library:-SortProducts(diff(O__H(t), t) = I*Physics[`*`](H, exp(I*H*t/`ℏ`), O__S, exp(-I*H*t/`ℏ`))/`ℏ`-I*Physics[`*`](exp(I*H*t/`ℏ`), O__S, H, exp(-I*H*t/`ℏ`))/`ℏ`, [exp(I*H*t/`ℏ`), H], usecommutator)

diff(O__H(t), t) = I*Physics:-`*`(exp(I*t*H/`ℏ`), H, O__S, exp(-I*t*H/`ℏ`))/`ℏ`-I*Physics:-`*`(exp(I*t*H/`ℏ`), Physics:-`*`(H, O__S)+Physics:-Commutator(O__S, H), exp(-I*t*H/`ℏ`))/`ℏ`


Normal(diff(O__H(t), t) = I*Physics[`*`](exp(I*H*t/`ℏ`), H, O__S, exp(-I*H*t/`ℏ`))/`ℏ`-I*Physics[`*`](exp(I*H*t/`ℏ`), Physics[`*`](H, O__S)+Physics[Commutator](O__S, H), exp(-I*H*t/`ℏ`))/`ℏ`)

diff(O__H(t), t) = -I*Physics:-`*`(exp(I*t*H/`ℏ`), Physics:-Commutator(O__S, H), exp(-I*t*H/`ℏ`))/`ℏ`


Finally, to express the right-hand-side in terms of  Physics:-Commutator(O__H(t), H) instead of Physics:-Commutator(O__S, H), we take the commutator of the equation (8) with the Hamiltonian

Commutator(O__S = Physics[`*`](exp(-I*H*t/`ℏ`), O__H(t), exp(I*H*t/`ℏ`)), H)

Physics:-Commutator(O__S, H) = Physics:-`*`(exp(-I*t*H/`ℏ`), Physics:-Commutator(O__H(t), H), exp(I*t*H/`ℏ`))


Combining these two expressions, we arrive at the expected result for b), the evolution equation of a given observable O__H in Heisenberg's picture

eval(diff(O__H(t), t) = -I*Physics[`*`](exp(I*H*t/`ℏ`), Physics[Commutator](O__S, H), exp(-I*H*t/`ℏ`))/`ℏ`, Physics[Commutator](O__S, H) = Physics[`*`](exp(-I*H*t/`ℏ`), Physics[Commutator](O__H(t), H), exp(I*H*t/`ℏ`)))

diff(O__H(t), t) = -I*Physics:-Commutator(O__H(t), H)/`ℏ`


Download:     Schrodinger_vs_Heisenberg_picture.pdf

Edgardo S. Cheb-Terrab
Physics, Differential Equations and Mathematical Functions, Maplesoft
Editor, Computer Physics Communications

I am trying to duplicate the answer to a problem of the nutation of a spinning top. The problem is number 6.1.3 in the book Computer Algebra Recipes for Classical Mechanics by Richard Ens and George McGuire.

TopSC := `<,>`(5*sin(z)*exp(-(Pi-z)/(1.5)), 0, z) for z = 0..Pi

Spacecurve TopSC, when rotated about the z axis, generates the surface of revolution of the top.

The solution requires knowing the moments of inertia (MI) of the top about the z axis and the x or y axis.

Does the following integral correctly calculate the MI about the z axis?

int(2*Pi*TopSC[1]*(TopSC[1]^2), [z = 0 .. z, z = 0 .. Pi])

What sequence of Maple 2016 commands will calculate the top's MI about the x or y axis?

Below is the worksheet with the whole material presented yesterday in the webinar, “Applying the power of computer algebra to theoretical physics”, broadcasted by the “Institute of Physics” (IOP, England). The material was very well received, rated 4.5 out of 5 (around 30 voters among the more than 300 attendants), and generated a lot of feedback. The webinar was recorded so that it is possible to watch it (for free, of course, click the link above, it will ask you for registration, though, that’s how IOP works).

Anyway, you can reproduce the presentation with the worksheet below (mw file linked at the end, or the corresponding pdf also linked with all the input lines executed). As usual, to reproduce the input/output you need to have installed the latest version of Physics, available in the Maplesoft R&D Physics webpage.

Why computer algebra?




... and why computer algebra?

We can concentrate more on the ideas instead of on the algebraic manipulations


We can extend results with ease


We can explore the mathematics surrounding a problem


We can share results in a reproducible way


Representation issues that were preventing the use of computer algebra in Physics



Notation and related mathematical methods that were missing:

coordinate free representations for vectors and vectorial differential operators,

covariant tensors distinguished from contravariant tensors,

functional differentiation, relativity differential operators and sum rule for tensor contracted (repeated) indices

Bras, Kets, projectors and all related to Dirac's notation in Quantum Mechanics


Inert representations of operations, mathematical functions, and related typesetting were missing:


inert versus active representations for mathematical operations

ability to move from inert to active representations of computations and viceversa as necessary

hand-like style for entering computations and textbook-like notation for displaying results


Key elements of the computational domain of theoretical physics were missing:


ability to handle products and derivatives involving commutative, anticommutative and noncommutative variables and functions

ability to perform computations taking into account custom-defined algebra rules of different kinds

(commutator, anticommutator and bracket rules, etc.)





The Maple computer algebra environment


Classical Mechanics


Inertia tensor for a triatomic molecule


Classical Field Theory


*The field equations for the lambda*Phi^4 model


*Maxwell equations departing from the 4-dimensional Action for Electrodynamics


*The Gross-Pitaevskii field equations for a quantum system of identical particles


Quantum mechanics


*The quantum operator components of  `#mover(mi("L",mathcolor = "olive"),mo("&rarr;",fontstyle = "italic"))` satisfy "[L[j],L[k]][-]=i `&epsilon;`[j,k,m] L[m]"


Quantization of the energy of a particle in a magnetic field


Unitary Operators in Quantum Mechanics


*Eigenvalues of an unitary operator and exponential of Hermitian operators


Properties of unitary operators



Consider two set of kets " | a[n] >" and "| b[n] >", each of them constituting a complete orthonormal basis of the same space.

One can always build an unitary operator U that maps one basis to the other, i.e.: "| b[n] >=U | a[n] >"

*Verify that "U=(&sum;) | b[k] >< a[k] |" implies on  "| b[n] >=U | a[n] >"


*Show that "U=(&sum;) | b[k] > < a[k] | "is unitary


*Show that the matrix elements of U in the "| a[n] >" and  "| b[n] >" basis are equal


Show that A and `&Ascr;` = U*A*`#msup(mi("U"),mo("&dagger;"))`have the same spectrum



Schrödinger equation and unitary transform



Consider a ket "| psi[t] > " that solves the time-dependant Schrödinger equation:


"i `&hbar;` (&PartialD;)/(&PartialD;t) | psi[t] >=H(t) | psi[t] >"

and consider

"| phi[t] > =U(t) | psi[t] >",


where U(t) is a unitary operator.


Does "| phi[t] >" evolves according a Schrödinger equation

 "i*`&hbar;` (&PartialD;)/(&PartialD;t) | phi[t] >=`&Hscr;`(t) | phi[t] >"

and if yes, which is the expression of `&Hscr;`(t)?




Translation operators using Dirac notation


In this section, we focus on the operator T[a] = exp((-I*a*P)*(1/`&hbar;`))



The Action (translation) of the operator T[a]"=(e)^(-i (a P)/(`&hbar;`))" on a ket


Action of T[a] on an operatorV(X)


General Relativity


*Exact Solutions to Einstein's Equations  Lambda*g[mu, nu]+G[mu, nu] = 8*Pi*T[mu, nu]


*"Physical Review D" 87, 044053 (2013)


Given the spacetime metric,

g[mu, nu] = (Matrix(4, 4, {(1, 1) = -exp(lambda(r)), (1, 2) = 0, (1, 3) = 0, (1, 4) = 0, (2, 1) = 0, (2, 2) = -r^2, (2, 3) = 0, (2, 4) = 0, (3, 1) = 0, (3, 2) = 0, (3, 3) = -r^2*sin(theta)^2, (3, 4) = 0, (4, 1) = 0, (4, 2) = 0, (4, 3) = 0, (4, 4) = exp(nu(r))}))

a) Compute the Ricci and Weyl scalars


b) Compute the trace of


"Z[alpha]^(beta)=Phi R[alpha]^(beta)+`&Dscr;`[alpha]`&Dscr;`[]^(beta) Phi+T[alpha]^(beta)"


where `&equiv;`(Phi, Phi(r)) is some function of the radial coordinate, R[alpha, `~beta`] is the Ricci tensor, `&Dscr;`[alpha] is the covariant derivative operator and T[alpha, `~beta`] is the stress-energy tensor


T[alpha, beta] = (Matrix(4, 4, {(1, 1) = 8*exp(lambda(r))*Pi, (1, 2) = 0, (1, 3) = 0, (1, 4) = 0, (2, 1) = 0, (2, 2) = 8*r^2*Pi, (2, 3) = 0, (2, 4) = 0, (3, 1) = 0, (3, 2) = 0, (3, 3) = 8*r^2*sin(theta)^2*Pi, (3, 4) = 0, (4, 1) = 0, (4, 2) = 0, (4, 3) = 0, (4, 4) = 8*exp(nu(r))*Pi*epsilon}))

c) Compute the components of "W[alpha]^(beta)"" &equiv;"the traceless part of  "Z[alpha]^(beta)" of item b)


d) Compute an exact solution to the nonlinear system of differential equations conformed by the components of  "W[alpha]^(beta)" obtained in c)


Background: paper from February/2013, "Withholding Potentials, Absence of Ghosts and Relationship between Minimal Dilatonic Gravity and f(R) Theories", by P. Fiziev.


a) The Ricci and Weyl scalars


b) The trace of "  Z[alpha]^(beta)=Phi R[alpha]^(beta)+`&Dscr;`[alpha]`&Dscr;`[]^(beta) Phi+T[alpha]^(beta)"


b) The components of "W[alpha]^(beta)"" &equiv;"the traceless part of " Z[alpha]^(beta)"


c) An exact solution for the nonlinear system of differential equations conformed by the components of  "W[alpha]^(beta)"


*The Equivalence problem between two metrics



From the "What is new in Physics in Maple 2016" page:


In the Maple PDEtools package, you have the mathematical tools - including a complete symmetry approach - to work with the underlying [Einstein’s] partial differential equations. [By combining that functionality with the one in the Physics and Physics:-Tetrads package] you can also formulate and, depending on the metrics also resolve, the equivalence problem; that is: to answer whether or not, given two metrics, they can be obtained from each other by a transformation of coordinates, as well as compute the transformation.

Example from: A. Karlhede, "A Review of the Geometrical Equivalence of Metrics in General Relativity", General Relativity and Gravitation, Vol. 12, No. 9, 1980


*Equivalence for Schwarzschild metric (spherical and Krustal coordinates)


Tetrads and Weyl scalars in canonical form



Generally speaking a canonical form is obtained using transformations that leave invariant the tetrad metric in a tetrad system of references, so that theWeyl scalars are fixed as much as possible (conventionally, either equal to 0 or to 1).


Bringing a tetrad in canonical form is a relevant step in the tackling of the equivalence problem between two spacetime metrics.

The implementation is as in "General Relativity, an Einstein century survey", edited by S.W. Hawking (Cambridge) and W. Israel (U. Alberta, Canada), specifically Chapter 7 written by S. Chandrasekhar, page 388:








Residual invariance

Petrov type I







Petrov type II







Petrov type III







Petrov type D






`&Psi;__2`  remains invariant under rotations of Class III

Petrov type N






`&Psi;__4` remains invariant under rotations of Class II



The transformations (rotations of the tetrad system of references) used are of Class I, II and III as defined in Chandrasekar's chapter - equations (7.79) in page 384, (7.83) and (7.84) in page 385. Transformations of Class I can be performed with the command Physics:-Tetrads:-TransformTetrad using the optional argument nullrotationwithfixedl_, of Class II using nullrotationwithfixedn_ and of Class III by calling TransformTetrad(spatialrotationsm_mb_plan, boostsn_l_plane), so with the two optional arguments simultaneously.


The determination of appropriate transformation parameters to be used in these rotations, as well as the sequence of transformations happens all automatically by using the optional argument, canonicalform of TransformTetrad .


restart; with(Physics); with(Tetrads)

`Setting lowercaselatin letters to represent tetrad indices `


0, "%1 is not a command in the %2 package", Tetrads, Physics


0, "%1 is not a command in the %2 package", Tetrads, Physics


[IsTetrad, NullTetrad, OrthonormalTetrad, PetrovType, SimplifyTetrad, TransformTetrad, e_, eta_, gamma_, l_, lambda_, m_, mb_, n_]


Petrov type I


Petrov type II


Petrov type III


Petrov type N


Petrov type D


Edgardo S. Cheb-Terrab
Physics, Differential Equations and Mathematical Functions, Maplesoft

Hi all

I want to analysis this beam  (a simple and a tip rotational spring) in Maple to find the first four Frequencies...

Please help me

after following a example , got error


               1   / d      \    1        2     2
               - m |--- x(t)|  - - m omega  x(t)
               2   \ dt     /    2              
Error, (in Mechanics:-LagrangeEqs) invalid input: subs received subst1, which is not valid for its 1st argument
Error, invalid input: Mechanics:-GeneralSol expects its 1st argument, eqs, to be of type list, but received eqs
Error, invalid input: rhs received sol, which is not valid for its 1st argument, expr


Mechanics := module()
export SetVariables, LagrangeEqs, GeneralSol;
option package;
local subst1, subst2, varN, t;

SetVariables = proc( vars:: list, time )
local i;
t := time;
varN := nops( vars );
subst1 := {};
subst2 := {};
for i from 1 to var N do
subst1 := subst1 union
{vars[i](t) = q[i], diff(vars[i](t), t) = v[i]};
subst2 := subst2 union
{q[i] = vars[i](t), v[i] = diff(vars[i](t), t)};
end do;
print( subst1 );
print( subst2 );
end proc;

LagrangeEqs := proc (L)
local i, l1, term1, term2;
l1 := subs(subst1, L):
for i to varN do
term1 := [seq(diff(subs(subst2, diff(l1, v[i])), t), i = 1..varN)]:
term2 := [seq(subs(subst2, diff(l1, q[i])), i = 1..varN)]:
end do;
[ seq(simplify(term1[i]-term2[i]) = 0, i = 1..varN) ];
end proc;

RayleighEqs := proc(L, R)
local i, l1, r1, term1, term2, term3;
l1 := subs( subst1, L ):
r1 := subs( subst1, R ):
for i from 1 to varN do
term1:=[seq(diff(subs(subst2, diff(l1, v[i])), t), i=1..varN)]:
term2:=[seq(subs(subst2, diff(l1, q[i])), i=1..varN)]:
term3:=[seq(subs(subst2, diff(r1, v[i])), i=1..varN)]:
end do:
[ seq(simplify(term1[i]-term2[i]+term3[i]), i=1..varN) ];
end proc;

LagrEqsII := proc( L, Q::list )
local i, l1, term1, term2;
l1 := subs(subst1, L):
for i to varN do
term1 := [seq(diff(subs(subst2, diff(l1, v[i])), t), i = 1 .. varN)]:
term2 := [seq(subs(subst2, diff(l1, q[i])), i = 1 .. varN)]:
end do;
[seq(simplify(term1[i]-term2[i]) = Q[i], i = 1 .. varN)];
end proc;

LagrEqsIII := proc (L, R, Q::list)
local i, l1, r1, term1, term2, term3;
l1 := subs(subst1, L):
r1 := subs(subst1, R):
for i to varN do
term1 := [seq(diff(subs(subst2, diff(l1, v[i])), t), i = 1 .. varN)]:
term2 := [seq(subs(subst2, diff(l1, q[i])), i = 1 .. varN)]:
term3 := [seq(subs(subst2, diff(r1, v[i])), i = 1 .. varN)]:
end do;
[seq(simplify(term1[i]-term2[i]+term3[i]) = Q[i], i = 1 .. varN)];
end proc;

GeneralSol := proc (eqs::list)
local i, initconds, eqs2;
initconds := NULL:
eqs2 := eqs[][]:
for i to varN do
initconds:=VarNames[i](0)=q[i], (D(VarNames[i]))(0)=v[i], initconds:
end do;
dsolve({initconds, eqs2});
end proc;

end module;

LibLocation := cat("c:\\Temp");
Save(Mechanics, LibLocation);
save(Mechanics, "c:\\Temp\\Mechanics.m");
read "c:\\Temp\\Mechanics.m";


SetVariables([x], t);
L := (1/2)*m*diff(x(t), t)^2 - (1/2)*m*omega^2 * x(t)^2;
eqs := LagrangeEqs(L);
sol := GeneralSol( eqs );
X := unapply( rhs(sol), t );




      I am a student doing some self study over the summer trying to work through some of the John Taylor computer problems from his classcial mechanics book. Currently I hit a snag that most likely comes from the fact I am not well acquinted with Maple for solving IVP and DE's (we used Matlab in my DE class). I just need to know how I remove the following error:

Error, (in dsolve/numeric/SC/IVPsetup) initial conditions must be numeric

Here is a copy of my code:

R := 5;
g := 9.8;
deq1 := {diff(x(t), [`$`(t, 2)]) = -g*sin(x(t))/R, x(0) = 20};
/ d / d \ \
{ --- |--- x(t)| = -1.960000000 sin(x(t)), x(0) = 20 }
\ dt \ dt / /
dsol1 := dsolve(deq1, numeric);
Error, (in dsolve/numeric/SC/IVPsetup) initial conditions must be numeric

My hunch is that I need to set x'(0)=0 or something like I do not have enough intial values to solve the problem, but I could be wrong. Anyway anyone who can point out my mistake feel free to do so! Thank you!

Developed and then implemented with open code components. It is very important to note this post is held for students of civil engineering and mechanics. Using advanced mathematical concepts to concepts in engineering.

(in spanish)









How to building a wall, I want to use the wall with something (ex:ball) make a elastic collision.

Thank you all.

the calculation is like the following command, the result in the picture

SetCoordinates(spherical[r, theta, phi]);
Fv := rho*VectorField(`<,>`(v[r](r, theta, phi), v[theta](r, theta, phi), v[phi](r, theta, phi)));



1) when the Divergence act on the Fv, then it will be expanded, which is lengthy and not like most book's formulation , especially when I want to continue for a Conversation law like in fluid mechanics, this will be too long and a messy for later check.

could there be a way to not expand this result, just as the eq(3) like.

2) when I want to calculate the Divergence of Fv, I must construct a VectorField at first, but this is in components way, is there a quick way for Vector Field Function


(Presentation in Spain a month ago with a full description of the project and its current status)

A computational environment for Physicists



"Algebraic manipulations in Physics and related numerical exploration and visualization come together within computer algebra systems"

Project background


Three reasons for the underuse of Computer Algebra Systems in Physics


The Physics project goals


Status of things in Maple 17








Edgardo S. Cheb-Terrab
Physics, Maplesoft

Page 1 of 1