Download Matrices Computation Homework and more Exercises Physics in PDF only on Docsity!
HOMEWORK 4, PHYS 562
Kam Modjtahedzadeh
January 22, 2021
Abstract
Linear algebra methods have been used to solve different different problems in physics and mathematics.
For the first problem, Gfortran and Python 3 are used separately to find Eigenvalues, Eigenvectors, and
exponentials of complex spin matrices. The physics was for this problem has largely been ignored 1 and
emphasis has been put on the computation. For the second problem, roots of a monic polynomial have
been found by finding the Eigenvalues of its companion matrix. The second problem was solved only
with GFortran (although I could have very easily done it with Python as well). The Lapack library was
used to solve each problem and it was imported in both languages.
1 Quantum Matrices
For this problem, the quantum spin number is j = 8. 5 , which means m varies starting from 8. 5 and going
down to − 8. 5 in increments of 1. The z component matrix of the spin J is given by the following relationship:
〈m
′ | Jz |m〉 = mδm′,m. (1)
The x and y matrices are constructed using the ± matrices, all given below;
〈m
′ | J± |m〉 =
j(j + 1) − m(m ± 1) δm,m± 1 (2)
Jx =
(J+ + J−) (3)
Jy =
2 i
(J+ − J−). (4)
The J 2 matrix is straight forward,
J
2 = J
2 x
+ J
2 y
+ J
2 z
Since the dimensions of the square matrices are m × m, which is 18 × 18 , they are too large to write by hand.
Hence, the following generating algorithms have been written in GFortran to construct (1)-(5).
1 Jz = n u l! i n i t i a l i z e t h e Jz m a t r i x
2 kk = j ∗ one! here , kk i s t h e e l e m e n t s o f Jz ( i. e. +m t o −m) 3 do i = 1 , ndim
4 Jz ( i , i ) = kk 5 kk = kk − 1! s t a r t from m, s u b t r a c t 1 t i l l r e a c h −m 6 end do! n o t e t h a t h e r e m’ i s t h e same a s m
7 8 Jp = n u l! i n i t i a l i z e t h e J+ m a t r i x
9 kk = j ∗ one! here , kk i s ’ supposed ’ t o be m; s e e l i n e 15 10 do i = 1 , ndim
11 i f ( i +1 > ndim ) then! e l i m i n a t e warning and e t c. 12 e x i t
13 e l s e
1 Mainly due to the fact that I am a little rusty on Quantum Mechanics.
14 Jp ( i , i +1) = s q r t ( j ∗ ( j +1) − ( kk −1) ∗ ( ( kk −1)+1) ) 15! ( kk −1) cuz kk = m’ but s h o u l d be kk = m, m = m’ −1 , i = m’ ; Jp (m’ ,m)
16 end i f 17 kk = kk − 1
18 end do 19
20 Jm = n u l! i n i t i a l i z e t h e J− m a t r i x 21 kk = j ∗ one! here , kk i s ’ supposed ’ t o be m; s e e l i n e 25
22 do i = 1 , ndim 23 i f ( i > 1 ) then
24 Jm( i , i −1) = s q r t ( j ∗ ( j +1) − ( kk+1) ∗ ( ( kk+1)−1) ) 25! ( kk+1) cuz kk = m’ but s h o u l d be kk = m, m = m’+1 , i = m’ ; Jp (m’ ,m)
26 end i f 27 kk = kk − 1 28 end do
29 30 Jx = 0. 5 _dp ∗ ( Jp + Jm)! c o n s t r u c t Jx
31 Jy = −(0.5_dp∗ i i c ) ∗ ( Jp − Jm)! c o n s t r u c t Jy 32 J s = matmul ( Jx , Jx ) + matmul ( Jy , Jy ) + matmul ( Jz , Jz )! c o n s t r u c t J^
Comments have been added to each line for explanation. Needless to say that each variable/parameter
such as Jz and kk have been defined in the preamble of the code. Note that ndim is the dimension of the
matrix, which is 18 here; one is equivalent to 1 + 0i, nul is a complex zero, and iic is the same as 0 + 1i.
The equivalent code in Python is listed below.
1 Jz = np. z e r o s ( shape =(ndim , ndim ) , dtype=complex )
2 kk = j ∗ one 3 f o r i i n r a n g e ( ndim ) :
4 Jz [ i ] [ i ] = kk 5 kk −= 1 6
7 Jp = np. z e r o s ( shape =(ndim , ndim ) , dtype=complex ) 8 kk = j ∗ one
9 f o r i i n r a n g e ( ndim ) : 10 i f ( i +1 > ndim−1) :
11 b r e a k 12 e l s e :
13 Jp [ i ] [ i +1] = np. s q r t ( j ∗ ( j +1) − ( kk −1) ∗ ( ( kk −1)+1) ) 14 kk −= 1
15 16 Jm = np. z e r o s ( shape =(ndim , ndim ) , dtype=complex )
17 kk = j ∗ one 18 f o r i i n r a n g e ( ndim ) :
19 i f i > 0 : 20 Jm [ i ] [ i −1] = np. s q r t ( j ∗ ( j +1) − ( kk+1) ∗ ( ( kk+1)−1) ) 21 kk −= 1
22 23 Jx = 0. 5 ∗ ( Jp + Jm)
24 Jy = −(0.5∗ i i c ) ∗ ( Jp − Jm) 25 J s = [ [ matrixmul ( Jx , Jx ) [ i ] [ j ] + matrixmul ( Jy , Jy ) [ i ] [ j ] + matrixmul ( Jz , Jz ) [ i ] [ j ]
f o r j i n r a n g e ( ndim ) ] f o r i i n r a n g e ( ndim ) ]
Note that matrixmul is my own defined matrix multiplication function (see Appendix A, can use np.matmul
more efficiently). Anyhow, the generated matrices Jx, Jy , Jz , and J 2 are included below. 2
Jx =
2 The 18 × 18 matrices are far too large to have any one of them fully included in this paper. Also, their non-zero entries only include 3 significant figures to save space as well as better visuals.
Table 1: Eigenvectors corresponding to Jx, Jy and J
2 respectively.
wx wy ws
− 8. 5 − 8. 5 80. 75
− 7. 5 − 7. 5 80. 75
− 6. 5 − 6. 5 80. 75
. . .
vx =
vy =
− 0. 003 i − 0. 011 i − 0. 032 i +0. 720 i · · · − 0. 011 i − 0. 003 i
− 0. 011 − 0. 041 − 0. 102 +0. 192 · · · +0. 041 +0. 011
+0. 032 i − 0. 102 i +0. 210 i − 0. 321 i · · · +0. 101 i +0. 032 i
. . .
− 0. 011 i +0. 041 i − 0. 102 i − 0. 192 i · · · +0. 041 i − 0. 011 i
− 0. 003 +0. 011 − 0. 032 − 0. 072 · · · − 0. 011 +0. 003
vs =
It can possibly be very important to note that the Eigenvector matrix output for vs was not as ordered as
equation (12); i.e. the 1’s were not necessarily in the main diagonal. Although since the Eigenvectors ws
shown in Table 1 are all equal to 80.75, the order of the 18 Eigenvectors should not matter.
1.2 Matrix Exponential Comparison
To compare the exponential matrices one must first take the exponential of the matrix. Fortran has no built
in module to do this; so I wrote the follow function is to be contained in the main program.
1 f u n c t i o n matrexp (A, m) r e s u l t ( Aexp )
2 u s e numtype 3 u s e numkam! f o r t h e f a c t o r i a l f u n c t i o n
4 u s e s e t u p 5 i m p l i c i t none
6 7 i n t e g e r : : m, i
8 complex ( dp ) , d i m e n s i o n (m, m) , i n t e n t ( i n ) : : A 9 complex ( dp ) , d i m e n s i o n (m, m) : : X, u n i t a r y , Aexp
10 11 u n i t a r y = n u l! i n i t i a l i u n i t a r y m a t r i x
12 do i = 1 , m 13 u n i t a r y ( i , i ) = one
14 end do! u n i t a r y m a t r i x I_{mxm} has been g e n e r a t e d 15
16 Aexp = u n i t a r y + A! i n i t i a l z i e t h e exp m a t r i x 17 X = A! i n i t i a l i z e dummy term i n exp sum d e f i n i t i o n
18 do i = 2 , 30! s t a r t from 2 cuz 0 , 1 terms a l r e a d y i n Aexp 19 X = matmul (X, A)
20 Aexp = Aexp + X/ ( f a c t o r i a l ( i ) ) 21 end do
22 end f u n c t i o n matrexp
Note that setup is where the aforementioned concepts such as nul, one, and iic have been defined. numkam
is the module where the used factorial function has been defined. Anyhow, the matrexp function approx-
imates the matrix exponential using the definition of the exponential:
exp(A) ≡
∞ ∑
n=
A
n
n!
= I + A +
AA
AAA
Obviously no programming language has the ability to sum infinite terms; and since matrexp can only sum
up to 30 terms (see line 18), it is nothing but an approximation. Very surprisingly Python 3 can sum up to
170 terms, over four times that of Fortran 95. Nevertheless it is still just an approximation that breaks down
for larger terms. That being said, Python also has a predefined matrix exponential function from the linalg
section of the SciPy library. To import the function the following line must be included in the Python script:
from scipy.linalg import expm
Although expm is also an approximation, it is a very good one, as it accurately applies the Pade approxima-
tion.
Now that the matrix exponential function has been defined and imported for GFortran and Python 3 re-
spectively, I will continue with solving the problem of comparing exp(iJy δ) exp(iJxδ) exp(−iJy δ) exp(−iJxδ)
with exp
iJz δ
2
for δ = π/ 10 and δ = π/ 20. For the first term the matrexp function has been used in Fortran
and expm in Python. However, since Jz is a diagonal matrix no approximation is needed; the exponential of
the matrix can be found exactly by simply taking the exponential of each term.
1 Ex = matrexp ( i i c ∗ Jx ∗ d e l t a , ndim )! c o n s t r u c t 1 s t exp m a t r i x
2 Ey = matrexp ( i i c ∗ d e l t a ∗Jy , ndim )! c o n s t r u c t 2nd exp m a t r i x 3 Exm = matrexp(− i i c ∗ d e l t a ∗Jx , ndim )! c o n s t r u c t 3 rd exp m a t r i x
4 Eym = matrexp(− i i c ∗ d e l t a ∗Jy , ndim )! c o n s t r u c t 4 th exp m a t r i x 5 mule = matmul ( matmul ( Ey , Ex ) , matmul (Eym, Exm) )! m u l t i p l y i n g a l l t h e above
6 7 Ez = 0
8 do i = 1 , ndim 9 Ez ( i , i ) = exp ( i i c ∗ Jz ( i , i ) ∗ ( d e l t a ∗ ∗ 2 ) )! d i a g o n a l m a t r i x e a s y 10 end do
The above is the Fortran code; where,
mule = exp(iJy δ) exp(iJxδ) exp(−iJy δ) exp(−iJxδ), (14)
and,
Ez = exp
iJz δ
2
The Python version of this is listed below.
1 Ez = 1 j ∗ Jz ∗ d e l t a ∗∗ 2 2 f o r i i n r a n g e ( ndim ) : 3 Ez [ i ] [ i ] = np. exp ( Ez [ i ] [ i ] )
4 5 Ex = expm ( 1 j ∗ Jx ∗ d e l t a )
6 Ey = expm ( 1 j ∗ Jy ∗ d e l t a ) 7 Exm = expm(−1 j ∗ Jx ∗ d e l t a )
8 Eym = expm(−1 j ∗ Jy ∗ d e l t a ) 9 mule = matrixmul ( matrixmul ( Ey , Ex ) , matrixmul (Eym, Exm) )
To compare mule with Ez, I will subtract the two matrices using the following line of code:
Using (20) the following function is written to store the coefficients.
8
1 c o n t a i n s
2 f u n c t i o n p o l y ( n ) r e s u l t ( q ) 3 u s e numtype
4 i m p l i c i t none 5
6 r e a l ( dp ) : : p , t 7 i n t e g e r : : i , n , m
8 r e a l ( dp ) , d i m e n s i o n ( n ) : : q! dont c a r e f o r q ( n ) i n m a t r i x 9
10 m = n+ 11! p = 0! i n i t i a l i z e t h e p o l y n o m i a l
12 do i = 0 , (m−2)! i f want q ( n ) can t o ( n−1) 13 q ( i +1) = m − i 14! p = p + q ∗ ( t ∗∗ i )! r e s u l t o f p o l y n o m i a l
15 end do 16 end f u n c t i o n
Now that the coefficients have been stored, an algorithm cab be written to generate the companion matrix
(19). The algorithm is listed below.
1 q = p o l y ( n )! s t o r i n g t h e c o e f f i c i e n t s i n q 2 A = 0! i n i t i a l i z e t h e companion m a t r i x
3 do i = 1 , ( n ) 4 i f ( ( i +1) > ( n ) ) then
5 e x i t! i g n o r i n g ’ warning ’ + bad c o d i n g 6 e l s e 7 A( i +1 , i ) = 1
8 A( i , n ) = −q ( i ) 9 end i f
10 end do ; A( n , n ) = −q ( n )
Now that the companion matrix (19) has bee created the roots of the polynomial pn can be found. These
roots are the Eigenvalues of the companion matrix [3]. The companion matrix A is real; yet it can definitely
has complex Eigenvalues. So the dgeev subroutine from the Lapack library should be a viable choice to
compute the Eigenvalues.
call dgeev(’v’, ’v’, n, e, n, wr, wi, vl, ldvl, vr, ldvr, work, lwork, info)
Note that e is the same as A; i.e. e=A was set to avoid making changes to A. Moreover, similar to the
zheev subroutine used in section 1.1, dgeev takes in many parameters to compute the Eigenvalues and
Eigenvectors. wr stores the real part of the Eigenvalues and wi stored the imaginary part. For n = 20, 25,
and 30 the Eigenvalues of the matrix (19) and roots of the polynomial (18) are in the tables below.
Table 2: Roots of pn(t)/Eigenvalues of A for various values of n
n = 20 n = 25 n = 30
04 ± 0. 37 i 1. 04 ± 0. 30 i 1. 04 ± 0. 25 i
90 ± 0. 68 i 0. 96 ± 0. 55 i 0. 98 ± 0. 46 i
68 ± 0. 93 i 0. 81 ± 0. 77 i 0. 88 ± 0. 66 i
. . .
− 0. 61 ± 1. 01 i − 0. 88 ± 0. 76 i − 0. 87 ± 0. 73 i
− 0. 88 ± 0. 79 i − 0. 67 ± 0. 94 i − 0. 71 ± 0. 89 i
Where for n = 20 the polynomial has 20 roots, n = 25 has 25 roots, and n = 30 has 30 roots.
8 Unlike the previous contained function matrexp it probably does not matter too much whether poly is contained in the main program or is an external function.
3 Discussion
For the first problem in section 1 the matrices which represent the angular momentum J have been generated
using both Fortran and Python. Then, the eigenvalues and Eigenvectors were computed with the help of the
Lapack library and its subroutine zheev. After that the exponential matrices mule and Ez were compared
for angles δ = π/ 20 and π/ 10. For the smaller angle both languages performed fine in computing them to
be equal. For the larger angle δ = π/ 10 only Python correctly computed Ez and mule to be the same. The
matrix exponential approximation matrexp written in GFortran was not sufficient enough to to calculate
mule (it did however almost approximate it be zero). Nevertheless, the correct results were computed and
included in section 1.2.
For the second problem in section 2 the roots of the monic polynomial pn were found for various values
of n. To do so the companion matrix A was constructed and it Eigenvalues were found using Lapack’s dgeev
and the results were included in Table 2. This problem was only solved in Fortran.
References
[1] Zoltan Papp, Mastering Computational Physics, Phys 562 Lecture Notes, CSU Long Beach, 2020.
[2] Todd Rowland, Companion Matrix, ‘MathWorld’–a Wolfram web recource, retrieved on 11/5/2020 from https:
//mathworld.wolfram.com/CompanionMatrix.html.
[3] Eigenvalue-Polynomials, Massachusetts Institute of Technology, 2017.
A matrixmul
def matrixmul(X, Y):
result = np.zeros(shape=(ndim, ndim), dtype=complex)
iterate through rows of X
for i in range(len(X)):
iterate through columns of Y
for j in range(len(Y[0])):
iterate through rows of Y
for k in range(len(Y)):
result[i][j] += X[i][k] * Y[k][j]
return result
Compare with the NumPy library’s np.matmul, and np.dot. np.matmul is probably the best, but matrixmul
allows me to have more control as it is my own defined function.
