CC10 PRACTICAL FILE

QUANTUM MECHANICS
CU REGISTRATION No: 223-1111-0336-22
CU ROLL NUMBER: 223223-21-0048
SEMESTER: 4

1
Contents
1 Problems and Solutions: 3
1.1 Problem 1: . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3
1.1.1 Source Code: . . . . . . . . . . . . . . . . . . . . . . . 3
1.1.2 Plots: . . . . . . . . . . . . . . . . . . . . . . . . . . . 5
1.2 Problem 2: . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6
1.2.1 Source Code: . . . . . . . . . . . . . . . . . . . . . . . 6
1.2.2 Plot: . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8
1.3 Problem 3: . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9
1.3.1 Source Code: . . . . . . . . . . . . . . . . . . . . . . . 9
1.3.2 Plot: . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12
1.4 Problem 4: . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13
1.4.1 Source Code: . . . . . . . . . . . . . . . . . . . . . . . 13
1.4.2 Plots: . . . . . . . . . . . . . . . . . . . . . . . . . . . 16
1.5 Problem 5: . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17
1.5.1 Source Code: . . . . . . . . . . . . . . . . . . . . . . . 17
1.5.2 Plot: . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19
1.6 Problem 6: . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20
1.6.1 Source Code: . . . . . . . . . . . . . . . . . . . . . . . 20
1.6.2 Plot: . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22

2
 1 Problems and Solutions:
1.1 Problem 1:

1.1.1 Source Code:

1 import numpy as np
2 import matplotlib . pyplot as plt
3 from scipy . optimize import newton
4 from scipy . integrate import odeint , simps
5

6 # Transendental Functions
7 f1 = lambda v : np . sqrt (( u0 ) **2 - v **2)
8 f2 = lambda v : -v / np . tan ( v )
9
10 # Parameters
11 a , V0 , m_by_hbar2 = 1.0 ,50.0 ,1.0
12 u0 = np . sqrt ( m_by_hbar2 * a **2* V0 /2)
13
14 # Scaling
15 v = np . linspace (1 e -6 , u0 ,1000 , endpoint = False )
16
17 # Plotting
18 plt . plot (v , f1 ( v ) , ls = " -" , label = " $ \ sqrt { u_0 ^2 - v ^2} $ " )
19 plt . plot (v , f2 ( v ) , ls = " --" , label = "$ - v \ cot { v } $ " )
20 plt . xlim ([0 , u0 ])
21 plt . ylim ([ -10 ,10])

3
22 plt . grid ()
23 plt . xlabel ( " $v$ " )
24 plt . legend ()
25 plt . savefig ( " Transendental_2 . png " )
26 plt . show ()
27
28 Guess = eval ( input ( " Enter a guess root from the graph = " ) )
29
30 # Transendental solution function
31 f = lambda v : np . sqrt ( u0 **2 - v **2) + v / np . tan ( v )
32
33 root = newton (f , Guess )
34

35 E =(2* root / a ) **2/(2* m_by_hbar2 )

36 print ( " Eigen - energy = " ,E )
37
38 # Potential
39 V = lambda x : V0 if (x <0 or x > a ) else 0.0
40

41 # Schrodinger solution
42 def SE ( ini ,x , E ) :
43 psi , psidot = ini
44 dpsi2_dx2 =(2.0* m_by_hbar2 ) *( V ( x ) -E ) * psi
45 return [ psidot , dpsi2_dx2 ]
46

47 b=a
48 x = np . linspace ( -b , a +b ,1000)
49 psi = odeint ( SE ,[0 ,1] , x , args =( E ,) ) [: ,0]
50
51 # Normalize
52 psi = psi / np . sqrt ( simps ( psi **2 , x ) )
53 prob = psi **2
54
55 plt . plot (x , psi , ’ - ’ , label = " $ \ psi ( x ) $ " )
56 plt . plot (x , prob , ’ -- ’ , label = " $ |\ psi ( x ) |^2 $ " )
57 plt . grid ()
58 plt . xlabel ( " x " )
59 plt . legend ( loc = ’ best ’)
60 plt . savefig ( " Transdendental . png " )
61 plt . show ()
62
63
64 ’’’ OUTPUT :
65
66 Enter a guess root from the graph = 2.59

4
67 Eigen - energy = 1 3. 47 57 22 73 92 42 37 3 ’’’

1.1.2 Plots:

Figure 1: Solution to Transcendental equation

Figure 2: Plot of the Wavefunction and corresponding probability

5
 1.2 Problem 2:

1.2.1 Source Code:

1 from scipy import optimize as opt

2 import numpy as np
3 import matplotlib . pyplot as plt
4 from scipy . integrate import simps
5
6 # Parameters
7 a , V0 , m_by_hbar2 = 1 ,40 ,1
8
9 # Euler method
10 def Solver (f , z0 ,x ,E , dx ) :
11 zs =[ z0 ]
12 z = z0
13 for i in range ( len ( x ) -1) :
14 z = z + dx * np . array ( f (z , x [ i ] , E ) )
15 zs . append ( z )
16 return np . array ( zs )
17
18 # Schrodinger Equation
19 def SE ( psi_psidot ,x , E ) :
20 psi , psidot = psi_psidot
21 psiddot =2.0* m_by_hbar2 *( V ( x ) -E ) * psi
22 return [ psidot , psiddot ]
23
24 # Potential
25 V = lambda x : np . where (( x >=0) & (x <= a ) ,-V0 , 0.0)
26
27 # Shooting algorithm

6
28 def shoot ( E ) :
29 psi = Solver ( SE , psi0 ,x ,E , dx ) [: ,0]
30 return psi [ -1]
31
32 b=a
33 x , dx = np . linspace ( -b , a +b ,1000 , retstep = True )
34 psi0 = np . array ([0 ,1 e -6])
35
36 # Ground State
37 Emin , Emax , dE =0 ,100 ,0.1
38 energies = np . arange ( Emin , Emax , dE )
39 shoots =[ shoot ( E ) for E in energies ]
40 Ev_guess = energies [ np . where ( np . diff ( np . signbit ( shoots ) ) ) ]
41 Ev = [ opt . newton ( shoot , i ) for i in Ev_guess [0:2]]
42
43 print ( " First two bound state energies = " , Ev )
44
45 def wavefunction ( E ) :
46 psi = Solver ( SE , psi0 ,x ,E , dx ) [: ,0]
47 psi = psi / np . sqrt ( simps ( psi * psi , x ) )
48 return psi
49
50 psi1 = wavefunction ( Ev [0])
51 psi2 = wavefunction ( Ev [1])
52

