Advanced Search

Applications, Examples and Libraries

Share your work here
Create a new Post in this category

Calculating components of acceleration in Maple

by: Lenin Araujo Castillo 1790 Maple 18
engineering embedded-components normal
April 08 2018
2  2

Using Maple's native syntax, we can calculate the components of acceleration. That is, the tangent and normal scalar component with its respective units of measure. Now the difficult calculations were in the past because with Maple we solved it and we concentrated on the interpretation of the results for engineering. In spanish.

Calculo_Componentes_Aceleracion_Curvilínea.mw

Uso_de_comandos_y_operadores_para_calculos_de_componentes_de_la_aceleración.mw

Lenin Araujo Castillo

Ambassador of Maple

 

 

Maple 2018: Updates to Physics Differential Equations...

by: ecterrab 13431 Maple
maplecloud
March 28 2018
7  0

The first update to the Maple 2018 Physics, Differential Equations and Mathematical Functions packages is available. As has been the case since 2013, this update contains fixes, enhancements to existing functionality, and new developments in the three areas. 

The webpage for these updates will continue being the Maplesoft R&D Physics webpage. Starting with Maple 2018, however, this update is also available from the MapleCloud.

To install the update: open Maple and click the Cloud icon (upper-right corner), select "Packages" and search for "Physics Updates". Then, in the corresponding "Actions" column, click the third icon (install pop-up).

NOTE May/1: the "Updates" icon of the MapleCloud toolbar (that opens when you click the upper-right icon within a Maple document / worksheet), now works fine, after having installed the Physics Updates version 32 or higher.

These first updates include:

  • New Physics functionality regarding Tensor Products of Quantum States; and Coherent States.
  • Updates to pdsolve regarding PDE & Boundary Conditions (exact solutions);
  • A change in notation: d_(x), the differential of a coordinate in the Physics package, is now displayed as shown in this Mapleprimes post.


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

Tensor Product of Quantum State Spaces

by: ecterrab 13431 Maple
quantum
March 23 2018
4  0

This is about the recent implementation of tensor products of quantum state spaces in the Physics package, in connection with an exchange with the Physics of Information Lab of the University of Waterloo. As usual this development is available to everybody from the Maplesoft R&D Physics webpage. This is the last update for Maple 2017. The updates for Maple 2018, starting with this same material, will begin being distributed through the MapleCloud next week.

Tensor Product of Quantum State Spaces

 

Basic ideas and design

 

 

Suppose A and B are quantum operators and Ket(A, n), et(B, m) are, respectively, their eigenkets. The following works since the introduction of the Physics package in Maple

> 

with(Physics)

> 

Setup(op = {A, B})

`* Partial match of 'op' against keyword 'quantumoperators'`

 

[quantumoperators = {A, B}]

(1)
> 

A*Ket(A, alpha) = A.Ket(A, alpha)

Physics:-`*`(A, Physics:-Ket(A, alpha)) = alpha*Physics:-Ket(A, alpha)

(2)
> 

B*Ket(B, beta) = B.Ket(B, beta)

Physics:-`*`(B, Physics:-Ket(B, beta)) = beta*Physics:-Ket(B, beta)

(3)

where on the left-hand sides the product operator `*` is used as a sort of inert form (it has all the correct mathematical properties but does not perform the contraction) of the dot product operator `.`, used on the right-hand sides.

 

Suppose now that A and B act on different, disjointed, Hilbert spaces.

 

1) To represent that, a new keyword in Setup , is introduced, to indicate which spaces are disjointed, say as in disjointedhilbertspaces = {A, B}.  We want this notation to pop up at some point as {`ℋ`[A], `ℋ`[B]} where the indexation indicates all the operators acting on that Hilbert space. The disjointedspaces keyword has for synonyms disjointedhilbertspaces and hilbertspaces. The display `ℋ`[A] is not yet implemented.

 

NOTE: noncommutative quantum operators acting on disjointed spaces commute between themselves, so after setting  - for instance - disjointedspaces = {A, B}, automatically, A, B become quantum operators satisfying (see comment (ii) on page 156 of ref. [1])

 

"[A,B][-]=0"

 

2) Product of Kets and Bras (KK, BB, KB and BK) where K and B belong to disjointed spaces, are understood as tensor products satisfying, for instance with disjointed spaces A and B (see footnote on page 154 of ref. [1]),

 

`⊗`(Ket(A, alpha), Ket(B, beta)) = `⊗`(Ket(B, beta), Ket(A, alpha)) 

 

`⊗`(Bra(A, alpha), Ket(B, beta)) = `⊗`(Ket(B, beta), Bra(A, alpha)) 

 

while of course

Bra(A, alpha)*Ket(A, alpha) <> Bra(A, alpha)*Ket(A, alpha)

 

Details

   

 

3) All the operators of one disjointed space act transparently over operators, Bras and Kets of the other disjointed spaces, for example

 

A*Ket(B, n) = A*Ket(B, n)

and the same for the Dagger of this equation, that is

Bra(B, n)*Dagger(A) = Bra(B, n)*Dagger(A)

 

And this happens automatically. Hence, when we write the left-hand sides and press enter, they are automatically rewritten (returned) as the right-hand sides.

 

Note that for the product of an operator times a Bra or a Ket we are not using the notation that expresses the product with the symbol 5.

 

Regarding the display of Bras and Kets and their tensor products, two enhancements are happening:

 

• 

A new Setup option hideketlabel makes all the labels in Kets and Bras to be hidden when displaying Kets, Bras and Bracket(s), with the indices presented one level up, as if they were a sequence of labels, so that:

 "Ket(A,m,n,l"  

is displayed as

Ket(A, m, n, l)

 

  

This is the notation used more frequently when working in quantum information. This hideketlabel option is already implemented entering Setup(hideketlabel = true)
• 

Tensor products formed with operators, or Bras and Kets, that belong to disjointed spaces (set as such using Setup ), are displayed with the symbol 5 in between, as in Ket(A, n)*Ket(B, n) instead of Ket(A, n)*Ket(B, n), and `&otimes;`(A, B) instead of A*B.

Tensor product notation and the hideketlabel option

   

The implementation of tensor products using `*` and `.`

   

Basic exercising with the new functionality

   

Related functionality already in place before these changes

   

Reference

 

[1] Cohen-Tannoudji, Diue, Laloe, Quantum Mechanics, Chapter 2, section F.

[2] Griffiths Robert B., Hilbert Space Quantum Mechanics, Quantum Computation and Quantum Information Theory Course, Physics Department, Carnegie Mellon University, 2014.

See also

   

 


 

Download TensorProductDesign.mw
 

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

Carmichael's lambda (n) function (reduced Euler...

by: kenrubenstein 30
animation totient
March 22 2018
1  0

Carmichael's lambda(n) function (as it relates to Euler's Totient Function).....this is just one of 8 stages of animation. 

 

Here's the complete animation with supporting music by the mighty Stormtroopers Of Death....

https://youtu.be/QN-s3EpZICs

 

 

 

 

A Boggle board

by: Christopher2222 5785 Maple
imagetools boggle rotate
March 16 2018
1  2

I've created a worksheet that outputs a boggle board.  I think it could be more efficient than the method I came up with but the idea is there.  The only way I could figure to rotate the letters was to output them to a bmp format then read them back in and use imagetools for rotation.  I used Times Roman font but the font Boggle uses I think is Tunga, Latha or Mangal.  

Note - remove the colon in the last line to produce the output.  One other thing I believe, in Tools->Options-> (uncheck)Limit Expression Length to 1000000  

Saving the file with the output would have produced a file in the tens of Megabytes and may have caused error loading.


 

restart; gc()

with(plots); with(ImageTools)

a := [seq(k, k = "A" .. "Z")]

["A", "B", "C", "D", "E", "F", "G", "H", "I", "J", "K", "L", "M", "N", "O", "P", "Q", "R", "S", "T", "U", "V", "W", "X", "Y", "Z"]

 (1)

for i in a do plotsetup(bmp, plotoutput = cat("c:/test/", i, ".bmp")); img || i := textplot([0, 1, i, font = ["times", "roman", 200]], axes = none, scaling = constrained) end do
NULL

plotsetup = default

 (2)

plotsetup(default)

with(combinat)``

``

Setting up the 16 boggle cubes

 

cube1 := ["H", "E", "E", "N", "W", "G"]

cube2 := ["T", "M", "I", "O", "C", "U"]

cube3 := ["D", "E", "X", "L", "R", "I"]

cube4 := ["S", "P", "F", "A", "K", "F"]

cube5 := ["T", "O", "E", "S", "I", "S"]

cube6 := ["H", "N", "L", "N", "Z", "R"]

cube7 := ["R", "L", "T", "Y", "T", "E"]

cube8 := ["D", "E", "Y", "L", "R", "V"]

cube9 := ["C", "A", "O", "S", "P", "H"]

cube10 := ["Qu", "U", "M", "H", "I", "N"]

cube11 := ["D", "Y", "I", "S", "T", "T"]

cube12 := ["S", "N", "I", "E", "E", "U"]

cube13 := ["T", "O", "O", "W", "A", "T"]

cube14 := ["W", "H", "E", "V", "R", "T"]

cube15 := ["J", "B", "O", "O", "A", "B"]

cube16 := ["N", "A", "E", "A", "E", "G"]

cubes := [seq(cat("cube", i), i = 1 .. 16)]

["cube1", "cube2", "cube3", "cube4", "cube5", "cube6", "cube7", "cube8", "cube9", "cube10", "cube11", "cube12", "cube13", "cube14", "cube15", "cube16"]

 (3)

c := randperm(cubes)

["cube13", "cube14", "cube11", "cube6", "cube9", "cube1", "cube16", "cube7", "cube3", "cube2", "cube12", "cube5", "cube4", "cube10", "cube8", "cube15"]

 (4)

cc := map(parse, c)

[cube13, cube14, cube11, cube6, cube9, cube1, cube16, cube7, cube3, cube2, cube12, cube5, cube4, cube10, cube8, cube15]

 (5)

ccf := [seq(op(randcomb(cc[i], 1)), i = 1 .. 16)]

["A", "H", "I", "N", "A", "N", "E", "T", "D", "M", "E", "S", "K", "N", "E", "O"]

 (6)

with(ArrayTools)

g := Reshape(Array(ccf), [4, 4])

Array(%id = 18446744074360417206)

 (7)

rr := proc () randcomb([0, 90, 180, 270], 1) end proc

Reshape(Array([seq(display(Preview(Rotate(Read(cat("c:/test/", ccf[i], ".bmp")), op(rr()))), axes = none), i = 1 .. 16)]), [4, 4])
 

````

 

 

 

NULL


 

Download Boggle3-6final.mw

Quantum Teleportation in Dirac notation with...

by: Louis Lamarche 105 Maple
physics tensor quantum-mechanics
March 08 2018
4  4

This is an application of the previous posts
https://www.mapleprimes.com/posts/209057-Procedure-For-Expanding-Tensor-Product

I have a fourth version of the ExpandQop that will expand automaticaly the power of
quantum tensor product. This is just a minor change to the procedure.

Now here is an application for all this that will help understanding a little about
quantum computing. This is the classical concept of quantum teleportation.

You will need to run the above mentionned file and uncomment the save line in the file
before running the example.

LL
 

> 

######################################################################
# NOTICE                                                             #
# Author: Louis Lamarche                                             #
#         Institute of Research of Hydro-Quebec (IREQ)               #
#         Science des données et haute performance                   #
#         2018, March 7                                              #
#                                                                    #
# Function name: ExpandQop (x)                                       #
#       Purpose: Compute the tensor product of two quantum           #
#                operators in Dirac notations                        #
#      Argument: x: a quantum operator                               #
#  Improvements: Manage all +, -, *, /, ^, mod  operations           #
#                in the argument. Manages multiple tensor products   #
#                like A*B*C*F                                        #
#       Version: 3                                                   #
#                                                                    #
#  Copyrigth(c) Hydro-Quebec.                                        #
#        Note 1: Permission to use this softwate is granted if you   #
#                acknowledge its author and copyright                #
#        Note 2: Permission to copy this softwate is granted if you  #
#                leave this 21 lines notice intact. Thank you.       #
######################################################################
restart;
> 

