Python Assignment 3
June 15, 2026
Name: Tapasmit Naskar
Dept: Physics
Roll No. 002420701014
Class: UG II (Sem II)
Applications of numpy/scipy modules
1. Projectile Motion Problem (Root Finding): A ball is launched with an initial speed v = 50 m/s
at an angle θ = 45◦ from the ground. Find the time when the ball hits the ground. Equation
to Solve: y(t) = vt sin(θ)0.5gt2 = 0 with g = 9.8m/s2
[1]: import numpy as np
from [Link] import fsolve
def eq(t):
v = 50
th = np.deg2rad(45)
g = 9.8
return v*t*[Link](th) - 0.5*g*t**2
time = fsolve(eq, x0=100)
print("Time of hitting ground (in sec): ",time[0])
Time of hitting ground (in sec): 7.215375318230077
2. Radioactive Decay Problem (ODE Solving): A radioactive substance A decays into a stable
substance B with a decay constant λ = 0.5 per day. If initially there are 1000 atoms of A,
nd the number of atoms of A after 1, 5 and 10 days.
[2]: from [Link] import solve_ivp
N = 1000
def decay(t, n):
l = 0.5
return -l*n
sol1 = solve_ivp(decay, t_span=(0,10), y0=[1000], t_eval=[Link](0,1,100))
sol2 = solve_ivp(decay, t_span=(0,10), y0=[1000], t_eval=[Link](0,5,500))
sol3 = solve_ivp(decay, t_span=(0,10), y0=[1000], t_eval=[Link](0,10,1000))
print("Approximate no. of atoms of A")
print("after 1 day,", round(sol1.y[0][-1]))
print("after 5 days,", round(sol2.y[0][-1]))
1
print("after 10 days,", round(sol3.y[0][-1]))
Approximate no. of atoms of A
after 1 day, 607
after 5 days, 82
after 10 days, 7
3. Integration Problem: A particle moves along a straight line under a variable force F (x) =
5x2 N . Calculate the work done by the force as the particle moves from x = 0 to x = 4 meters.
Hint: Work done
[3]: from sympy import symbols, integrate
from [Link] import quad
def force(x):
return 5*x**2
#symbolic integration
x = symbols('x')
res1 = integrate(force(x), (x,0,4))
#numerical integration
res2, error = quad(force,0,4)
print("Work done is", res1,"or",res2,"joules")
Work done is 320/3 or 106.66666666666666 joules
4. Compute and plot the spectral radiance of a blackbody as a function of wavelength at dif-
ferent temperatures (3000, 5000, 10000)K using Planck's law (use numpy modules such as
[Link], [Link], [Link] etc)
[4]: import [Link] as plt
def wav_dist(l,T):
h = 6.626*10**(-34)
c = 2.998*10**8
k_B = 1.381*10**(-23)
return (2*h*c*c/(l**5))/([Link](h*c/(l*k_B*T)) - 1)
l_list = [Link](0,5*10**(-6),10000)
radiance1 = wav_dist(l_list,3000)
radiance2 = wav_dist(l_list,5000)
radiance3 = wav_dist(l_list,10000)
fig, axes = [Link](1, 3, figsize=(15,4))
axes[0].plot(l_list, radiance1, color='red')
axes[0].set_title("at 3000 K")
axes[0].set_xlabel("wavelength")
axes[0].set_ylabel("Spectral radiance")
2
axes[1].plot(l_list, radiance2, color='green')
axes[1].set_title("at 5000 K")
axes[1].set_xlabel("wavelength")
axes[1].set_ylabel("Spectral radiance")
axes[2].plot(l_list, radiance3, color='blue')
axes[2].set_title("at 10000 K")
axes[2].set_xlabel("wavelength")
axes[2].set_ylabel("Spectral radiance")
[Link]("Spectral radiance vs Wavelength")
[Link]()
[Link]()
C:\Users\tapas\AppData\Local\Temp\ipykernel_31124\[Link]:
RuntimeWarning: divide by zero encountered in divide
return (2*h*c*c/(l**5))/([Link](h*c/(l*k_B*T)) - 1)
C:\Users\tapas\AppData\Local\Temp\ipykernel_31124\[Link]:
RuntimeWarning: overflow encountered in exp
return (2*h*c*c/(l**5))/([Link](h*c/(l*k_B*T)) - 1)
C:\Users\tapas\AppData\Local\Temp\ipykernel_31124\[Link]:
RuntimeWarning: invalid value encountered in divide
return (2*h*c*c/(l**5))/([Link](h*c/(l*k_B*T)) - 1)
C:\Users\tapas\AppData\Local\Temp\ipykernel_31124\[Link]: UserWarning:
No artists with labels found to put in legend. Note that artists whose label
start with an underscore are ignored when legend() is called with no argument.
[Link]()
5. Use the scipy module [Link].curve_ t to t a data with a given curve. [Choose your
own data set and the function]
[5]: from [Link] import curve_fit
t_dat, temp = [Link]("[Link]",unpack=True)
3
def coolmodel(t,T0,k):
T_env=20
return T_env + (T0 - T_env)*[Link](-k*t)
popt, pcov = curve_fit(coolmodel,t_dat,temp,p0=(80,0.05))
dense_t = [Link](0,50,2000)
temp_dense = coolmodel(dense_t,*popt)
print(popt)
[Link]("Cooling Curve")
[Link]("Time")
[Link]("Temperature")
[Link](t_dat,temp,c='r',label="real data")
[Link](dense_t,temp_dense,label="fitted curve")
[Link]()
[Link]()
[8.44182403e+01 3.69683679e-02]
[ ]: