Advanced Search

Applications, Examples and Libraries

Share your work here
Create a new Post in this category

add, floats, and Kahan sum

by: vv 11619 Maple
evalf add
April 27 2018
3  2

add, floats, and Kahan sum

 

I found an intresting fact about the Maple command add for floating point values.
It seems that add in this case uses a summation algorithm in order to reduce the numerical error.
It is probably the Kahan summation algorithm (see wiki), but I wonder why this fact is not documented.

Here is a simple Maple procedure describing and implementing the algorithm.

 

 

> 

restart;
> 

Digits:=15;

15

 (1)
> 

KahanSum := proc(f::procedure, ab::range)  
local S,c,y,t, i;      # https://en.wikipedia.org/wiki/Kahan_summation_algorithm
S := 0.0;              # S = result (final sum: add(f(n), n=a..b))
c := 0.0;              # c = compensation for lost low-order bits.
for i from lhs(ab) to rhs(ab) do
    y := f(i) - c;     
    t := S + y;              
    c := (t - S) - y;        
    S := t;                  
od;                         
return S
end proc:

 

Now, a numerical example.

 

 

> 

f:= n ->  evalf(1/(n+1/n^3+1) - 1/(n+1+1/(n+1)^3+1));

proc (n) options operator, arrow; evalf(1/(n+1/n^3+1)-1/(n+2+1/(n+1)^3)) end proc

 (2)
> 

n := 50000;
K := KahanSum(f, 1..n);

50000

  

.333313334133301

 (3)
> 

A := add(f(k),k=1..n);

.333313334133302

 (4)
> 

s:=0.0:  for i to n do s:=s+f(i) od:
's' = s;

s = .333313334133413

 (5)
> 

exact:=( 1/3 - 1/(n+1+1/(n+1)^3+1) );

6250249999999900000/18751875067501050009

 (6)
> 

evalf( [errK = K-exact, errA = A-exact, err_for=s-exact] );

[errK = 0., errA = 0.1e-14, err_for = 0.112e-12]

 (7)
> 

evalf[20]( [errK = K-exact, errA = A-exact, err_for=s-exact] );

[errK = -0.33461e-15, errA = 0.66539e-15, err_for = 0.11166539e-12]

 (8)
> 

 


Download KahanSum.mw

STEM Education: Application generator in Maple

by: Lenin Araujo Castillo 1790 Maple 2018
April 18 2018
2  0

Due to the mechanistic process of our students and little creativity in analysis in schools and universities to be professionally trained is that STEM education appears (science, technology, engineering and mathematics) is a new model that is being considered in other countries and with very slow step in our city. In this work the methods with STEM will be visualized but using computational tools provided by Maplesoft which is a company that leads online education for adolescents and adults in the current market. In Spanish.

ECI_UNT_2018.pdf

ECI_UNT_2018.mw

Lenin Araujo Castillo

Ambassador of Maple

Bug corrected in inttrans[hilbert]

by: vv 11619 Maple
inttrans hilbert
April 10 2018
2  2

There is a bug in inttrans:-hilbert:

> 

restart;
> 

inttrans:-hilbert(sin(a)*sin(t+b), t, s);
# should be:
sin(a)*cos(s+b);   expand(%);

sin(a)*cos(s)

  

sin(a)*cos(s+b)

  

sin(a)*cos(s)*cos(b)-sin(a)*sin(s)*sin(b)

 (1)
> 

########## correction ##############
> 

`inttrans/expandc` := proc(expr, t)
local xpr, j, econst, op1, op2;
      xpr := expr;      
      for j in indets(xpr,specfunc(`+`,exp)) do
          econst := select(type,op(j),('freeof')(t));
          if 0 < nops(econst) and econst <> 0 then
              xpr := subs(j = ('exp')(econst)*combine(j/('exp')(econst),exp),xpr)
          end if
      end do;
      for j in indets(xpr,{('cos')(linear(t)), ('sin')(linear(t))}) do
          if type(op(j),`+`) then
              op1:=select(has, op(j),t); ##
              op2:=op(j)-op1;            ##
              #op1 := op(1,op(j));
              #op2 := op(2,op(j));
              if op(0,j) = sin then
                  xpr := subs(j = cos(op2)*sin(op1)+sin(op2)*cos(op1),xpr)
              else
                  xpr := subs(j = cos(op1)*cos(op2)-sin(op1)*sin(op2),xpr)
              end if
          end if
      end do;
      return xpr
end proc:
> 

#######################################
> 

inttrans:-hilbert(sin(a)*sin(t+b), t, s); expand(%);

-(1/2)*cos(a-b)*sin(s)+(1/2)*sin(a-b)*cos(s)+(1/2)*cos(a+b)*sin(s)+(1/2)*sin(a+b)*cos(s)

  

sin(a)*cos(s)*cos(b)-sin(a)*sin(s)*sin(b)

 (2)
> 

 


Download hilbert.mw

 

Recognizing Handwritten Digits with Machine Learni...

by: Samir Khan 1744 Maple 2018
April 09 2018
8  5

To demonstrate Maple 2018’s new Python connectivity, we wanted to integrate a large Python library. The result is the DeepLearning package - this offers an interface to a subset of the Tensorflow framework for machine learning.

I thought I’d share an application that demonstrates how the DeepLearning package can be used to recognize the numbers in images of handwritten digits.

The application employs a very small subset of the MNIST database of handwritten digits. Here’s a sample image for the digit 0.

This image can be represented as a matrix of pixel intensities.        

The application generates weights for each digit by training a two-layer neural network using multinomial logistic regression. When visualized, the weights for each digit might look like this.

Let’s say that we’re comparing an image of a handwritten digit to the weights for the digit 0. If a pixel with a high intensity lands in

  • an intensely red area, the evidence is high that the number in the image is 0
  • an intensely blue area, the evidence is low that the number in the image is 0

While this explanation is technically simplistic, the application offers more detail.

Get the application here

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 12095 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 12095 Maple
products 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 {`&Hscr;`[A], `&Hscr;`[B]} where the indexation indicates all the operators acting on that Hilbert space. The disjointedspaces keyword has for synonyms disjointedhilbertspaces and hilbertspaces. The display `&Hscr;`[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]),

 

`&otimes;`(Ket(A, alpha), Ket(B, beta)) = `&otimes;`(Ket(B, beta), Ket(A, alpha)) 

 

`&otimes;`(Bra(A, alpha), Ket(B, beta)) = `&otimes;`(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 5560 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 12095 Maple
products 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 12095 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)