with(Physics):
interface(imaginaryunit=i):
Setup(mathematicalnotation=true);

[mathematicalnotation = true]

 (1)
> 

Setup(unitaryoperators={I,U,X,Y,Z,H,HI,CNOT,CnotI});
Setup(noncommutativeprefix={q,beta,psi});

[unitaryoperators = {CNOT, CnotI, H, HI, I, U, X, Y, Z}]

  

[noncommutativeprefix = {beta, psi, q}]

 (2)
> 

Setup(bracketrules= { %Bracket(%Bra(q0), %Ket(q0))=1,
                      %Bracket(%Bra(q1), %Ket(q1))=1,
                      %Bracket(%Bra(q1), %Ket(q0))=0,
                      %Bracket(%Bra(q0), %Ket(q1))=0
                    });

[bracketrules = {%Bracket(%Bra(q0), %Ket(q0)) = 1, %Bracket(%Bra(q0), %Ket(q1)) = 0, %Bracket(%Bra(q1), %Ket(q0)) = 0, %Bracket(%Bra(q1), %Ket(q1)) = 1}]

 (3)
> 

####################################################################################
# Load the procedure and set the required global variables
#
read "ExpandQop.m": optp:=op(0,Ket(q0)*Ket(q1)): optpx:= op(0,(Ket(q0)+Ket(q1))^2):
#
####################################################################################
> 

#
# Pauli operators
#
print("Pauli gates");
I:=Ket(q0)*Bra(q0)+Ket(q1)*Bra(q1);        # = sigma[0]
X:=Ket(q1)*Bra(q0)+Ket(q0)*Bra(q1);        # = sigma[1] = sigma[x]
Y:=-i*Ket(q1)*Bra(q0)+i*Ket(q0)*Bra(q1);   # = sigma[2] = sigma[y]
Z:=Ket(q0)*Bra(q0)-Ket(q1)*Bra(q1);        # = sigma[3] = sigma[z]

"Pauli gates"

  

Physics:-`*`(Physics:-Ket(q0), Physics:-Bra(q0))+Physics:-`*`(Physics:-Ket(q1), Physics:-Bra(q1))

  

Physics:-`*`(Physics:-Ket(q1), Physics:-Bra(q0))+Physics:-`*`(Physics:-Ket(q0), Physics:-Bra(q1))

  

-I*Physics:-`*`(Physics:-Ket(q1), Physics:-Bra(q0))+I*Physics:-`*`(Physics:-Ket(q0), Physics:-Bra(q1))

  

Physics:-`*`(Physics:-Ket(q0), Physics:-Bra(q0))-Physics:-`*`(Physics:-Ket(q1), Physics:-Bra(q1))

 (4)
> 

##############################
# Defining the Hadamard gate #
##############################
print("Hadamard gate");
H:= Ket(q0)*Bra(q0)/sqrt(2)+Ket(q0)*Bra(q1)/sqrt(2)+Ket(q1)*Bra(q0)/sqrt(2)-Ket(q1)*Bra(q1)/sqrt(2);

"Hadamard gate"

  

(1/2)*2^(1/2)*Physics:-`*`(Physics:-Ket(q0), Physics:-Bra(q0))+(1/2)*2^(1/2)*Physics:-`*`(Physics:-Ket(q0), Physics:-Bra(q1))+(1/2)*2^(1/2)*Physics:-`*`(Physics:-Ket(q1), Physics:-Bra(q0))-(1/2)*2^(1/2)*Physics:-`*`(Physics:-Ket(q1), Physics:-Bra(q1))

 (5)
> 

# This is usefull to represent a 2 qubits system
# A more general approach is needed for a n qubit system.
DefineStates:=proc()
    Ket(q00):=Ket(q0)*Ket(q0);  Ket(q01):=Ket(q0)*Ket(q1);
    Ket(q10):=Ket(q1)*Ket(q0);  Ket(q11):=Ket(q1)*Ket(q1);
    Bra(q00):=Dagger(Ket(q00)); Bra(q01):=Dagger(Ket(q01));
    Bra(q10):=Dagger(Ket(q10)); Bra(q11):=Dagger(Ket(q11));
    return;
    end proc:
UndefineStates:=proc()
    Ket(q00):='Ket(q00)'; Ket(q01):='Ket(q01)';
    Ket(q10):='Ket(q10)'; Ket(q11):='Ket(q11)';
    Bra(q00):='Bra(q00)'; Bra(q01):='Bra(q01)';
    Bra(q10):='Bra(q10)'; Bra(q11):='Bra(q11)';
    return;
    end proc:
> 

####################################
# Defining the CNOT gate (2 qubits)
####################################
print("CNOT gate");
CNOT:=Ket(q00)*Bra(q00)+ Ket(q01)*Bra(q01)+ Ket(q11)*Bra(q10)+Ket(q10)*Bra(q11);
DefineStates();
'CNOT'=CNOT;

"CNOT gate"

  

Physics:-`*`(Physics:-Ket(q00), Physics:-Bra(q00))+Physics:-`*`(Physics:-Ket(q01), Physics:-Bra(q01))+Physics:-`*`(Physics:-Ket(q11), Physics:-Bra(q10))+Physics:-`*`(Physics:-Ket(q10), Physics:-Bra(q11))

  

CNOT = Physics:-`*`(Physics:-Ket(q0), Physics:-Ket(q0), Physics:-Bra(q0), Physics:-Bra(q0))+Physics:-`*`(Physics:-Ket(q0), Physics:-Ket(q1), Physics:-Bra(q1), Physics:-Bra(q0))+Physics:-`*`(Physics:-Ket(q1), Physics:-Ket(q1), Physics:-Bra(q0), Physics:-Bra(q1))+Physics:-`*`(Physics:-Ket(q1), Physics:-Ket(q0), Physics:-Bra(q1), Physics:-Bra(q1))

 (6)
> 

###########################
# Defining the Bell states
###########################
Ket(beta,x,y)='CNOT.(((H.Ket(x)))*Ket(y))';
Ket(beta00):=CNOT.(Expand((H.Ket(q0)))*Ket(q0));
Ket(beta01):=CNOT.(Expand((H.Ket(q0)))*Ket(q1));
Ket(beta10):=CNOT.(Expand((H.Ket(q1)))*Ket(q0));
Ket(beta11):=CNOT.(Expand((H.Ket(q1)))*Ket(q1));

Physics:-Ket(beta, x, y) = Physics:-`.`(CNOT, Physics:-`*`(Physics:-`.`(H, Physics:-Ket(x)), Physics:-Ket(y)))

  

(1/2)*2^(1/2)*(Physics:-`*`(Physics:-Ket(q0), Physics:-Ket(q0))+Physics:-`*`(Physics:-Ket(q1), Physics:-Ket(q1)))

  

(1/2)*2^(1/2)*(Physics:-`*`(Physics:-Ket(q0), Physics:-Ket(q1))+Physics:-`*`(Physics:-Ket(q1), Physics:-Ket(q0)))

  

-(1/2)*2^(1/2)*(-Physics:-`*`(Physics:-Ket(q0), Physics:-Ket(q0))+Physics:-`*`(Physics:-Ket(q1), Physics:-Ket(q1)))

  

(1/2)*2^(1/2)*(Physics:-`*`(Physics:-Ket(q0), Physics:-Ket(q1))-Physics:-`*`(Physics:-Ket(q1), Physics:-Ket(q0)))

 (7)
> 

##########################################################
# Quantum teleportation
# Reference: Quantum Computation and Quantum Information
#            10th Anniversary Edition
#            Michael A. Nielsen & Isaac L. Chuang
#            Cambridge University Press, Cambridge 2010
#            pp 25-28
##########################################################
print("State to be teleported");
Ket(psi) := a*Ket(q0)+b*Ket(q1);
print("Step 1: Compute the tensor product of the state to be teleported with ", 'Ket(beta00)');
Ket(psi[0])='Ket(psi)'*'Ket(beta00)';
Ket(psi[0]):=Expand(Ket(psi)*Ket(beta00));
print("This is a 3 qubits state");
#######
print("Step 2: Pass these 3 qubits through a  CNOT*I  operator");
'CnotI'='CNOT*I';
CnotI:=ExpandQop(Expand(CNOT*I)):
#
# To see what the CNOTI operator looks like
#
# print("CNOTI=");
# print(op(1,CNOTI)+op(2,CNOTI)+op(3,CNOTI)+op(4,CNOTI));
# print(op(5,CNOTI)+op(6,CNOTI)+op(7,CNOTI)+op(8,CNOTI));
'Ket(psi[1])'='CnotI.Ket(psi[0])';
Ket(psi[1]):=Expand(CnotI.Ket(psi[0]));
#######
print("Step 3: Pass these 3 qubits through an Haldamard*I  operator");
'HalI'='H*I';
HalI:=ExpandQop(Expand(H*I)):
#
# To see what the Haldamard*I operator looks like
#
# print("HalI=");
# print(op(1,HalI)+op(2,HalI)+op(3,HalI)+op(4,HalI));
# print(op(5,HalI)+op(6,HalI)+op(7,HalI)+op(8,HalI));
'Ket(psi[2])'='HalI.Ket(psi[1])';
Ket(psi[2]):=Expand(HalI.Ket(psi[1]));
 

"State to be teleported"

  

a*Physics:-Ket(q0)+b*Physics:-Ket(q1)

  

"Step 1: Compute the tensor product of the state to be teleported with ", Physics:-Ket(beta00)

  

Physics:-Ket(psi[0]) = Physics:-`*`(Physics:-Ket(psi), Physics:-Ket(beta00))

  

(1/2)*2^(1/2)*a*Physics:-`*`(Physics:-Ket(q0), Physics:-Ket(q0), Physics:-Ket(q0))+(1/2)*2^(1/2)*a*Physics:-`*`(Physics:-Ket(q0), Physics:-Ket(q1), Physics:-Ket(q1))+(1/2)*2^(1/2)*b*Physics:-`*`(Physics:-Ket(q1), Physics:-Ket(q0), Physics:-Ket(q0))+(1/2)*2^(1/2)*b*Physics:-`*`(Physics:-Ket(q1), Physics:-Ket(q1), Physics:-Ket(q1))

  

"This is a 3 qubits state"

  

"Step 2: Pass these 3 qubits through a CNOT*I operator"

  

CnotI = Physics:-`*`(CNOT, I)

  

Physics:-Ket(psi[1]) = Physics:-`.`(CnotI, Physics:-Ket(psi[0]))

  

(1/2)*2^(1/2)*a*Physics:-`*`(Physics:-Ket(q0), Physics:-Ket(q0), Physics:-Ket(q0))+(1/2)*2^(1/2)*a*Physics:-`*`(Physics:-Ket(q0), Physics:-Ket(q1), Physics:-Ket(q1))+(1/2)*2^(1/2)*b*Physics:-`*`(Physics:-Ket(q1), Physics:-Ket(q1), Physics:-Ket(q0))+(1/2)*2^(1/2)*b*Physics:-`*`(Physics:-Ket(q1), Physics:-Ket(q0), Physics:-Ket(q1))

  

"Step 3: Pass these 3 qubits through an Haldamard*I operator"

  

HalI = Physics:-`*`(H, I)

  

Physics:-Ket(psi[2]) = Physics:-`.`(HalI, Physics:-Ket(psi[1]))

  