53 plt . plot (x , psi1 , ls = " -" , label = " First Bound state " )
54 plt . plot (x , psi2 , ls = " --" , label = " Second Bound state " )
55 plt . plot (x , V ( x ) /20 , ls = " -" , lw =0.8 , color = " r " , alpha =0.6 , label = "
Potential ( Scaled ) " )
56 plt . title ( " First two bound states of Finite potential well " )
57 plt . legend ( loc = " best " )
58 plt . xlabel ( " x " )
59 plt . grid ()
60 plt . savefig ( " finite_well . png " )
61 plt . show ()
62
63 ’’’ OUTPUT :
64
65 First two bound state energies = [2.0872570796673386 ,
4. 90 91 78 40 04 60 91 2] ’’’

7
1.2.2 Plot:

8
 1.3 Problem 3:

1.3.1 Source Code:

1 import numpy as np
2 import scipy as sp
3 import matplotlib . pyplot as plt
4
5 hb , m =0.1 ,1
6 def Solver (V ,x , psi0 , dpsi ,E , dx ) :
7 mh =2* m / hb **2
8
9 # Numerov
10 psi1 = psi0 + dpsi * dx
11 psi = np . array ([ psi0 , psi1 ])
12
13 for i in range (2 , len ( x ) ) :
14 a = 2*(1+ mh *5/12* dx **2*( V ( x [i -1]) -E ) )
15 b = -(1+ mh *1/12* dx **2*( V ( x [i -2]) -E ) )
16 d = (1+ mh *1/12* dx **2*( V ( x [ i ]) -E ) )
17
18 ppsi = a / d * psi [i -1] + b / d * psi [i -2]
19 psi = np . append ( psi , ppsi )
20
21 return psi
22
23
24
25 # Bisection
26 def bisect_TISE (V , Emin , Emax ,x , dx , nodes , tol =1 e -6) :
27 mxItr =2000 # Maximum number of iterations
28 N = len ( x )
29 psi0 , psiN =0 ,0
30 dpsi =( -1) ** nodes *1 e -3 # 2 nd initial value is set accordingly
31

32 i =0 # Iteration
33
34 while abs ( Emax - Emin ) > tol or i < mxItr :
35 E =( Emax + Emin ) /2

9
36 psi = Solver (V ,x , psi0 , dpsi ,E , dx )
37 C_nodes =0
38
39 for k in range (1 ,N -2) :
40 if psi [ k ]* psi [ k +1] <0:
41 C_nodes +=1
42

43 if C_nodes > nodes : # Counting nodes

44 Emax = E
45 elif C_nodes < nodes :
46 Emin = E
47 else :
48 if psi [N -1] > psiN : # Matching RHS boundary condition
49 Emin = E
50 else :
51 Emax = E
52
53 i +=1
54

55 return E , psi
56
57
58
59 # Potential
60 def V ( xx ) :
61 w =0.25
62 k = m * w **2
63 Vv = 0.5* k * xx **2
64 return Vv
65
66

67
68 # Main
69 x0 , xN = -2.0 ,2.0
70 x , dx = np . linspace ( x0 , xN ,1000 , retstep = True )
71 Vx = V ( x )
72

73 Emin , Emax = min ( Vx ) , max ( Vx )

74 node0 =0 # Number of nodes
75 node1 =1
76
77 # Bisection Calling
78 E0 , psi_0 = bisect_TISE (V , Emin , Emax ,x , dx , node0 )
79 E1 , psi_1 = bisect_TISE (V , Emin , Emax ,x , dx , node1 )
80

10
 81

82 # Normalization
83 from scipy . integrate import simpson
84 def norm (x , psi ) :
85 psimod2 = psi **2
86 inte = simpson ( psimod2 , x = x )
87 Psi = psi / inte **0.5
88 return Psi
89
90 psi0 = norm (x , psi_0 ) # Normalized psi_0
91 psi1 = norm (x , psi_1 ) # Normalized psi_1
92
93 # Output and plotting
94 print ( f " Eigenvalues of Ground and 1 st excited state = { round (
E0 ,5) } ,{ round ( E1 ,5) } " )
95 plt . title ( " The Harmonic Oscillator Eigenfunctions for $n =0 ,1 $
")
96 plt . plot (x , psi_0 , ls = " -" , label = f " Excited state { node0 } " )
97 plt . plot (x , psi_1 , ls = " --" , label = f " Excited state { node1 } " )
98 plt . grid ()
99 plt . axhline ( y =0 , color = " k " )
100 plt . axvline ( x =0 , color = " k " )
101 plt . legend ( loc = " best " )
102 plt . savefig ( " H a rm o n ic _ o sc i l la t o r . png " )
103 plt . show ()
104
105
106
107 ’’’ OUTPUT :
108
109 Eigenvalues of Ground and 1 st excited state = 0.01534 ,0.04629
110
111 ’’’

11
1.3.2 Plot:

12
 1.4 Problem 4:

1.4.1 Source Code:

1 from scipy import optimize as opt

2 import numpy as np
3 from scipy . integrate import simps
4 import matplotlib . pyplot as plt
5
6 # Euler method
7 def Solver (f , z0 ,s , En ) :
8 zs =[ z0 ]
9 z = z0
10 for i in range ( len ( s ) -1) :
11 z = z +( s [ i +1] - s [ i ]) * np . array ( f (z , s [ i ] , En ) )
12 zs . append ( z )
13 return np . array ( zs )
14
15
16 def SEr (y ,r , En ) :
17 (u , up ) = y
18 upp = ( l *( l +1) / r **2 -2/ r - En ) * u
19 return np . array ([ up , upp ])
20
21 # Parameters
22 n , l =1 ,0 # 1 s state

13
23 En = -1.0/ n **2
24
25 def SolveSEr ( En ) :
26 rb = r [:: -1]
27 u0 = np . array ([0.0 , -1 e -3])
28 ub = Solver ( SEr , u0 , rb , En )
29 u = ub [: ,0][:: -1]
30 u = u / np . sqrt ( simps ( u **2 , r ) )
31 return u
32 r = np . logspace ( -4.0 ,2 ,1000)
33 u = SolveSEr ( En )
34
35 # Actual Radial Solution
36 R=u/r
37
38 plt . plot (r ,R , " -k " , label = " Radial Wavefunction $R ( r ) $ " )
39 plt . xlabel ( ’ Radial Distance ( r ) in Bohr Radius ( $a_0$ ) ’)
40 plt . title ( " Plot of Radial Wavefunction " )
41 plt . ylabel ( " R ( r ) " )
42 plt . grid ()
43 plt . xlim (0 ,6)
44 plt . ylim (0 ,2.5)
45 plt . legend ()
46 plt . savefig ( " radial_wf . png " )
47 plt . show ()
48
49 # Probability density
50 PD = u **2
51
52 plt . plot (r , PD , " -k " , label = " Probability Density $P ( r ) $ " )
53 plt . xlabel ( ’ Radial Distance ( r ) in Bohr Radius ( $a_0$ ) ’)
54 plt . title ( " Plot of Probability Density " )
55 plt . xlim (0 ,5)
56 plt . ylim (0 ,1)
57 plt . grid ()
58 plt . legend ()
59 plt . savefig ( " prob_density . png " )
60 plt . show ()
61
62
63 # Computational Average Value
64 r_avg_comp = simps ( u * r *u , r )
65

66 # Computational Most Probable value

67 r_mp_comp = r [ np . argmax ( PD ) ]

14
 68

69 print ( " COMPUTATIONAL : " )

70 print ( f " Average distance = { r_avg_comp } a_0 " )
71 print ( f " Most probable distance = { r_mp_comp } a_0 \ n \ n " )
72
73 # Theoritical Values
74 r_avg_th = 1.0
75 r_mp_th = 1.5
76
77 print ( " THEORETICAL : " )
78 print ( f " Theoretical Average distance = { r_avg_th } a_0 " )
79 print ( f " Theoretical Most probable distance = { r_mp_th } a_0 \ n \ n
")
80
81 print ( " ERRORS : " )
82 print ( f " Average distance = { round ( abs ( r_avg_th - r_avg_comp ) ,4) }
a_0 " )
83 print ( f " Most probable distance = { round ( abs ( r_mp_th - r_mp_comp )
,4) } a_0 \ n " )
84
85 print ( " Where a_0 is Bohr Radius " )
86
87
88 ’’’ OUTPUT :
89

90 COMPUTATIONAL :
91 Average distance = 1. 48 970 00 35 65 16 83 9 a_0
92 Most probable distance = 0.9 86 26 58 46 13 12 82 1 a_0
93
94
95 THEORETICAL :
96 Theoretical Average distance = 1.0 a_0
97 Theoretical Most probable distance = 1.5 a_0
98
99
100 ERRORS :
101 Average distance = 0.4897 a_0
102 Most probable distance = 0.5137 a_0
103
104 Where a_0 is Bohr Radius ’’’

15
1.4.2 Plots:

16
 1.5 Problem 5:

1.5.1 Source Code:

1 from scipy import optimize as opt

2 import numpy as np
3 import matplotlib . pyplot as plt
4 from scipy . integrate import simps
5
6 # Parameters
7 a = 10
8
9 # Euler method
10 def Solver (f , z0 ,x ,E , dx ) :
11 zs =[ z0 ]
12 z = z0
13 for i in range ( len ( x ) -1) :
14 z = z + dx * np . array ( f (z , x [ i ] , E ) )
15 zs . append ( z )
16 return np . array ( zs )
17
18 # Schrodinger Equation
19 def SE ( psi_psidot ,x , E ) :
20 psi , psidot = psi_psidot
21 psiddot =( V ( x ) -E ) * psi
22 return [ psidot , psiddot ]
23
24 # Potential
25 V = lambda x : x
26
27 # Shooting algorithm
28 def shoot ( E ) :

17
29 psi = Solver ( SE , psi0 ,x ,E , dx ) [: ,0]
30 return psi [ -1]
31
32 x , dx = np . linspace (0 ,a ,1000 , retstep = True )
33 psi0 = np . array ([0 ,1 e -6])
34
35 Emin , Emax , dE =0 ,100 ,0.1
36 energies = np . arange ( Emin , Emax , dE )
37
38 # Shoot adjustment
39 shoots =[ shoot ( E ) for E in energies ]
40 Ev_guess = energies [ np . where ( np . diff ( np . signbit ( shoots ) ) ) ]
41 Ev = [ opt . newton ( shoot , i ) for i in Ev_guess [0:2]]
42
43 print ( " First two bound state energies = " , Ev )
44
45 def wavefunction ( E ) :
46 psi = Solver ( SE , psi0 ,x ,E , dx ) [: ,0]
47 psi = psi / np . sqrt ( simps ( psi * psi , x ) )
48 return psi
49
50 psi1 = wavefunction ( Ev [0])
51 psi2 = wavefunction ( Ev [1])
52
53 plt . plot (x , psi1 + Ev [0] , ls = " -" , label = " First Bound state " )
54 plt . plot (x , psi2 + Ev [1] , ls = " --" , label = " Second Bound state " )
55 plt . plot (x , V ( x ) , ’ - ’ , lw =0.7 , alpha =0.6 , color = " r " , label = "
Potential " )
56 plt . title ( " First two bound states of Triangular potential
well " )
57 plt . axhline ( y = Ev [0] , color = " k " , lw =0.6)
58 plt . axhline ( y = Ev [1] , color = " k " , lw =0.6)
59 plt . legend ( loc = " best " )
60 plt . ylabel ( " $ \ psi ( x ) $ " )
61 plt . xlabel ( " $x$ " )
62 plt . savefig ( " Triangular_well . png " )
63 plt . show ()
64
65
66
67 ’’’ OUTPUT :
68
69 First two bound state energies = [2.3331754390428237 ,
4. 08 31 67 69 24 07 03 8]
70 ’’’

18
1.5.2 Plot:

Figure 3: Horizontal black lines represent the eigenvalues of the correspond-

ing states

19
 1.6 Problem 6:

1.6.1 Source Code:

1 import numpy as np
2 import scipy as sp
3 import matplotlib . pyplot as plt
4

5 l =40
6 dx =0.01
7 ti , tf , dt =0 ,10 ,0.1
8 N =1000
9 x , dx = np . linspace ( - l /2 , l /2 ,N , retstep = True )
10

11 t = np . arange ( ti , tf + dt , dt )
12
13 x0 , sig , k = -10.0 ,2.0 , np . pi /10
14
15 gwp = 1/((2* np . pi ) **0.5* sig ) **0.5* np . exp ( -(x - x0 ) **2/(4* sig
**2) +1 j * k * x )
16
17 line , = plt . plot (x , np . abs ( gwp **2) )
18
19 # Potential
20 V =0* x [1: -1]
21

22 alpha = 1 j * dt /(2* dx **2)

23 beta = 1+2* alpha +(1 j * dt * V ) *0.5
24 gamma = 1 -2* alpha -(1 j * dt * V ) *0.5
25

20
26 # Defining transformation matrices
27 U1 = - alpha * np . eye (N -2 , k = -1) + np . diag ( beta ) - alpha * np . eye (N -2 , k
=1)
28 U2 = alpha * np . eye (N -2 , k = -1) + np . diag ( gamma ) + alpha * np . eye (N -2 , k
=1)
29
30 U1_inv = np . linalg . inv ( U1 )
31
32 # Transformation matrix P
33 P = np . dot ( U1_inv , U2 )
34
35 plt . xlabel ( " x ( a . u .) " )
36 plt . ylabel ( " $ |\ psi ( x ) |^{2} $ " )
37 plt . title ( " Time evolution of a Gaussian wave - packet " )
38 tt =0
39 T_inst =[0 ,2 ,4]
40
41 for i , T in enumerate ( T_inst ) :
42 while tt < T :
43 gwp [1: -1]= np . dot (P , gwp [1: -1])
44 tt += dt
45 plt . plot (x , np . abs ( gwp ) **2 , ls =(0 ,(3 , i ) ) , label = f ’t ={ T_inst [ i
]} ( a . u .) ’)
46 plt . legend ( loc = ’ best ’)
47 plt . grid ()
48 plt . savefig ( " G a us s i an _ w av e p ac k e t . png " )
49 plt . show ()

21
1.6.2 Plot:

22