DSC – 11 (PRACTICAL)
PYTHON PROGRAMMING
1|Page
[Link] TOPIC NAME DATE TEACHER
SIGNATURE
1. Eigen Values 18-11-25
2. Shooting Method 02-12-25
3.
4.
5.
6.
CONTENTS
2|Page
1. FINDING EIGEN VALUES.
QUESTION 1.–
CODE-
import [Link] as plt
import numpy as np
from [Link] import newton
from [Link] import odeint, simps
V0 = 20
m_h_cut2 = 1
a=1
z0 = a * [Link](2 * m_h_cut2 * V0)
z = [Link](0.01, z0, 100)
f1 = lambda z: [Link]((z0 / z) ** 2 - 1)
3|Page
f2 = lambda z: -1 / [Link](z)
f3 = lambda z: [Link]((z0 / z) ** 2 - 1) + 1 / [Link](z)
[Link](figsize=(8, 5))
[Link](-5, 5)
[Link](z, f1(z), label=r'$f_1(z) = \sqrt{(z_0/z)^2 - 1}$', color='blue')
[Link](z, f2(z), label=r'$f_2(z) = -1/\tan(z)$', color='red', linestyle='--')
[Link]("PLOT 1: Transcendental Functions", fontsize=13)
[Link]("z", fontsize=12)
[Link]("f(z)", fontsize=12)
[Link]()
[Link](alpha=0.4)
[Link]()
root = newton(f3, 2.7)
root1 = newton(f3, 5.2)
E = lambda z: (z ** 2 / (2 * m_h_cut2 * a ** 2)) - V0
r = [Link]([root, root1])
eig = E(r)
print(f"Root 1 = {root:.4f}")
print(f"Root 2 = {root1:.4f}")
print(f"Eigenvalues (E1, E2) = {eig[0]:.4f}, {eig[1]:.4f}")
def V(x):
return -V0 if abs(x) <= a else 0
x = [Link](-5, 5, 100)
def se(m, x, E):
psi, dpsi = m
d2psi = -2 * m_h_cut2 * (E - V(x)) * psi
4|Page
return [dpsi, d2psi]
e = eig[0]
sol = odeint(se, [0, 1], x, args=(e,))[:, 0]
norm_psi = sol / [Link](simps(sol ** 2, x))
[Link](figsize=(8, 5))
[Link](x, norm_psi, color='purple')
[Link]("PLOT 2: Normalized Wavefunction ψ(x)", fontsize=13)
[Link]("x", fontsize=12)
[Link](r'$\psi(x)$', fontsize=12)
[Link](alpha=0.4)
[Link]()
OUTPUT-
5|Page
6|Page
QUESTION 2. –
CODE-
import [Link] as plt
import numpy as np
from [Link] import newton
from [Link] import odeint, simps
V0 = 20
m_h_cut2 = 1
a=2
u0 = [Link](m_h_cut2 * a**2 * V0 / 2)
v = [Link](0.01, u0, 100)
7|Page
f1 = lambda v: [Link](u0**2 - v**2)
f2 = lambda v: v * [Link](v)
f3 = lambda v: [Link](u0**2 - v**2) - v * [Link](v)
[Link](figsize=(8,5))
[Link](v, f1(v), label=r'$\sqrt{u_0^2 - v^2}$', color='blue')
[Link](v, f2(v), label=r'$v tan(v)$', color='red', linestyle='--')
[Link]("PLOT 1: Transcendental Functions (Even Parity)", fontsize=13)
[Link]("v")
[Link]("f(v)")
[Link]()
[Link](alpha=0.4)
root1 = newton(f3,1.308 )
root2 = newton(f3,3.952 )
E = lambda v: (2 * v**2) / a**2
eig = [Link]([E(root1), E(root2)])
print(f"Root 1 = {root1:.4f}")
print(f"Root 2 = {root2:.4f}")
print(f"Eigenvalues (E1, E2) = {eig[0]:.4f}, {eig[1]:.4f}")
def V(x):
return 0 if abs(x) <= a/2 else V0
x = [Link](-5, 5, 500)
8|Page
def se(m, x, E):
psi, dpsi = m
d2psi = -2 * m_h_cut2 * (E - V(x)) * psi
return [dpsi, d2psi]
for e in eig:
sol = odeint(se, [1.0, 0.0], x, args=(e,))[:,0]
norm_psi = sol / [Link](simps(sol**2, x))
[Link](figsize=(8,5))
[Link](x, norm_psi, color='purple')
[Link](f"PLOT 2: Normalized Wavefunction ψ(x) for E = {e:.3f}")
[Link]("x")
[Link]("ψ(x)")
[Link](alpha=0.4)
[Link]()
Output-
9|Page
10 | P a g e
2. SHOOTING METHOD
QUESTION 1-
CODE-
from [Link] import newton
from [Link] import odeint
import numpy as np
import [Link] as plt
x = [Link](-5, 5, 100)
m = h2 = 1
def V(x):
if abs(x) < 0:
return 100
else:
return x
def se(psi, x, E):
psi, dpsi = psi
d2psi = -h2 * (E - V(x)) * psi
11 | P a g e
return [dpsi, d2psi]
E = [Link](-5, 0, 100)
def shoot(e):
sol = odeint(se, [0, 1], x, args=(e,))
return sol[-1, 0]
shoots = [Link]([shoot(e) for e in E])
root1 = newton(shoot, -0.9)
root2 = newton(shoot, -2.6)
print("Roots are :", root1, root2)
[Link](E, shoots)
[Link]("Energy (E)")
[Link]("ψ(x=L)")
[Link](-1, 1)
[Link]()
[Link]()
OUTPUT-
12 | P a g e
13 | P a g e