(1/2)*a*Physics:-`*`(Physics:-Ket(q0), Physics:-Ket(q0), Physics:-Ket(q0))+(1/2)*a*Physics:-`*`(Physics:-Ket(q0), Physics:-Ket(q1), Physics:-Ket(q1))+(1/2)*b*Physics:-`*`(Physics:-Ket(q0), Physics:-Ket(q0), Physics:-Ket(q1))+(1/2)*b*Physics:-`*`(Physics:-Ket(q0), Physics:-Ket(q1), Physics:-Ket(q0))+(1/2)*a*Physics:-`*`(Physics:-Ket(q1), Physics:-Ket(q0), Physics:-Ket(q0))+(1/2)*a*Physics:-`*`(Physics:-Ket(q1), Physics:-Ket(q1), Physics:-Ket(q1))-(1/2)*b*Physics:-`*`(Physics:-Ket(q1), Physics:-Ket(q0), Physics:-Ket(q1))-(1/2)*b*Physics:-`*`(Physics:-Ket(q1), Physics:-Ket(q1), Physics:-Ket(q0))

 (8)
> 

UndefineStates();
print("Using contracted names for the first two qubits");
Ket(q00)*Bra(q0)*Bra(q0)='I';
Ket(q01)*Bra(q0)*Bra(q1)='I';
Ket(q10)*Bra(q1)*Bra(q0)='I';
Ket(q11)*Bra(q1)*Bra(q1)='I';
'Ket(psi[2])'=Ket(q00)*Bra(q0)*Bra(q0).Ket(psi[2])+
              Ket(q01)*Bra(q0)*Bra(q1).Ket(psi[2])+
              Ket(q10)*Bra(q1)*Bra(q0).Ket(psi[2])+
              Ket(q11)*Bra(q1)*Bra(q1).Ket(psi[2]);

"Using contracted names for the first two qubits"

  

Physics:-`*`(Physics:-Ket(q00), Physics:-Bra(q0), Physics:-Bra(q0)) = I

  

Physics:-`*`(Physics:-Ket(q01), Physics:-Bra(q0), Physics:-Bra(q1)) = I

  

Physics:-`*`(Physics:-Ket(q10), Physics:-Bra(q1), Physics:-Bra(q0)) = I

  

Physics:-`*`(Physics:-Ket(q11), Physics:-Bra(q1), Physics:-Bra(q1)) = I

  

Physics:-Ket(psi[2]) = (1/2)*a*Physics:-`*`(Physics:-Ket(q00), Physics:-Ket(q0))+(1/2)*b*Physics:-`*`(Physics:-Ket(q00), Physics:-Ket(q1))+(1/2)*a*Physics:-`*`(Physics:-Ket(q01), Physics:-Ket(q0))-(1/2)*b*Physics:-`*`(Physics:-Ket(q01), Physics:-Ket(q1))+(1/2)*a*Physics:-`*`(Physics:-Ket(q10), Physics:-Ket(q1))+(1/2)*b*Physics:-`*`(Physics:-Ket(q10), Physics:-Ket(q0))+(1/2)*a*Physics:-`*`(Physics:-Ket(q11), Physics:-Ket(q1))-(1/2)*b*Physics:-`*`(Physics:-Ket(q11), Physics:-Ket(q0))

 (9)
> 

print("Rewriting this result by hand");
'Ket(psi[2])'=(Ket(q00)*(a*Ket(q0)+b*Ket(q1))+
               Ket(q01)*(a*Ket(q0)-b*Ket(q1))+
               Ket(q10)*(a*Ket(q1)+b*Ket(q0))+
               Ket(q11)*(a*Ket(q1)-b*Ket(q0)))/2;

"Rewriting this result by hand"

  

Physics:-Ket(psi[2]) = (1/2)*Physics:-`*`(Physics:-Ket(q00), a*Physics:-Ket(q0)+b*Physics:-Ket(q1))+(1/2)*Physics:-`*`(Physics:-Ket(q01), a*Physics:-Ket(q0)-b*Physics:-Ket(q1))+(1/2)*Physics:-`*`(Physics:-Ket(q10), a*Physics:-Ket(q1)+b*Physics:-Ket(q0))+(1/2)*Physics:-`*`(Physics:-Ket(q11), a*Physics:-Ket(q1)-b*Physics:-Ket(q0))

 (10)
> 

DefineStates();
print("If Alice measures 00 Bob does noting");
''I'.   '2*Bra(q00).Ket(psi[2])'' =  I.   2*Bra(q00).Ket(psi[2]);
print("If Alice measures 01 Bob applies the X gate");
''X'.   '2*Bra(q01).Ket(psi[2])'' =  X.   2*Bra(q01).Ket(psi[2]);
print("If Alice measures 10 Bob applies the Z gate");
''Z'.   '2*Bra(q10).Ket(psi[2])'' =  Z.   2*Bra(q10).Ket(psi[2]);
print("If Alice measures 11 Bob applies the X gate and then the Z gate");
''Z'.'X'. '2*Bra(q11).Ket(psi[2])'' =  Z.X. 2*Bra(q11).Ket(psi[2]);

"If Alice measures 00 Bob does noting"

  

Physics:-`.`('I', 'Physics:-`.`(Physics:-`*`(2, Physics:-Bra(q00)), Physics:-Ket(psi[2]))') = a*Physics:-Ket(q0)+b*Physics:-Ket(q1)

  

"If Alice measures 01 Bob applies the X gate"

  

Physics:-`.`('X', 'Physics:-`.`(Physics:-`*`(2, Physics:-Bra(q01)), Physics:-Ket(psi[2]))') = a*Physics:-Ket(q0)+b*Physics:-Ket(q1)

  

"If Alice measures 10 Bob applies the Z gate"

  

Physics:-`.`('Z', 'Physics:-`.`(Physics:-`*`(2, Physics:-Bra(q10)), Physics:-Ket(psi[2]))') = a*Physics:-Ket(q0)+b*Physics:-Ket(q1)

  

"If Alice measures 11 Bob applies the X gate and then the Z gate"

  

Physics:-`.`('Z', 'X', 'Physics:-`.`(Physics:-`*`(2, Physics:-Bra(q11)), Physics:-Ket(psi[2]))') = a*Physics:-Ket(q0)+b*Physics:-Ket(q1)

 (11)
> 

 


 

Download QuantumTeleportation.mw

 

 

Procedure for expanding tensor product of quantum...

by: Louis Lamarche 105 Maple
physics tensor quantum-mechanics
March 08 2018
1  1

Version 2 do not enable to expand multiple product like A*A*B*E
Version 3 will now do that.
I just forgot to add this feature.

LL.
 

> 

######################################################################
# NOTICE                                                             #
# Author: Louis Lamarche                                             #
#         Institute of Research of Hydro-Quebec (IREQ)               #
#         Science des données et haute performance                   #
#         2018, March 7                                              #
#                                                                    #
# Function name: ExpandQop (x)                                       #
#       Purpose: Compute the tensor product of two quantum           #
#                operators in Dirac notations                        #
#      Argument: x: a quantum operator                               #
#  Improvements: Manage all +, -, *, /, ^, mod  operations           #
#                in the argument. Manages multiple tensor products   #
#                like A*B*C*F                                        #
#       Version: 3                                                   #
#                                                                    #
#  Copyrigth(c) Hydro-Quebec.                                        #
#        Note 1: Permission to use this softwate is granted if you   #
#                acknowledge its author and copyright                #
#        Note 2: Permission to copy this softwate is granted if you  #
#                leave this 21 lines notice intact. Thank you.       #
######################################################################
restart;
> 

with(Physics):
interface(imaginaryunit=i):
Setup(mathematicalnotation=true);

[mathematicalnotation = true]

 (1)
> 

Setup(quantumoperators={A,B,C,Cn});
Setup(noncommutativeprefix={a,b,q});

[quantumoperators = {A, B, C, Cn}]

  

[noncommutativeprefix = {a, b, q}]

 (2)
> 

opexp:= op(0,10^x):            # exponentiation id
opnp := op(0,10*x):            # normal product id
optp := op(0,Ket(q0)*Ket(q1)): # tensor product id
opdiv:= `Fraction`:            # fraction       id          
opsym:= op(0,x):               # symbol         id
opint:= op(0,10):              # integer        id
opflt:= op(0,10.0):            # float          id
opcpx:= op(0,i):               # complex        id
opbra:= op(0,Bra(q)):          # bra            id
opket:= op(0,Ket(q)):          # ket            id
opmod:= op(0, a mod b):        # mod            id
ExpandQop:=proc(x)
    local nx,ret,j,lkb,cbk,rkb,no,lop,success;
    lop:=op(0,x);
    no:=nops(x);
    if lop = opsym or lop = opint or lop = opflt or
       lop = opbra or lop = opket or lop = opcpx then
         ret:=x;
    else
    if lop = opexp then
        ret:=x;
    else       
    if lop = opnp then
        ret:=1;
        for j from 1 to no do
            ret:=ret*ExpandQop(op(j,x));
        end do;        
    else
    if lop = `+` then
        ret:=0;
        for j from 1 to no do
            ret:=ret+ExpandQop(op(j,x));
        end do;
    else
    if lop = `-` then
        ret:=0;
        for j from 1 to no do
            ret:=ret-ExpandQop(op(j,x));
        end do;
    else
    if lop = opdiv then
       ret:=1;
       for j from 1 to no do
           ret:=ret/ExpandQop(op(j,x));
       end do;
    else
    if lop = opmod then
       ret:=x;
    else
    if lop = optp then
       if (no > 3 ) then
           success:=false;
           nx:=x;
           while not success do
             lkb:=0; cbk:=0; rkb:=0;ret:=1;
             for j from 1 to no do
                 if (j>1) then
                      if(lkb=0) then
                          if( type(op(j-1,nx),Ket) and
                              type(op(j,nx),Bra) ) then lkb:=j-1; fi;
                      else
                          if( type(op(j-1,nx),Ket) and
                              type(op(j,nx),Bra) ) then rkb:=j;   fi;
                      fi;
                      if( type(op(j-1,nx),Bra) and type(op(j,nx),Ket) )
                                                   then cbk:=j;   fi;
                 fi;
             end do;
             if ( (lkb < cbk) and (cbk<rkb) ) then
                 for j from 1     to lkb   do ret := ret*op(j,nx); end do;
                 for j from cbk   to no    do ret := ret*op(j,nx); end do;
                 for j from lkb+1 to cbk-1 do ret := ret*op(j,nx); end do;
             else
               ret:=nx;
             fi;
             
             if nx = ret then
                success := true;
             else
                nx := ret;
             fi
           end do;
       else
           ret:=x;
       fi;
    else ret:=x;
    fi; # optp
    fi; # opmod
    fi; # opdiv
    fi; # `-`
    fi; # `+`
    fi; # `opnp
    fi; # `opexp`
    fi; # opsym, opint, opflt, opbra, opket, opcpx

    return ret;
end proc:

# For saving
# save opexp,opnp,optp,opdiv,opint,opflt,opcpx,opbra,opket,opmod, ExpandQop,"ExpandQop.m"
> 

# Let A be an operator in a first Hilbert space of dimension n
#  using the associated orthonormal ket and bra vectors
#
#
kets1:=Ket(a1)*Ket(a2)*Ket(a3)*Ket(a4)*Ket(a5):
A:=kets1*Dagger(kets1);


# Let B be an operator in a second Hilbert (Ketspace of dimension m
#  using the associated orthonormal ket and bra vectors
#
#
kets2:=Ket(b1)*Ket(b2)*Ket(b3):
B:=kets2*Dagger(kets2);


# The tensor product of the two operators acts on a n+m third
# Hilbert space   unsing the appropriately ordered ket
# and bra  vectors of the two preceding spaces. The rule for
# building this operator in Dirac notation is as follows,
#
#


print("Maple do not compute the tensor product of operators,");
print("C=A*B gives:");
C:=A*B;

print("ExpandQop(C) gives the expected result:");
Cn:=ExpandQop(C);

Physics:-`*`(Physics:-Ket(a1), Physics:-Ket(a2), Physics:-Ket(a3), Physics:-Ket(a4), Physics:-Ket(a5), Physics:-Bra(a5), Physics:-Bra(a4), Physics:-Bra(a3), Physics:-Bra(a2), Physics:-Bra(a1))

  

Physics:-`*`(Physics:-Ket(b1), Physics:-Ket(b2), Physics:-Ket(b3), Physics:-Bra(b3), Physics:-Bra(b2), Physics:-Bra(b1))

  

"Maple do not compute the tensor product of operators,"

  

"C=A*B gives:"

  

Physics:-`*`(Physics:-Ket(a1), Physics:-Ket(a2), Physics:-Ket(a3), Physics:-Ket(a4), Physics:-Ket(a5), Physics:-Bra(a5), Physics:-Bra(a4), Physics:-Bra(a3), Physics:-Bra(a2), Physics:-Bra(a1), Physics:-Ket(b1), Physics:-Ket(b2), Physics:-Ket(b3), Physics:-Bra(b3), Physics:-Bra(b2), Physics:-Bra(b1))

  

"ExpandQop(C) gives the expected result:"

  

Physics:-`*`(Physics:-Ket(a1), Physics:-Ket(a2), Physics:-Ket(a3), Physics:-Ket(a4), Physics:-Ket(a5), Physics:-Ket(b1), Physics:-Ket(b2), Physics:-Ket(b3), Physics:-Bra(b3), Physics:-Bra(b2), Physics:-Bra(b1), Physics:-Bra(a5), Physics:-Bra(a4), Physics:-Bra(a3), Physics:-Bra(a2), Physics:-Bra(a1))

 (3)
> 

kets3:=kets1*kets2;
bras3:=Dagger(kets3);
print("Matrix elements computed with C appears curious");
'bras3.C. kets3'="...";
bras3.C.kets3;
print("Matrix elements computed with Cn as expected");
'bras3.Cn.kets3'=bras3.Cn.kets3;

Physics:-`*`(Physics:-Ket(a1), Physics:-Ket(a2), Physics:-Ket(a3), Physics:-Ket(a4), Physics:-Ket(a5), Physics:-Ket(b1), Physics:-Ket(b2), Physics:-Ket(b3))

  

Physics:-`*`(Physics:-Bra(b3), Physics:-Bra(b2), Physics:-Bra(b1), Physics:-Bra(a5), Physics:-Bra(a4), Physics:-Bra(a3), Physics:-Bra(a2), Physics:-Bra(a1))

  

"Matrix elements computed with C appears curious"

  

Physics:-`.`(bras3, C, kets3) = "..."

  

Physics:-Bracket(Physics:-Bra(b1), Physics:-Ket(a1))*Physics:-Bracket(Physics:-Bra(b2), Physics:-Ket(a2))*Physics:-Bracket(Physics:-Bra(b3), Physics:-Ket(a3))*Physics:-Bracket(Physics:-Bra(a1), Physics:-Ket(b1))*Physics:-Bracket(Physics:-Bra(a2), Physics:-Ket(b2))*Physics:-Bracket(Physics:-Bra(a3), Physics:-Ket(b3))*Physics:-Bracket(Physics:-Bra(a4), Physics:-Ket(a4))^2*Physics:-Bracket(Physics:-Bra(a5), Physics:-Ket(a5))^2*Physics:-Bracket(Physics:-Bra(a1), Physics:-Ket(a1))*Physics:-Bracket(Physics:-Bra(a2), Physics:-Ket(a2))*Physics:-Bracket(Physics:-Bra(a3), Physics:-Ket(a3))*Physics:-Bracket(Physics:-Bra(b1), Physics:-Ket(b1))*Physics:-Bracket(Physics:-Bra(b2), Physics:-Ket(b2))*Physics:-Bracket(Physics:-Bra(b3), Physics:-Ket(b3))

  

"Matrix elements computed with Cn as expected"

  

Physics:-`.`(bras3, Cn, kets3) = Physics:-Bracket(Physics:-Bra(a1), Physics:-Ket(a1))^2*Physics:-Bracket(Physics:-Bra(a2), Physics:-Ket(a2))^2*Physics:-Bracket(Physics:-Bra(a3), Physics:-Ket(a3))^2*Physics:-Bracket(Physics:-Bra(a4), Physics:-Ket(a4))^2*Physics:-Bracket(Physics:-Bra(a5), Physics:-Ket(a5))^2*Physics:-Bracket(Physics:-Bra(b1), Physics:-Ket(b1))^2*Physics:-Bracket(Physics:-Bra(b2), Physics:-Ket(b2))^2*Physics:-Bracket(Physics:-Bra(b3), Physics:-Ket(b3))^2

 (4)
> 

print("Example");
En:=ExpandQop(10*(1-x+y+z)*i*(1/sqrt(2))*A*B);

"Example"

  

-(5*I)*2^(1/2)*(-1+x-y-z)*Physics:-`*`(Physics:-Ket(a1), Physics:-Ket(a2), Physics:-Ket(a3), Physics:-Ket(a4), Physics:-Ket(a5), Physics:-Ket(b1), Physics:-Ket(b2), Physics:-Ket(b3), Physics:-Bra(b3), Physics:-Bra(b2), Physics:-Bra(b1), Physics:-Bra(a5), Physics:-Bra(a4), Physics:-Bra(a3), Physics:-Bra(a2), Physics:-Bra(a1))

 (5)
> 

print("Another example");
'F'='A*B/sqrt(2)+B*A/sqrt(2)';
F:=A*B/sqrt(2)+B*A/sqrt(2):
'op(1,F)'=op(1,F);
'op(2,F)'=op(2,F);

'Fn'='ExpandQop(F)';
Fn:=ExpandQop(F):
'op(1,Fn)'=op(1,Fn);
'op(2,Fn)'=op(2,Fn);

"Another example"

  

F = Physics:-`*`(Physics:-`*`(A, B), Physics:-`^`(sqrt(2), -1))+Physics:-`*`(Physics:-`*`(B, A), Physics:-`^`(sqrt(2), -1))

  

op(1, F) = (1/2)*2^(1/2)*Physics:-`*`(Physics:-Ket(a1), Physics:-Ket(a2), Physics:-Ket(a3), Physics:-Ket(a4), Physics:-Ket(a5), Physics:-Bra(a5), Physics:-Bra(a4), Physics:-Bra(a3), Physics:-Bra(a2), Physics:-Bra(a1), Physics:-Ket(b1), Physics:-Ket(b2), Physics:-Ket(b3), Physics:-Bra(b3), Physics:-Bra(b2), Physics:-Bra(b1))

  

op(2, F) = (1/2)*2^(1/2)*Physics:-`*`(Physics:-Ket(b1), Physics:-Ket(b2), Physics:-Ket(b3), Physics:-Bra(b3), Physics:-Bra(b2), Physics:-Bra(b1), Physics:-Ket(a1), Physics:-Ket(a2), Physics:-Ket(a3), Physics:-Ket(a4), Physics:-Ket(a5), Physics:-Bra(a5), Physics:-Bra(a4), Physics:-Bra(a3), Physics:-Bra(a2), Physics:-Bra(a1))

  

Fn = ExpandQop(F)

  

op(1, Fn) = (1/2)*2^(1/2)*Physics:-`*`(Physics:-Ket(a1), Physics:-Ket(a2), Physics:-Ket(a3), Physics:-Ket(a4), Physics:-Ket(a5), Physics:-Ket(b1), Physics:-Ket(b2), Physics:-Ket(b3), Physics:-Bra(b3), Physics:-Bra(b2), Physics:-Bra(b1), Physics:-Bra(a5), Physics:-Bra(a4), Physics:-Bra(a3), Physics:-Bra(a2), Physics:-Bra(a1))

  

op(2, Fn) = (1/2)*2^(1/2)*Physics:-`*`(Physics:-Ket(b1), Physics:-Ket(b2), Physics:-Ket(b3), Physics:-Ket(a1), Physics:-Ket(a2), Physics:-Ket(a3), Physics:-Ket(a4), Physics:-Ket(a5), Physics:-Bra(a5), Physics:-Bra(a4), Physics:-Bra(a3), Physics:-Bra(a2), Physics:-Bra(a1), Physics:-Bra(b3), Physics:-Bra(b2), Physics:-Bra(b1))

 (6)
> 

print("Final example, multiple products");
G:=B*B*B;
'G'=ExpandQop(G);

"Final example, multiple products"

  

Physics:-`*`(Physics:-Ket(b1), Physics:-Ket(b2), Physics:-Ket(b3), Physics:-Bra(b3), Physics:-Bra(b2), Physics:-Bra(b1), Physics:-Ket(b1), Physics:-Ket(b2), Physics:-Ket(b3), Physics:-Bra(b3), Physics:-Bra(b2), Physics:-Bra(b1), Physics:-Ket(b1), Physics:-Ket(b2), Physics:-Ket(b3), Physics:-Bra(b3), Physics:-Bra(b2), Physics:-Bra(b1))

  

G = Physics:-`*`(Physics:-Ket(b1), Physics:-Ket(b2), Physics:-Ket(b3), Physics:-Ket(b1), Physics:-Ket(b2), Physics:-Ket(b3), Physics:-Ket(b1), Physics:-Ket(b2), Physics:-Ket(b3), Physics:-Bra(b3), Physics:-Bra(b2), Physics:-Bra(b1), Physics:-Bra(b3), Physics:-Bra(b2), Physics:-Bra(b1), Physics:-Bra(b3), Physics:-Bra(b2), Physics:-Bra(b1))

 (7)
> 

 


 

Download ExpandQopV3.mw

Automatic handling of collision of tensor indices...

by: ecterrab 13431 Maple
rule
March 07 2018
3  0
 

Automatic handling of collision of tensor indices in products

 

 

The design of products of tensorial expressions that have contracted indices got enhanced. The idea: repeated indices in certain subexpressions are actually dummies. So suppose T[a, b] and B[b] are tensors, then in T[trace] = T[a, `~a`], a is just a dummy, therefore T[a, `~a`]*B[a] = T[b, `~b`]*B[a] is a well defined object. The new design automatically maps input like T[a, `~a`]*B[a] into T[b, `~b`]*B[a]. This significantly simplifies the manipulation of indices when working with products of tensorial expressions: each tensorial expression being multiplied has its repeated indices automatically transformed into different ones when they would collide with the free or repeated indices of the other expressions being multiplied.

 

This functionality is available within the Physics update distributed at the Maplesoft R&D Physics webpage (but for what you see under Preview that makes use of a new kernel feature of the Maple version under development).

 

> 

restart

> 

with(Physics); Setup(spacetimeindices = lowercaselatin, quiet)

[spacetimeindices = lowercaselatin]

(1)
> 

Define(T[a, b], B[b])

`Defined objects with tensor properties`

 

{B[b], Physics:-Dgamma[a], Physics:-Psigma[a], T[a, b], Physics:-d_[a], Physics:-g_[a, b], Physics:-KroneckerDelta[a, b], Physics:-LeviCivita[a, b, c, d]}

(2)

This shows the automatic handling of collision of indices

> 

T[a, a]*B[a]

T[b, `~b`]*B[a]

(3)
> 

T[a, a]^2

T[a, `~a`]*T[b, `~b`]

(4)
> 

``

Preview only in the upcomming version under development

 

Consider now the case of three tensors

> 

Define(A[a], C[a])

`Defined objects with tensor properties`

 

{A[a], B[b], C[a], Physics:-Dgamma[a], Physics:-Psigma[a], T[a, b], Physics:-d_[a], Physics:-g_[a, b], Physics:-KroneckerDelta[a, b], Physics:-LeviCivita[a, b, c, d]}

(5)
> 

A[a]*B[a]*C[a]

A[a]*B[a]*C[a]

(6)

The product above has indeed the index a repeated more than once, therefore none of its occurrences got automatically transformed into contravariant in the output, and Check  detects the problem interrupting with an error  message

> 

Check(A[a]*B[a]*C[a])

Error, (in Physics:-Check) wrong use of the summation rule for repeated indices: `a repeated 3 times`, in A[a]*B[a]*C[a]

 

 

However, it is now also possible to indicate, using parenthesis, that the product of two of these tensors actually form a subexpression, so that the following two tensorial expressions are well defined, and the colliding dummy is automatically replaced making that explicit

> 

A[a]*B[a]*C[a]

A[b]*B[`~b`]*C[a]

(7)
> 

A[a]*B[a]*C[a]

A[a]*B[b]*C[`~b`]

(8)

 

 

This change in design makes concretely simpler the use of indices in that it eliminates the need for manually replacing dummies. For example, consider the tensorial expression for the angular momentum in terms of the coordinates and momentum vectors, in 3 dimensions

 

> 

Setup(coordinates = cartesian, dimension = 3, metric = euclidean, quiet)

[coordinatesystems = {X}, dimension = 3, metric = {(1, 1) = 1, (2, 2) = 1, (3, 3) = 1}]

(9)

Define L[j], p[k] respectively representing angular and linear momentum

> 

Define(L[j], p[k])

`Defined objects with tensor properties`

 

{Physics:-Dgamma[a], L[j], Physics:-Psigma[a], Physics:-d_[a], Physics:-g_[a, b], p[k], Physics:-KroneckerDelta[a, b], Physics:-LeviCivita[a, b, c], Physics:-SpaceTimeVector[a](X)}

(10)

Introduce the tensorial expression for L[a]

> 

L[a] = LeviCivita[a, b, c]*X[b]*p[c]

L[a] = Physics:-LeviCivita[a, b, c]*Physics:-SpaceTimeVector[b](X)*p[c]

(11)

The left-hand side has one free index, a, while the right-hand side has two dummy indices b and c

> 

Check(L[a] = Physics[LeviCivita][a, b, c]*Physics[SpaceTimeVector][b](X)*p[c], all)

`The repeated indices per term are: `[{`...`}, {`...`}, `...`]*`; the free indices are: `*{`...`}

 

([{}], {a}) = ([{b, c}], {a})

(12)

If we want to compute`#mrow(msup(mfenced(mover(mi("L"),mo("&rarr;")),open = "&Vert;",close = "&Vert;"),mn("2")),mo("&equals;"),msubsup(mi("L"),mi("a"),mn("2")))`we can now take the square of (11) directly, and the dummy indices on the right-hand side are automatically handled, there is now no need to manually substitute the repeated indices to avoid their collision

> 

(L[a] = Physics[LeviCivita][a, b, c]*Physics[SpaceTimeVector][b](X)*p[c])^2

L[a]^2 = Physics:-LeviCivita[a, b, c]*Physics:-SpaceTimeVector[b](X)*p[c]*Physics:-LeviCivita[a, d, e]*Physics:-SpaceTimeVector[d](X)*p[e]

(13)

The repeated indices on the right-hand side are now a, b, c, d, e

> 

Check(L[a]^2 = Physics[LeviCivita][a, b, c]*Physics[SpaceTimeVector][b](X)*p[c]*Physics[LeviCivita][a, d, e]*Physics[SpaceTimeVector][d](X)*p[e], all)

`The repeated indices per term are: `[{`...`}, {`...`}, `...`]*`; the free indices are: `*{`...`}

 

([{a}], {}) = ([{a, b, c, d, e}], {})

(14)
> 

NULL


 

Download AutomaticHandlingCollisionOfTensorIndices.mw

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

Procedure for expanding tensor product of quantum...

by: Louis Lamarche 105 Maple
physics tensor quantum-mechanics
March 07 2018
2  8

Here is a major upgrade of the procedure I submitted in february.

https://www.mapleprimes.com/posts/209030-Procedure-For-Computing-The-Tensor-Product

There is a line after the procedure to save it in the file "ExpandQop.m"
In future post I will use it in order to minimize the size of the examples.

Louis Lamarche
 

> 

######################################################################
# NOTICE                                                             #
# Author: Louis Lamarche                                             #
#         Institute of Research of Hydro-Quebec (IREQ)               #
#         Science des données et haute performance                   #
#         2018, March 7                                              #
#                                                                    #
# Function name: ExpandQop (x)                                       #
#       Purpose: Compute the tensor product of two quantum           #
#                operators in Dirac notations                        #
#      Argument: x: a quantum operator                               #
#  Improvements: Manage all +, -, *, /, ^, mod  operations           #
#                in the argument                                     #
#       Version: 2                                                   #
#                                                                    #
#  Copyrigth(c) Hydro-Quebec.                                        #
#        Note 1: Permission to use this softwate is granted if you   #
#                acknowledge its author and copyright                #
#        Note 2: Permission to copy this softwate is granted if you  #
#                leave this 21 lines notice intact. Thank you.       #
######################################################################
restart;
> 

with(Physics):
interface(imaginaryunit=i):
Setup(mathematicalnotation=true);

[mathematicalnotation = true]

 (1)
> 

Setup(quantumoperators={A,B,C,Cn});
Setup(noncommutativeprefix={a,b,q});

[quantumoperators = {A, B, C, Cn}]

  

[noncommutativeprefix = {a, b, q}]

 (2)
> 

opexp:= op(0,10^x):            # exponentiation id
opnp := op(0,10*x):            # normal product id
optp := op(0,Ket(q0)*Ket(q1)): # tensor product id
opdiv:= `Fraction`:            # fraction       id          
opsym:= op(0,x):               # symbol         id
opint:= op(0,10):              # integer        id
opflt:= op(0,10.0):            # float          id
opcpx:= op(0,i):               # complex        id
opbra:= op(0,Bra(q)):          # bra            id
opket:= op(0,Ket(q)):          # ket            id
opmod:= op(0, a mod b):        # mod            id
ExpandQop:=proc(x)
    local ret,j,lkb,cbk,rkb,no,lop;
    lkb:=0; cbk:=0; rkb:=0;
    lop:=op(0,x);
    no:=nops(x);
    if lop = opsym or lop = opint or lop = opflt or
       lop = opbra or lop = opket or lop = opcpx then
         ret:=x;
    else
    if lop = opexp then
        ret:=x;
    else       
    if lop = opnp then
        ret:=1;
        for j from 1 to no do
            ret:=ret*ExpandQop(op(j,x));
        end do;        
    else
    if lop = `+` then
        ret:=0;
        for j from 1 to no do
            ret:=ret+ExpandQop(op(j,x));
        end do;
    else
    if lop = `-` then
        ret:=0;
        for j from 1 to no do
            ret:=ret-ExpandQop(op(j,x));
        end do;
    else
    if lop = opdiv then
       ret:=1;
       for j from 1 to no do
           ret:=ret/ExpandQop(op(j,x));
       end do;
    else
    if lop = opmod then
       ret:=x;
    else
    if lop = optp then
        ret:=1;
       if (no > 3 ) then
           for j from 1 to no do
               if (j>1) then
                    if(lkb=0) then
                        if( type(op(j-1,x),Ket) and
                            type(op(j,x),Bra) ) then lkb:=j-1; fi;
                    else
                        if( type(op(j-1,x),Ket) and
                            type(op(j,x),Bra) ) then rkb:=j;   fi;
                    fi;
                    if( type(op(j-1,x),Bra) and type(op(j,x),Ket) )
                                                then cbk:=j;   fi;
               fi;
           end do;
           if ( (lkb < cbk) and (cbk<rkb) ) then
               for j from 1     to lkb   do ret := ret*op(j,x); end do;
               for j from cbk   to no    do ret := ret*op(j,x); end do;
               for j from lkb+1 to cbk-1 do ret := ret*op(j,x); end do;
           else
               ret:=x;
           fi;
       else
           ret:=x;
       fi;
    else ret:=x;
    fi; # optp
    fi; # opmod
    fi; # opdiv
    fi; # `-`
    fi; # `+`
    fi; # `opnp
    fi; # `opexp`
    fi; # opsym, opint, opflt, opbra, opket, opcpx

    return ret;
end proc:

# For saving
# save opexp,opnp,optp,opdiv,opint,opflt,opcpx,opbra,opket,opmod, ExpandQop,"ExpandQop.m"
> 

# Let A be an operator in a first Hilbert space of dimension n
#  using the associated orthonormal ket and bra vectors
#
#
kets1:=Ket(a1)*Ket(a2)*Ket(a3)*Ket(a4)*Ket(a5):
A:=kets1*Dagger(kets1);


# Let B be an operator in a second Hilbert (Ketspace of dimension m
#  using the associated orthonormal ket and bra vectors
#
#
kets2:=Ket(b1)*Ket(b2)*Ket(b3):
B:=kets2*Dagger(kets2);


# The tensor product of the two operators acts on a n+m third
# Hilbert space   unsing the appropriately ordered ket
# and bra  vectors of the two preceding spaces. The rule for
# building this operator in Dirac notation is as follows,
#
#


print("Maple do not compute the tensor product of operators,");
print("C=A*B gives:");
C:=A*B;

print("ExpandQop(C) gives the expected result:");
Cn:=ExpandQop(C);

Physics:-`*`(Physics:-Ket(a1), Physics:-Ket(a2), Physics:-Ket(a3), Physics:-Ket(a4), Physics:-Ket(a5), Physics:-Bra(a5), Physics:-Bra(a4), Physics:-Bra(a3), Physics:-Bra(a2), Physics:-Bra(a1))

  

Physics:-`*`(Physics:-Ket(b1), Physics:-Ket(b2), Physics:-Ket(b3), Physics:-Bra(b3), Physics:-Bra(b2), Physics:-Bra(b1))

  

"Maple do not compute the tensor product of operators,"

  

"C=A*B gives:"

  

Physics:-`*`(Physics:-Ket(a1), Physics:-Ket(a2), Physics:-Ket(a3), Physics:-Ket(a4), Physics:-Ket(a5), Physics:-Bra(a5), Physics:-Bra(a4), Physics:-Bra(a3), Physics:-Bra(a2), Physics:-Bra(a1), Physics:-Ket(b1), Physics:-Ket(b2), Physics:-Ket(b3), Physics:-Bra(b3), Physics:-Bra(b2), Physics:-Bra(b1))

  

"ExpandQop(C) gives the expected result:"

  

Physics:-`*`(Physics:-Ket(a1), Physics:-Ket(a2), Physics:-Ket(a3), Physics:-Ket(a4), Physics:-Ket(a5), Physics:-Ket(b1), Physics:-Ket(b2), Physics:-Ket(b3), Physics:-Bra(b3), Physics:-Bra(b2), Physics:-Bra(b1), Physics:-Bra(a5), Physics:-Bra(a4), Physics:-Bra(a3), Physics:-Bra(a2), Physics:-Bra(a1))

 (3)
> 

kets3:=kets1*kets2;
bras3:=Dagger(kets3);
print("Matrix elements computed with C appears curious");
'bras3.C. kets3'="...";
bras3.C.kets3;
print("Matrix elements computed with Cn as expected");
'bras3.Cn.kets3'=bras3.Cn.kets3;

Physics:-`*`(Physics:-Ket(a1), Physics:-Ket(a2), Physics:-Ket(a3), Physics:-Ket(a4), Physics:-Ket(a5), Physics:-Ket(b1), Physics:-Ket(b2), Physics:-Ket(b3))

  

Physics:-`*`(Physics:-Bra(b3), Physics:-Bra(b2), Physics:-Bra(b1), Physics:-Bra(a5), Physics:-Bra(a4), Physics:-Bra(a3), Physics:-Bra(a2), Physics:-Bra(a1))

  

"Matrix elements computed with C appears curious"

  

Physics:-`.`(bras3, C, kets3) = "..."

  

Physics:-Bracket(Physics:-Bra(b1), Physics:-Ket(a1))*Physics:-Bracket(Physics:-Bra(b2), Physics:-Ket(a2))*Physics:-Bracket(Physics:-Bra(b3), Physics:-Ket(a3))*Physics:-Bracket(Physics:-Bra(a1), Physics:-Ket(b1))*Physics:-Bracket(Physics:-Bra(a2), Physics:-Ket(b2))*Physics:-Bracket(Physics:-Bra(a3), Physics:-Ket(b3))*Physics:-Bracket(Physics:-Bra(a4), Physics:-Ket(a4))^2*Physics:-Bracket(Physics:-Bra(a5), Physics:-Ket(a5))^2*Physics:-Bracket(Physics:-Bra(a1), Physics:-Ket(a1))*Physics:-Bracket(Physics:-Bra(a2), Physics:-Ket(a2))*Physics:-Bracket(Physics:-Bra(a3), Physics:-Ket(a3))*Physics:-Bracket(Physics:-Bra(b1), Physics:-Ket(b1))*Physics:-Bracket(Physics:-Bra(b2), Physics:-Ket(b2))*Physics:-Bracket(Physics:-Bra(b3), Physics:-Ket(b3))

  

"Matrix elements computed with Cn as expected"

  

Physics:-`.`(bras3, Cn, kets3) = Physics:-Bracket(Physics:-Bra(a1), Physics:-Ket(a1))^2*Physics:-Bracket(Physics:-Bra(a2), Physics:-Ket(a2))^2*Physics:-Bracket(Physics:-Bra(a3), Physics:-Ket(a3))^2*Physics:-Bracket(Physics:-Bra(a4), Physics:-Ket(a4))^2*Physics:-Bracket(Physics:-Bra(a5), Physics:-Ket(a5))^2*Physics:-Bracket(Physics:-Bra(b1), Physics:-Ket(b1))^2*Physics:-Bracket(Physics:-Bra(b2), Physics:-Ket(b2))^2*Physics:-Bracket(Physics:-Bra(b3), Physics:-Ket(b3))^2

 (4)
> 

print("Example");
En:=ExpandQop(10*(1-x+y+z)*i*(1/sqrt(2))*A*B);

"Example"

  

-(5*I)*2^(1/2)*(-1+x-y-z)*Physics:-`*`(Physics:-Ket(a1), Physics:-Ket(a2), Physics:-Ket(a3), Physics:-Ket(a4), Physics:-Ket(a5), Physics:-Ket(b1), Physics:-Ket(b2), Physics:-Ket(b3), Physics:-Bra(b3), Physics:-Bra(b2), Physics:-Bra(b1), Physics:-Bra(a5), Physics:-Bra(a4), Physics:-Bra(a3), Physics:-Bra(a2), Physics:-Bra(a1))

 (5)
> 

print("Another example");
'F'='A*B/sqrt(2)+B*A/sqrt(2)';
F:=A*B/sqrt(2)+B*A/sqrt(2):
'op(1,F)'=op(1,F);
'op(2,F)'=op(2,F);

'Fn'='ExpandQop(F)';
Fn:=ExpandQop(F):
'op(1,Fn)'=op(1,Fn);
'op(2,Fn)'=op(2,Fn);

"Another example"

  

F = Physics:-`*`(Physics:-`*`(A, B), Physics:-`^`(sqrt(2), -1))+Physics:-`*`(Physics:-`*`(B, A), Physics:-`^`(sqrt(2), -1))

  

op(1, F) = (1/2)*2^(1/2)*Physics:-`*`(Physics:-Ket(a1), Physics:-Ket(a2), Physics:-Ket(a3), Physics:-Ket(a4), Physics:-Ket(a5), Physics:-Bra(a5), Physics:-Bra(a4), Physics:-Bra(a3), Physics:-Bra(a2), Physics:-Bra(a1), Physics:-Ket(b1), Physics:-Ket(b2), Physics:-Ket(b3), Physics:-Bra(b3), Physics:-Bra(b2), Physics:-Bra(b1))

  

op(2, F) = (1/2)*2^(1/2)*Physics:-`*`(Physics:-Ket(b1), Physics:-Ket(b2), Physics:-Ket(b3), Physics:-Bra(b3), Physics:-Bra(b2), Physics:-Bra(b1), Physics:-Ket(a1), Physics:-Ket(a2), Physics:-Ket(a3), Physics:-Ket(a4), Physics:-Ket(a5), Physics:-Bra(a5), Physics:-Bra(a4), Physics:-Bra(a3), Physics:-Bra(a2), Physics:-Bra(a1))

  

Fn = ExpandQop(F)

  

op(1, Fn) = (1/2)*2^(1/2)*Physics:-`*`(Physics:-Ket(a1), Physics:-Ket(a2), Physics:-Ket(a3), Physics:-Ket(a4), Physics:-Ket(a5), Physics:-Ket(b1), Physics:-Ket(b2), Physics:-Ket(b3), Physics:-Bra(b3), Physics:-Bra(b2), Physics:-Bra(b1), Physics:-Bra(a5), Physics:-Bra(a4), Physics:-Bra(a3), Physics:-Bra(a2), Physics:-Bra(a1))

  

op(2, Fn) = (1/2)*2^(1/2)*Physics:-`*`(Physics:-Ket(b1), Physics:-Ket(b2), Physics:-Ket(b3), Physics:-Ket(a1), Physics:-Ket(a2), Physics:-Ket(a3), Physics:-Ket(a4), Physics:-Ket(a5), Physics:-Bra(a5), Physics:-Bra(a4), Physics:-Bra(a3), Physics:-Bra(a2), Physics:-Bra(a1), Physics:-Bra(b3), Physics:-Bra(b2), Physics:-Bra(b1))

 (6)
> 

 


 

Download ExpandQopV2.mw

Minimize the number of tensor components according...

by: ecterrab 13431 Maple
tensor physics
March 05 2018
5  0

 

Minimize the number of tensor components according to its symmetries
(and relabel, redefine or count the number of independent tensor components)

 

 

The nice development described below is work in collaboration with Pascal Szriftgiser from Laboratoire PhLAM, Université Lille 1, France, used in the Mapleprimes post Magnetic traps in cold-atom physics

 

A new keyword in Define  and Setup : minimizetensorcomponents, allows for automatically minimizing the number of tensor components taking into account the tensor symmetries. For example, if a tensor with two indices in a 4D spacetime is defined as antisymmetric using Define with this new keyword, the number of different tensor components will be exactly 6, and the elements of the diagonal are automatically set equal to 0. After setting this keyword to true with Setup , all subsequent definitions of tensors automatically minimize the number of components while using this keyword with Define  makes this minimization only happen with the tensors being defined in the call to Define .

 

Related to this new functionality, 4 new Library routines were added: MinimizeTensorComponents, NumberOfIndependentTensorComponents, RelabelTensorComponents and RedefineTensorComponents

 

Example:

> 

restart; with(Physics)

 

Define an antisymmetric tensor with two indices

> 

Define(F[mu, nu], antisymmetric)

`Defined objects with tensor properties`

 

{Physics:-Dgamma[mu], F[mu, nu], Physics:-Psigma[mu], Physics:-d_[mu], Physics:-g_[mu, nu], Physics:-KroneckerDelta[mu, nu], Physics:-LeviCivita[alpha, beta, mu, nu]}

(1.1)

Although the system knows that F[mu, nu] is antisymmetric, you need to use Simplify to apply the (anti)symmetry

> 

F[mu, nu]+F[nu, mu]

F[mu, nu]+F[nu, mu]

(1.2)

 

> 

Simplify(F[mu, nu]+F[nu, mu])

0

(1.3)

so by default the components of F[mu, nu] do not automatically reflect the (anti)symmetry; likewise

> 

F[1, 2]+F[2, 1]

F[1, 2]+F[2, 1]

(1.4)
> 

Simplify(F[1, 2]+F[2, 1])

0

(1.5)

and computing the array form of F[mu, nu]we do not see the elements of the diagonal equal to zero nor the lower-left triangle equal to the upper-right triangle but for a different sign:

> 

TensorArray(F[mu, nu])

Matrix(%id = 18446744078270093062)

(1.6)

 

On the other hand, this new functionality, here called minimizetensorcomponents, makes the symmetries of the tensor be explicitly reflected in its components.

 

There are three ways to use it. First, one can minimize the number of tensor components of a tensor previously defined. For example

 

> 

Library:-MinimizeTensorComponents(F)

Matrix(%id = 18446744078270064630)

(1.7)

After this, both (1.2) and (1.3) are automatically equal to 0 without having to use Simplify

> 

F[mu, nu]+F[nu, mu]

0

(1.8)
> 

0

0

(1.9)

And the output of TensorArray  in (1.6) becomes equal to (1.7).

 

NOTE: in addition, after using minimizetensorcomponents in the definition of a tensor, say F, all the keywords implemented for Physics tensors are available for F:

 

> 

F[]

F[mu, nu] = Matrix(%id = 18446744078247910206)

(1.10)
> 

F[trace]

0

(1.11)
> 

F[nonzero]

F[mu, nu] = {(1, 2) = F[1, 2], (1, 3) = F[1, 3], (1, 4) = F[1, 4], (2, 1) = -F[1, 2], (2, 3) = F[2, 3], (2, 4) = F[2, 4], (3, 1) = -F[1, 3], (3, 2) = -F[2, 3], (3, 4) = F[3, 4], (4, 1) = -F[1, 4], (4, 2) = -F[2, 4], (4, 3) = -F[3, 4]}

 (1.12)
> 

"F[~1,mu,matrix]"

F[`~1`, mu] = Vector[row](%id = 18446744078247885990)

(1.13)

Alternatively, one can define a tensor, specifying that the symmetries should be taken into account to minimize the number of its components passing the keyword minimizetensorcomponents to Define .

 

Example:

 

Define a tensor with the symmetries of the Riemann  tensor, that is, a tensor of 4 indices that is symmetric with respect to interchanging the positions of the 1st and 2nd pair of indices and antisymmetric with respect to interchanging the position of its 1st and 2nd indices, or 3rd and 4th indices, and define it minimizing the number of tensor components

 

> 

Define(R[alpha, beta, mu, nu], symmetric = {[[1, 2], [3, 4]]}, antisymmetric = {[1, 2], [3, 4]}, minimizetensorcomponents)

`Defined objects with tensor properties`

 

{Physics:-Dgamma[mu], F[mu, nu], Physics:-Psigma[mu], R[mu, nu, alpha, beta], Physics:-d_[mu], Physics:-g_[mu, nu], Physics:-KroneckerDelta[mu, nu], Physics:-LeviCivita[alpha, beta, mu, nu]}

(1.14)

We now have

> 

R[1, 2, 3, 4]+R[2, 1, 3, 4]

0

(1.15)
> 

R[alpha, beta, mu, nu]-R[mu, nu, alpha, beta]

0

(1.16)
• 

One can always retrieve the symmetry properties in the abstract notation used by the Define command using the new Library:-GetTensorSymmetryProperties, its output is ordered, first the symmetric then the antisymmetric properties

 

> 

Library:-GetTensorSymmetryProperties(R)

{[[1, 2], [3, 4]]}, {[1, 2], [3, 4]}

(1.17)
• 

After making the symmetries explicit (and also before that), it is frequently useful to know the number of independent components of a given tensor. For this purpose you can use the new Library:-NumberOfIndependentTensorComponents

 

> 

Library:-NumberOfIndependentTensorComponents(R)

21

(1.18)

and besides taking into account the symmetries, in the case of the Riemann  tensor, after taking into account the first Bianchi identity this number of components is further reduced to 20.

 

A third way of using the new minimizetensorcomponents functionality is using Setup , so that, automatically, every subsequent definition of tensors with symmetries is performed minimizing the number of its components using the indicated symmetries

 

Example:

> 

Setup(minimizetensorcomponents = true)

[minimizetensorcomponents = true]

(1.19)

So from hereafter you can define tensors taking into account their symmetries explicitly and without having to include the keyword minimizetensorcomponents at each definition

 

> 

Define(C[alpha, beta], antisymmetric)

`Defined objects with tensor properties`

 

{C[mu, nu], Physics:-Dgamma[mu], F[mu, nu], Physics:-Psigma[mu], R[mu, nu, alpha, beta], Physics:-d_[mu], Physics:-g_[mu, nu], Physics:-KroneckerDelta[mu, nu], Physics:-LeviCivita[alpha, beta, mu, nu]}

(1.20)

 

> 

C[]

C[mu, nu] = Matrix(%id = 18446744078408747598)

(1.21)
• 

Two new related functionalities are provided via Library:-RelabelTensorComponents and Library:-RedefineTensorComponent, the first one to have the number of tensor components directly reflected in the names of the components, the second one to redefine only one of these components

> 

Library:-RelabelTensorComponents(C)

Matrix(%id = 18446744078408729774)

(1.22)

 

Suppose now we want to make one of these components equal to 1, say C__2

> 

Library:-RedefineTensorComponent(C[1, 2] = 1)

C[mu, nu] = Matrix(%id = 18446744078270104390)

(1.23)

This nice development is work in collaboration with Pascal Szriftgiser from Laboratoire PhLAM, UMR CNRS 8523, Université Lille 1, F-59655, France.

> 

``

 

Download MinimizeTensorComponents.mw

 

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

An example of the kinds of worksheets that i make...

by: Adam Ledger 360 Maple
number_theory
March 03 2018
1  0

 

Hello, 

I study mainly subjects that fall under umbrella of number theory, but i have specified a little further in the worksheet. This is really a request for assistance, because in as much as i have met so many brilliant people online via social media etc,  I would always love to meet more, and especially ones who are more experienced in this field. 

 

Basically i am too cheap and old to think about going to a good university, so I am trying to get free advice from the people who have probably completed doctorates in the relevant field. Got to be honest I say.

 

Anyway my contact email is at the top of the attached worksheet.

 

First thing that stood out to me about the distributions produced in this worksheet is how sparse the number of points is for N=17 relative to all the other values of N.

EXAMPLE_FOR_MAPLE3.mw

 

Edit: Another example worksheet added.

MAPLE_EXAMPLE_13.mw

Procedure for computing the tensor product of...

by: Louis Lamarche 105 Maple
procedure product physics
February 22 2018
2  0

I wanted to use MAPLE to preform symbolic quantum computations. The role
of quantum operators and their tensor product is very important in simplying
the understating of such new calculus at least for the beginners. For instance,
(using "o" for the tensor product and "." for the scalar product, H being the Hadamard
operator on a qubit, I the identity operator, and CNOT the 2 qubit controled not
operator)
1) generating the Bells states |Bxy> two stages of operators are needed
     (CNOT) .  (H o  I)  . |x> o |y>

2) performing quantum teleportation of |psi>
     (H o I o I) . (CNOT o I ) . |psi>o |B00>
    followed by a measurements on the first two qubits for driving the application of
    quantum gates to the third qubit.

All these tensor products of operators can be easily written in MAPLE.

Here is a first version of the ExpandQop procedure that will be usefull the purpose of
expanding correctly the tensor product of two quantum operator expressed in Dirac notation.

I hope this is usefull.

LL 

 

> 

######################################################################
# Author: Louis Lamarche                                             #
#         Institute of Research of Hydro-Quebec (IREQ)               #
#         Science des données et haute performance                   #
#         2018/02/20                                                 #
#                                                                    #
#         Function name: ExpandQop (x)                               #
#               Purpose: Compute the tensor product of two quantum   #
#                        operators in Dirac notations                #
#              Argument: x - a simple quantum operator               #
#   Future improvements: Manage all +, -, *, /, ^, mod  operations   #
#                        in the argument                             #
#               Version: 1.0                                         #
######################################################################
restart;
> 

with(Physics):
interface(imaginaryunit=i):
Setup(mathematicalnotation=true);

[mathematicalnotation = true]

 (1)
> 

Setup(quantumoperators={A,B,C,Cn});
Setup(noncommutativeprefix={a,b});

[quantumoperators = {A, B, C, Cn}]

  

[noncommutativeprefix = {a, b}]

 (2)
> 

ExpandQop:=proc(x)
    local ret,j,lkb,cbk,rkb,no;
    ret:=1; lkb:=0; cbk:=0; rkb:=0; no:=nops(x);
    if (no > 3 ) then
        for j from 1 to no do
            if (j>1) then
                 if(lkb=0) then
                     if( type(op(j-1,x),Ket) and
                         type(op(j,x),Bra) ) then lkb:=j-1; fi;
                 else
                     if( type(op(j-1,x),Ket) and
                         type(op(j,x),Bra) ) then rkb:=j;   fi;
                 fi;
                 if( type(op(j-1,x),Bra) and type(op(j,x),Ket) )
                                             then cbk:=j;   fi;
            fi;
        end do;
        if ( (lkb < cbk) and (cbk<rkb) ) then
            for j from 1     to lkb   do ret := ret*op(j,x); end do;
            for j from cbk   to no    do ret := ret*op(j,x); end do;
            for j from lkb+1 to cbk-1 do ret := ret*op(j,x); end do;
        else
            ret:=x;
        fi;
    else
        ret:=x;
    fi;
    return ret;
end proc:
> 

# Let A be an operator in a first Hilbert space of dimension n
#  using the associated orthonormal ket and bra vectors
#
#
kets1:=Ket(a1)*Ket(a2)*Ket(a3)*Ket(a4)*Ket(a5):
A:=kets1*Dagger(kets1);


# Let B be an operator in a second Hilbert (Ketspace of dimension m
#  using the associated orthonormal ket and bra vectors
#
#
kets2:=Ket(b1)*Ket(b2)*Ket(b3):
B:=kets2*Dagger(kets2);


# The tensor product of the two operators acts on a n+m third
# Hilbert space   unsing the appropriately ordered ket
# and bra  vectors of the two preceding spaces. The rule for
# building this operator in Dirac notation is as follows,
#
#


print("Maple do not compute the tensor product of operators,");
print("C=A*B gives:");
C:=A*B;

print("ExpandQop(C) gives the expected result:");
Cn:=ExpandQop(C);

Physics:-`*`(Physics:-Ket(a1), Physics:-Ket(a2), Physics:-Ket(a3), Physics:-Ket(a4), Physics:-Ket(a5), Physics:-Bra(a5), Physics:-Bra(a4), Physics:-Bra(a3), Physics:-Bra(a2), Physics:-Bra(a1))

  

Physics:-`*`(Physics:-Ket(b1), Physics:-Ket(b2), Physics:-Ket(b3), Physics:-Bra(b3), Physics:-Bra(b2), Physics:-Bra(b1))

  

"Maple do not compute the tensor product of operators,"

  

"C=A*B gives:"

  

Physics:-`*`(Physics:-Ket(a1), Physics:-Ket(a2), Physics:-Ket(a3), Physics:-Ket(a4), Physics:-Ket(a5), Physics:-Bra(a5), Physics:-Bra(a4), Physics:-Bra(a3), Physics:-Bra(a2), Physics:-Bra(a1), Physics:-Ket(b1), Physics:-Ket(b2), Physics:-Ket(b3), Physics:-Bra(b3), Physics:-Bra(b2), Physics:-Bra(b1))

  

"ExpandQop(C) gives the expected result:"

  

Physics:-`*`(Physics:-Ket(a1), Physics:-Ket(a2), Physics:-Ket(a3), Physics:-Ket(a4), Physics:-Ket(a5), Physics:-Ket(b1), Physics:-Ket(b2), Physics:-Ket(b3), Physics:-Bra(b3), Physics:-Bra(b2), Physics:-Bra(b1), Physics:-Bra(a5), Physics:-Bra(a4), Physics:-Bra(a3), Physics:-Bra(a2), Physics:-Bra(a1))

 (3)
> 

kets3:=kets1*kets2;
bras3:=Dagger(kets3);
print("Matrix element computed with C appears curious");
'bras3.C. kets3'="...";
bras3.C.kets3;
print("Matrix element computed with Cn as expected");
'bras3.Cn.kets3'=bras3.Cn.kets3;

Physics:-`*`(Physics:-Ket(a1), Physics:-Ket(a2), Physics:-Ket(a3), Physics:-Ket(a4), Physics:-Ket(a5), Physics:-Ket(b1), Physics:-Ket(b2), Physics:-Ket(b3))

  

Physics:-`*`(Physics:-Bra(b3), Physics:-Bra(b2), Physics:-Bra(b1), Physics:-Bra(a5), Physics:-Bra(a4), Physics:-Bra(a3), Physics:-Bra(a2), Physics:-Bra(a1))

  

"Matrix element computed with C"

  

Physics:-`.`(bras3, C, kets3) = "..."

  

Physics:-Bracket(Physics:-Bra(b1), Physics:-Ket(a1))*Physics:-Bracket(Physics:-Bra(b2), Physics:-Ket(a2))*Physics:-Bracket(Physics:-Bra(b3), Physics:-Ket(a3))*Physics:-Bracket(Physics:-Bra(a1), Physics:-Ket(b1))*Physics:-Bracket(Physics:-Bra(a2), Physics:-Ket(b2))*Physics:-Bracket(Physics:-Bra(a3), Physics:-Ket(b3))*Physics:-Bracket(Physics:-Bra(a4), Physics:-Ket(a4))^2*Physics:-Bracket(Physics:-Bra(a5), Physics:-Ket(a5))^2*Physics:-Bracket(Physics:-Bra(a1), Physics:-Ket(a1))*Physics:-Bracket(Physics:-Bra(a2), Physics:-Ket(a2))*Physics:-Bracket(Physics:-Bra(a3), Physics:-Ket(a3))*Physics:-Bracket(Physics:-Bra(b1), Physics:-Ket(b1))*Physics:-Bracket(Physics:-Bra(b2), Physics:-Ket(b2))*Physics:-Bracket(Physics:-Bra(b3), Physics:-Ket(b3))

  

"Matrix element computed with Cn as expected"

  

Physics:-`.`(bras3, Cn, kets3) = Physics:-Bracket(Physics:-Bra(a1), Physics:-Ket(a1))^2*Physics:-Bracket(Physics:-Bra(a2), Physics:-Ket(a2))^2*Physics:-Bracket(Physics:-Bra(a3), Physics:-Ket(a3))^2*Physics:-Bracket(Physics:-Bra(a4), Physics:-Ket(a4))^2*Physics:-Bracket(Physics:-Bra(a5), Physics:-Ket(a5))^2*Physics:-Bracket(Physics:-Bra(b1), Physics:-Ket(b1))^2*Physics:-Bracket(Physics:-Bra(b2), Physics:-Ket(b2))^2*Physics:-Bracket(Physics:-Bra(b3), Physics:-Ket(b3))^2

 (4)
> 

 


 

Download ExpandQop.mw

 

 

Different types of 'functions' -

by: das1404 130 Maple 7
eval evalf operator
February 22 2018
2  4

Lesson_on_functions.mws

As the title says, a lesson on functions:  eg the -> operator, f(2), eval, evalf etc 

Balance of a body with learning sequence with...

by: Lenin Araujo Castillo 1790 Maple 18
engineering programming
February 14 2018
5  0

Using the learning sequence as an alternative to learn problems related to "balance of a body" is shown in this video; thanks to the kindness that Maple offers us in its fundamental programming syntax.

Balance_of_a_body_with_learning_sequence.mw

Lenin Araujo

Ambassador Of Maple

 

 

Hyperasymptotics

by: _Maxim_ 724
hypergeom meijerg
February 14 2018
7  3

One case where an "expansion beyond all orders" may be needed is investigating the asymptotic behavior of the difference of two functions with coinciding dominant series.

We are interested in the asymptotic behavior of F(z) for large positive z:

> 

h1 := proc (z) options operator, arrow; hypergeom([1/2], [5/4, 3/2, 7/4], z) end proc; h2 := proc (z) options operator, arrow; (3/4)*sqrt(2*Pi)*hypergeom([1/4], [3/4, 5/4, 3/2], z)/(GAMMA(3/4)*z^(1/4)) end proc; F := proc (z) options operator, arrow; h1(z)-h2(z)+(3/8)*sqrt(Pi)/sqrt(z) end proc

series does not succeed:

> 

series(F(z), z = infinity, 20)

O((1/z)^(23/3))*exp(3/(1/z)^(1/3))

 (1)

The reason is that the dominant terms containing the factor exp(3*z^(1/3)) in the two hypergeometric functions cancel exactly, and we have to look for the subdominant terms.

The order of the leading terms can be found from DETools:-formal_sol:

> 

deq1 := FunctionAdvisor(DE, h1(z), f(z))[2, 1]

diff(diff(diff(diff(f(z), z), z), z), z) = -(15/2)*(diff(diff(diff(f(z), z), z), z))/z-(195/16)*(diff(diff(f(z), z), z))/z^2+(1/32)*(32*z-105)*(diff(f(z), z))/z^3+(1/2)*f(z)/z^3

 (2)
> 

DETools:-formal_sol(deq1, f(z), z = infinity, order = 0)

[(1/z)^(1/2), -exp(-3/(-1/z)^(1/3))/z, -exp(3*(-1)^(1/3)/(-1/z)^(1/3))/z, -exp(-3*(-1)^(2/3)/(-1/z)^(1/3))/z]

 (3)

As expected, one of the solutions (the third one for positive z) contains the exp(3*z^(1/3)) factor, the leading term being of the order exp(3*z^(1/3))/z.

Another, subdominant, solution is algebraic and, in fact, is a series containing only one term, as 1/z^(1/2) is an exact solution. It will turn out that the algebraic part in F(z) also cancels out.

Thus we have to look for the subsubdominant terms, which contain decaying exponentials. We will accomplish this by applying the steepest descent method to the integral representations of h1(z) and h2(z).

> 

ms := convert([h1(z), h2(z)], MeijerG)

[(3/32)*Pi*2^(1/2)*MeijerG([[1/2], []], [[0], [-1/4, -1/2, -3/4]], -z), (3/32)*2^(1/2)*Pi*MeijerG([[3/4], []], [[0], [1/4, -1/4, -1/2]], -z)/z^(1/4)]

 (4)
> 

m2g := proc (m, y) local a, b, c, d; a, b := op(op(1, m)); c, d := op(op(2, m)); -((1/2)*I)*mul(`~`[GAMMA](`~`[`-`](1+y, a)))*mul(`~`[GAMMA](`~`[`-`](c, y)))*op(3, m)^y/(Pi*mul(`~`[GAMMA](`~`[`-`](b, y)))*mul(`~`[GAMMA](`~`[`-`](1+y, d)))) end proc
> 

gs := applyrule(conditional(e::anything, _op(0, e) = MeijerG) = 'm2g(e, y)', ms)

[-((3/64)*I)*2^(1/2)*GAMMA(1/2+y)*GAMMA(-y)*(-z)^y/(GAMMA(5/4+y)*GAMMA(3/2+y)*GAMMA(7/4+y)), -((3/64)*I)*2^(1/2)*GAMMA(1/4+y)*GAMMA(-y)*(-z)^y/(z^(1/4)*GAMMA(3/4+y)*GAMMA(5/4+y)*GAMMA(3/2+y))]

 (5)
> 

gs[2] := combine(eval(gs[2], [1/z^(1/4) = exp(I*Pi*(1/4))/(-z)^(1/4), y = y+1/4]), power)

-((3/64)*I)*GAMMA(1/2+y)*((1/2)*2^(1/2)+((1/2)*I)*2^(1/2))*(-z)^y*GAMMA(-1/4-y)*2^(1/2)/(GAMMA(1+y)*GAMMA(3/2+y)*GAMMA(7/4+y))

 (6)

h1(z) and h2(z)are the integrals of gs[1] and of gs[2] over the same path, which is a loop encircling the poles ofGAMMA(-y) and of GAMMA(-1/4-y). Now a standard technique is to extend the integration contour far to the left, while still keeping both endpoints at "+infinity". Then the arguments of the gamma functions can be made large everywhere on the integration path, and the gamma functions can be replaced by their asymptotic approximations.

When moving the contour, we have to take into account the pole of the integrand at y = -1/2. The other poles of GAMMA(1/2+y) will be cancelled by the zeros of 1/GAMMA(3/2+y), which is why the algebraic part of the expansion will contain the single term of the order 1/z^(1/2).

This is the negative of the third term in F(z):

> 

`assuming`([simplify((2*Pi*I)*residue(gs[1]-gs[2], y = -1/2))], [z > 0])

-(3/8)*Pi^(1/2)/z^(1/2)

 (7)

Expanding the gamma functions produces terms containing exp(-I*Pi*y) and exp(I*Pi*y)

> 

`assuming`([simplify(convert(MultiSeries:-series((gs[1]-gs[2])/z^y, y = -infinity, 1), polynom))], [z > 0]); collect(convert(%, exp), exp)

(-3/64-(3/64)*I)*(-1/y)^(1/2)*exp(-3*(ln(-y)-1)*y)*exp(-I*Pi*y)/(y^3*Pi^(1/2))+(3/64-(3/64)*I)*(-1/y)^(1/2)*exp(-3*(ln(-y)-1)*y)*exp(I*Pi*y)/(y^3*Pi^(1/2))

 (8)

As we shall see, those terms have saddle points y0(z) = exp(`&+-`((1/3)*(2*Pi*I)))*z^(1/3) located in the left half-plane and contribute exponentially small factors exp(3*y0(z)). The terms for which the saddle point would be located at y = z^(1/3) have cancelled out, thus cancelling the exponentially large contributions. Another possible way to achieve the same result was to write h1(z)-h2(z) as a single Meijer G-function -(3/32)*MeijerG([[1/2], []], [[-1/4, 0], [-3/4, -1/2]], z).

We write the first term above in the form g(y)*exp(f(y)):

> 

f := proc (z, y) options operator, arrow; -3*y*(ln(-y)-1)-I*Pi*y+y*ln(z) end proc
> 

g := proc (y) options operator, arrow; (-3/64-(3/64)*I)*sqrt(-1/y)/(sqrt(Pi)*y^3) end proc
> 

diff(f(z, y), y)

-3*ln(-y)-I*Pi+ln(z)

 (9)

For this to become zero, we need argument(-y) = -(1/3)*Pi, and thus y = exp((1/3)*(2*I)*Pi)*z^(1/3). We can visualize the paths where the imaginary part of f(z, y) stays constant. The path of the steepest descent is the one that goes through the saddle point in the direction exp(I*Pi*(1/3)); the blue color indicates smaller values of the real part of f(z, y):

> 

y0 := proc (z) options operator, arrow; exp(((2/3)*I)*Pi)*z^(1/3) end proc
> 

(proc () local z; z := 2; plots:-display(plots:-contourplot(Re(f(z, u+I*v)), u = -5 .. 5, v = -5 .. 5, contours = ([seq])(Re(f(z, y0(z)))+i, i = -30 .. 6, 6), filledregions, coloring = [blue, red], grid = [100, 100]), plots:-implicitplot(Im(f(z, u+I*v)-f(z, y0(z))), u = -5 .. 5, v = -5 .. 5, gridrefine = 5, color = green), plot([cos((1/3)*Pi)*xi+Re(y0(z)), sin((1/3)*Pi)*xi+Im(y0(z)), xi = -3 .. 3], linestyle = dot, color = white), axes = boxed) end proc)()

 

The real part of f(z, y) has a maximum along this path at y0(z).

> 

`assuming`([(`@`(`@`(simplify, evalc), series))(f(z, y0(z)+exp(I*Pi*(1/3))*xi), xi = 0, 3)], [z > 0]); quad := convert(%, polynom)

series((3/2)*(-1+I*3^(1/2))*z^(1/3)-((3/2)/z^(1/3))*xi^2+O(xi^3),xi,3)

 (10)

Now we can compute the lead asymptotic term contributed by the saddle point y0(z):

> 

lt1 := `assuming`([(`@`(simplify, evalc))(g(y0(z))*exp(I*Pi*(1/3))*(int(exp(quad), xi = -infinity .. infinity)))], [z > 0])

-(1/64)*exp(-(3/2)*z^(1/3))*3^(1/2)*((-1+I)*cos((3/2)*z^(1/3)*3^(1/2))+(-1-I)*sin((3/2)*z^(1/3)*3^(1/2)))*2^(1/2)/z

 (11)

We repeat the same procedure for the second term of the integrand.

> 

f := proc (z, y) options operator, arrow; -3*y*(ln(-y)-1)+I*Pi*y+y*ln(z) end proc
> 

g := proc (y) options operator, arrow; (3/64-(3/64)*I)*sqrt(-1/y)/(sqrt(Pi)*y^3) end proc
> 

diff(f(z, y), y)

-3*ln(-y)+I*Pi+ln(z)

 (12)
> 

y0 := proc (z) options operator, arrow; exp(-((2/3)*I)*Pi)*z^(1/3) end proc

The direction should be chosen as exp((1/3)*(2*I)*Pi) to be consistent with the direction of the integration contour, which goes from the lower to the upper half-plane.

> 

lterm := proc (gy, fy, eq, dir) options operator, arrow; (eval(gy*exp(fy), eq))*dir*sqrt(-2*Pi/((eval(diff(fy, `$`(y, 2)), eq))*dir^2)) end proc
> 

lt2 := `assuming`([(`@`(simplify, evalc))(lterm(g(y), f(z, y), y = y0(z), exp((1/3)*(2*I)*Pi)))], [z > 0])

(1/64)*exp(-(3/2)*z^(1/3))*((1+I)*cos((3/2)*z^(1/3)*3^(1/2))+(1-I)*sin((3/2)*z^(1/3)*3^(1/2)))*3^(1/2)*2^(1/2)/z

 (13)

Combining the two results yields the leading term of F(z). The next terms can be obtained by expanding gs[1] and gs[2] to higher orders.

> 

Fasympt := unapply(simplify(lt1+lt2), z)

proc (z) options operator, arrow; (1/32)*exp(-(3/2)*z^(1/3))*3^(1/2)*2^(1/2)*(cos((3/2)*z^(1/3)*3^(1/2))+sin((3/2)*z^(1/3)*3^(1/2)))/z end proc

 (14)
> 

(proc () Digits := 50; plot(`~`[`*`](exp((3/2)*z^(1/3)), [F(z), Fasympt(z)]), z = 1000 .. 10000, linestyle = [solid, dot], thickness = [1, 5], axes = frame) end proc)()

 

Download steep.mw

First 21 22 23 24 25 26 27 Last Page 23 of 71