I’ll create a simulation of a Spring-Mass System that demonstrates both simple harmonic motion and damped oscillation with external forcing.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82

import numpy as np
from scipy.integrate import odeint
import matplotlib.pyplot as plt

def spring_mass_system(state, t, m, k, c, F0, omega_f):
    """
    Define the system of first-order ODEs for the spring-mass system
    state[0] = x (position)
    state[1] = v (velocity)
    """
    x, v = state
    
    # External forcing function
    F_external = F0 * np.sin(omega_f * t)
    
    # Acceleration from Newton's second law
    acceleration = (-k * x - c * v + F_external) / m
    
    return [v, acceleration]

# System parameters
m = 1.0  # mass (kg)
k = 10.0  # spring constant (N/m)
c = 0.5   # damping coefficient (Ns/m)
F0 = 1.0  # forcing amplitude (N)
omega_n = np.sqrt(k/m)  # natural frequency
omega_f = omega_n  # forcing frequency (at resonance)

# Initial conditions
x0 = 1.0  # initial displacement (m)
v0 = 0.0  # initial velocity (m/s)
state0 = [x0, v0]

# Time array
t = np.linspace(0, 20, 1000)

# Solve ODEs for three cases
# 1. Simple harmonic motion (no damping, no forcing)
solution_shm = odeint(spring_mass_system, state0, t, args=(m, k, 0, 0, omega_f))

# 2. Damped oscillation (with damping, no forcing)
solution_damped = odeint(spring_mass_system, state0, t, args=(m, k, c, 0, omega_f))

# 3. Forced oscillation (with damping and forcing)
solution_forced = odeint(spring_mass_system, state0, t, args=(m, k, c, F0, omega_f))

# Create plots
fig, (ax1, ax2, ax3) = plt.subplots(3, 1, figsize=(12, 12))

# Plot 1: Simple Harmonic Motion
ax1.plot(t, solution_shm[:, 0], 'b-', label='Position')
ax1.plot(t, solution_shm[:, 1], 'r--', label='Velocity')
ax1.set_title('Simple Harmonic Motion')
ax1.set_xlabel('Time (s)')
ax1.set_ylabel('Position (m) / Velocity (m/s)')
ax1.grid(True)
ax1.legend()

# Plot 2: Damped Oscillation
ax2.plot(t, solution_damped[:, 0], 'b-', label='Position')
ax2.plot(t, solution_damped[:, 1], 'r--', label='Velocity')
ax2.set_title('Damped Oscillation')
ax2.set_xlabel('Time (s)')
ax2.set_ylabel('Position (m) / Velocity (m/s)')
ax2.grid(True)
ax2.legend()

# Plot 3: Forced Oscillation
ax3.plot(t, solution_forced[:, 0], 'b-', label='Position')
ax3.plot(t, solution_forced[:, 1], 'r--', label='Velocity')
ax3.set_title('Forced Oscillation with Damping')
ax3.set_xlabel('Time (s)')
ax3.set_ylabel('Position (m) / Velocity (m/s)')
ax3.grid(True)
ax3.legend()

# Calculate and display system parameters
print(f"Natural frequency: {omega_n:.2f} rad/s")
print(f"Damping ratio: {c/(2*np.sqrt(m*k)):.2f}")

plt.tight_layout()

Let me explain the key components of this simulation:

1. System Parameters:

   • Mass $(m) = 1.0 kg$
   • Spring constant $(k) = 10.0 N/m$
   • Damping coefficient $(c) = 0.5 Ns/m$
   • Natural frequency $(ωn) = √(k/m) ≈ 3.16 rad/s$
   • Initial displacement = $1.0 m$
   • Initial velocity = $0.0 m/s$

2. Three Cases Simulated:

   • Simple Harmonic Motion (no damping, no forcing)
   • Damped Oscillation (with damping, no forcing)
   • Forced Oscillation (with damping and external force)

3. Mathematical Model:

   • Based on Newton’s Second Law: F = ma
   • Differential equation: $mẍ + cẋ + kx = F₀sin(ωt)$
   • Solved using scipy’s odeint function

4. Plot Analysis:

   • Top Plot (Simple Harmonic Motion):

     • Shows perfect sinusoidal oscillation
     • Energy continuously transfers between potential and kinetic

   • Middle Plot (Damped Oscillation):

     • Amplitude decreases exponentially due to damping
     • System eventually comes to rest

   • Bottom Plot (Forced Oscillation):

     • Initial transient response followed by steady-state oscillation
     • Shows resonance effects as forcing frequency matches natural frequency

5. Physical Insights:

   • Energy conservation in simple harmonic motion
   • Energy dissipation in damped systems
   • Resonance phenomenon in forced oscillation
   • Phase relationships between position and velocity

Result

Natural frequency: 3.16 rad/s
Damping ratio: 0.08

Example

The Particle in a Box is a fundamental problem in quantum mechanics that describes a particle confined to a one-dimensional box with infinite potential walls.

The particle’s wavefunction and energy levels can be solved using the Schrödinger equation.

Problem Statement

Given a particle of mass $( m )$ confined in a $1D$ box of length $( L )$, calculate and visualize:

1. The wavefunctions $( \psi_n(x) )$ for the first three quantum states $( n = 1, 2, 3 )$.
2. The corresponding energy levels $( E_n )$.

Equations

1. Wavefunction:
   $$
   \psi_n(x) = \sqrt{\frac{2}{L}} \sin\left(\frac{n \pi x}{L}\right)
   $$

   • $( n )$: Quantum number $( n = 1, 2, 3, \ldots )$.
   • $( x )$: Position within the box $( 0 \leq x \leq L $)$.

2. Energy Levels:
   $$
   E_n = \frac{n^2 h^2}{8 m L^2}
   $$

   • $( h )$: Planck’s constant $( 6.626 \times 10^{-34} \text{Js} )$.
   • $( m )$: Mass of the particle.
   • $( L )$: Length of the box.

Parameters

• Box length: $( L = 1 , \text{nm} = 1 \times 10^{-9} \text{m} )$.
• Particle mass:$ ( m = 9.11 \times 10^{-31} \text{kg} )$ (mass of an electron).
• Planck’s constant: $( h = 6.626 \times 10^{-34} \text{Js} )$.

Python Code

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39

import numpy as np
import matplotlib.pyplot as plt

# Constants
L = 1e-9  # Length of the box in meters (1 nm)
m = 9.11e-31  # Mass of the particle (electron) in kg
h = 6.626e-34  # Planck's constant in Js

# Quantum numbers
n_values = [1, 2, 3]  # First three quantum states

# Position array
x = np.linspace(0, L, 1000)  # 1000 points within the box

# Energy levels
E_n = [(n**2 * h**2) / (8 * m * L**2) for n in n_values]

# Wavefunctions
wavefunctions = [np.sqrt(2 / L) * np.sin(n * np.pi * x / L) for n in n_values]

# Plotting the wavefunctions
plt.figure(figsize=(12, 6))
for i, psi in enumerate(wavefunctions):
    plt.plot(x * 1e9, psi, label=f"n={n_values[i]}, E={E_n[i]:.2e} J")
plt.title("Wavefunctions of a Particle in a 1D Box")
plt.xlabel("Position (nm)")
plt.ylabel("Wavefunction $\\psi_n(x)$")
plt.legend()
plt.grid()
plt.show()

# Energy level plot
plt.figure(figsize=(6, 4))
plt.bar(n_values, E_n, width=0.5, color='skyblue')
plt.title("Energy Levels of a Particle in a 1D Box")
plt.xlabel("Quantum Number (n)")
plt.ylabel("Energy (J)")
plt.grid(axis='y')
plt.show()

Code Explanation

1. Constants:

   • $ L $, $ m $, and $ h $ are defined in $SI$ units.

2. Energy Calculation:

   • The energy levels $ E_n $ are calculated for $( n = 1, 2, 3 )$ using the energy formula.

3. Wavefunction Calculation:

   • Wavefunctions are computed for each $ n $ using the sinusoidal formula.

4. Visualization:

   • The first graph shows the spatial wavefunctions $ \psi_n(x) $.
   • The second graph represents the energy levels.

Graph Explanation

1. Wavefunctions:
   • The wavefunctions have $ n-1 $ nodes (zero crossings) for quantum state $ n $.
   • They represent the probability amplitude of finding the particle at position $ x $.

2. Energy Levels:
   • Energy levels increase quadratically with $ n $, as $ E_n \propto n^2 $.

Results

1. The particle has distinct wavefunctions and quantized energy levels.
2. Higher quantum states have more complex wavefunctions with additional nodes.
3. Energy levels grow rapidly, reflecting the confinement in a small space.

Example

In the theory of relativity, time dilation refers to the difference in elapsed time as measured by two observers due to a relative velocity between them or the presence of a gravitational field.

The formula for time dilation due to relative velocity is:

$$
\Delta t’ = \frac{\Delta t}{\sqrt{1 - \frac{v^2}{c^2}}}
$$

• $ \Delta t $:
  Proper time (time measured in the rest frame of the object).
• $ \Delta t’ $:
  Dilated time (time measured by a moving observer).
• $ v $:
  Velocity of the moving object.
• $ c $:
  Speed of light $( 3 \times 10^8 \ \text{m/s} )$.

Problem Statement

A spaceship is moving at different fractions of the speed of light $( v/c )$.

Compute the time dilation for a proper time of $1$ second as observed in the spaceship’s rest frame.

Visualize the relationship between the velocity fraction $( v/c )$ and the dilated time.

Python Code

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24

import numpy as np
import matplotlib.pyplot as plt

# Constants
c = 3e8  # Speed of light in m/s
proper_time = 1  # Proper time in seconds

# Generate velocity fractions (v/c)
velocity_fractions = np.linspace(0, 0.99, 100)  # Avoid v/c = 1 to prevent division by zero
velocities = velocity_fractions * c

# Calculate dilated time for each velocity
dilated_times = proper_time / np.sqrt(1 - (velocities / c)**2)

# Plotting the results
plt.figure(figsize=(10, 6))
plt.plot(velocity_fractions, dilated_times, label="Dilated Time", color="blue")
plt.axhline(y=proper_time, color='red', linestyle='--', label="Proper Time (1s)")
plt.title("Time Dilation vs Velocity Fraction (v/c)")
plt.xlabel("Velocity Fraction (v/c)")
plt.ylabel("Dilated Time (s)")
plt.legend()
plt.grid()
plt.show()

Code Explanation

1. Constants:

   • The speed of light $( c )$ is set to $ 3 \times 10^8 \text{m/s} $.
   • The proper time $( \Delta t )$ is set to $1$ second.

2. Velocity Fractions:

   • We create an array of velocity fractions $( v/c )$ ranging from $0$ to $0.99$.
     The value $0.99$ is chosen to avoid the singularity at $( v = c )$.

3. Dilated Time Calculation:

   • Using the time dilation formula, the dilated time is calculated for each velocity fraction.

4. Visualization:

   • A graph is plotted with the velocity fraction $( v/c )$ on the $x$-axis and the dilated time $( \Delta t’ )$ on the $y$-axis.
   • The proper time ($1$ second) is shown as a reference line.

Graph Explanation

• $x$-axis: Represents the velocity fraction $( v/c )$.
• $y$-axis: Represents the dilated time $( \Delta t’ )$.
• As the velocity approaches the speed of light $( v/c \to 1 )$, the dilated time $( \Delta t’ )$ increases dramatically, illustrating the effect of time dilation.
• The red dashed line represents the proper time of $1$ second.
  For low velocities $( v/c \approx 0 )$, the dilated time is almost equal to the proper time.

Results

1. At $ v/c = 0.1 $, the dilated time is approximately $ 1.005 \text{s} $.
2. At $ v/c = 0.9 $, the dilated time is approximately $ 2.294 \text{s} $.
3. As $ v/c \to 1 $, the dilated time tends to infinity, showing the extreme effects of relativistic speeds.

Example Problem

We will calculate and visualize the electric field generated by a uniform line charge along the $ x $-axis.

Problem Description

1. Uniform Line Charge:

   • A line charge with linear charge density $( \lambda )$ (charge per unit length) extends from $( x = -L )$ to $( x = +L )$.
   • The total charge is $( Q = 2L \lambda )$.

2. Electric Field Calculation:

   • The electric field at a point $ P $ $(0, y )$ is:
     $$
     \vec{E} = \frac{1}{4 \pi \epsilon_0} \int_{-L}^{L} \frac{\lambda , dx}{r^2} \hat{r}
     $$
   • Here, $ r = \sqrt{x^2 + y^2} $ and $ \hat{r} $ is the unit vector.

Python Code Implementation

Here’s the $Python$ code to compute and visualize the electric field.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50

import numpy as np
import matplotlib.pyplot as plt

# Constants
epsilon_0 = 8.854e-12  # Vacuum permittivity, F/m
lambda_ = 1e-6  # Linear charge density, C/m
L = 1.0  # Length of the line charge, m
k = 1 / (4 * np.pi * epsilon_0)  # Coulomb constant, N·m²/C²

# Field computation grid
x_points = np.linspace(-2, 2, 100)  # x-coordinates
y_points = np.linspace(-2, 2, 100)  # y-coordinates
X, Y = np.meshgrid(x_points, y_points)

# Initialize electric field components
Ex = np.zeros_like(X)
Ey = np.zeros_like(Y)

# Numerical integration for the electric field
num_segments = 100  # Number of segments for the line charge
dx = 2 * L / num_segments  # Segment length
x_line = np.linspace(-L, L, num_segments)  # Line charge x-coordinates

for x_l in x_line:
    r_x = X - x_l
    r_y = Y
    r = np.sqrt(r_x**2 + r_y**2)
    r_hat_x = r_x / r
    r_hat_y = r_y / r
    dE = k * lambda_ * dx / r**2  # Field contribution from each segment
    Ex += dE * r_hat_x
    Ey += dE * r_hat_y

# Normalize vectors for plotting
magnitude = np.sqrt(Ex**2 + Ey**2)
Ex_normalized = Ex / magnitude
Ey_normalized = Ey / magnitude

# Plotting the electric field
plt.figure(figsize=(10, 8))
plt.quiver(X, Y, Ex_normalized, Ey_normalized, magnitude, cmap='inferno', scale=50)
plt.plot([-L, L], [0, 0], color='blue', linewidth=4, label="Line Charge")
plt.colorbar(label="Electric Field Magnitude (N/C)")
plt.title("Electric Field of a Uniform Line Charge")
plt.xlabel("X (m)")
plt.ylabel("Y (m)")
plt.legend()
plt.grid()
plt.axis('equal')
plt.show()

Code Explanation

1. Grid and Initialization:

   • A $2D$ grid of points is created where the electric field will be calculated.
   • The field components $( E_x )$ and $( E_y )$ are initialized to zero.

2. Numerical Integration:

   • The line charge is divided into small segments $( dx )$.
   • For each segment, the contribution to the field at every grid point is calculated using Coulomb’s law:
     $$
     d\vec{E} = \frac{\lambda , dx}{4 \pi \epsilon_0 r^2} \hat{r}
     $$
   • The total electric field is obtained by summing the contributions.

3. Vector Normalization:

   • The field vectors are normalized for plotting, and their magnitude is used as the color map.

4. Plot:

   • The quiver plot shows the direction and relative magnitude of the electric field at various points.
   • The line charge is displayed in blue along the $ x $-axis.

Results Visualization

1. Electric Field Pattern:

   • The electric field vectors point away from the line charge for a positive charge density $( \lambda > 0 )$.
   • The field strength decreases as the distance from the line charge increases.

2. Color Mapping:

   • The color represents the magnitude of the electric field, with warmer colors indicating stronger fields.

This example demonstrates how numerical integration can be used to compute and visualize complex electric fields.

Example Problem

We will simulate the orbit of a planet around a star using Kepler’s laws of planetary motion and Newtonian mechanics.

The goal is to calculate the orbit of the planet and visualize it.

Problem Description

1. Kepler’s Laws:

   • Planets move in elliptical orbits with the star at one focus.
   • The line connecting a planet to the star sweeps out equal areas in equal times.

2. Newtonian Mechanics:

   • The gravitational force between the star and the planet dictates the motion:
     $$
     F = \frac{G M m}{r^2}
     $$
   • The planet’s position is updated using numerical integration of this force.

Python Code Implementation

Here’s the $Python$ implementation:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49

import numpy as np
import matplotlib.pyplot as plt

# Constants
G = 6.67430e-11  # Gravitational constant, m^3 kg^-1 s^-2
M = 1.989e30  # Mass of the star (e.g., the Sun), kg
m = 5.972e24  # Mass of the planet (e.g., the Earth), kg
AU = 1.496e11  # Astronomical Unit, m (average Earth-Sun distance)

# Initial conditions
r0 = AU  # Initial distance, m
v0 = 30000  # Initial tangential velocity, m/s
dt = 60 * 60  # Time step, seconds
time_steps = 10000  # Number of time steps

# Arrays to store position and velocity
pos = np.zeros((time_steps, 2))
vel = np.zeros((time_steps, 2))
acc = np.zeros((time_steps, 2))

# Initial position and velocity
pos[0] = [r0, 0]
vel[0] = [0, v0]

# Simulation loop
for t in range(1, time_steps):
    # Distance and gravitational force
    r = np.linalg.norm(pos[t - 1])
    F = -G * M * m / r**2
    acc[t - 1] = F / m * pos[t - 1] / r

    # Update velocity and position
    vel[t] = vel[t - 1] + acc[t - 1] * dt
    pos[t] = pos[t - 1] + vel[t] * dt

# Convert to AU for visualization
pos_au = pos / AU

# Plotting the orbit
plt.figure(figsize=(8, 8))
plt.plot(pos_au[:, 0], pos_au[:, 1], label="Planet Orbit")
plt.scatter(0, 0, color="orange", label="Star (e.g., Sun)", s=100)
plt.title("Planetary Orbit Simulation")
plt.xlabel("X Position (AU)")
plt.ylabel("Y Position (AU)")
plt.axis("equal")
plt.grid()
plt.legend()
plt.show()

Code Explanation

1. Initial Conditions:

   • The planet starts at a distance $( r_0 ) (1 AU)$ with an initial velocity perpendicular to the radius.

2. Numerical Integration:

   • At each time step:
     • Calculate the gravitational acceleration $( a = F / m )$.
     • Update the velocity using $( v = v + a \cdot dt )$.
     • Update the position using $( r = r + v \cdot dt )$.

3. Scaling:

   • The positions are scaled to astronomical units $(AU)$ for clear visualization.

4. Plot:

   • The orbit is plotted, showing the elliptical path of the planet around the star.

Results Visualization

1. Orbit Shape:

   • The planet’s trajectory forms an ellipse with the star at one focus.

2. Key Observations:

   • The planet moves faster when it is closer to the star (perihelion).
   • The planet moves slower when it is farther from the star (aphelion), illustrating Kepler’s second law.

This simulation demonstrates how planetary orbits emerge from gravitational interactions.

Example Problem

We will solve an example of both elastic and inelastic collisions in one dimension.

The simulation will calculate the velocities after the collision for two masses and illustrate their behavior.

Problem Description

1. Elastic Collision:

   • Two objects, $(m_1)$ and $(m_2)$, collide with initial velocities $(v_{1i})$ and $(v_{2i})$.
   • Conservation of momentum and kinetic energy holds:
     $$
     m_1 v_{1i} + m_2 v_{2i} = m_1 v_{1f} + m_2 v_{2f}
     $$
     $$
     \frac{1}{2} m_1 v_{1i}^2 + \frac{1}{2} m_2 v_{2i}^2 = \frac{1}{2} m_1 v_{1f}^2 + \frac{1}{2} m_2 v_{2f}^2
     $$

2. Inelastic Collision:

   • Two objects stick together after the collision.
   • Only momentum is conserved:
     $$
     m_1 v_{1i} + m_2 v_{2i} = (m_1 + m_2) v_f
     $$

Python Code Implementation

Here’s the code to calculate and visualize the results:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52

import numpy as np
import matplotlib.pyplot as plt

# Define the masses and initial velocities
m1 = 2.0  # mass of object 1
m2 = 1.0  # mass of object 2
v1i = 4.0  # initial velocity of object 1
v2i = -2.0  # initial velocity of object 2

# Elastic collision calculations
v1f_elastic = ((m1 - m2) * v1i + 2 * m2 * v2i) / (m1 + m2)
v2f_elastic = ((m2 - m1) * v2i + 2 * m1 * v1i) / (m1 + m2)

# Inelastic collision calculations
vf_inelastic = (m1 * v1i + m2 * v2i) / (m1 + m2)

# Time and positions for visualization
time = np.linspace(0, 2, 100)
pos1_elastic = v1i * time
pos2_elastic = v2i * time
pos1_elastic_after = v1f_elastic * time
pos2_elastic_after = v2f_elastic * time
pos_inelastic_after = vf_inelastic * time

# Plotting results
plt.figure(figsize=(12, 6))

# Elastic collision
plt.subplot(1, 2, 1)
plt.plot(time, pos1_elastic, label="Object 1 (Before)", linestyle="--")
plt.plot(time, pos2_elastic, label="Object 2 (Before)", linestyle="--")
plt.plot(time, pos1_elastic_after, label="Object 1 (After)")
plt.plot(time, pos2_elastic_after, label="Object 2 (After)")
plt.title("Elastic Collision")
plt.xlabel("Time (s)")
plt.ylabel("Position (m)")
plt.legend()
plt.grid()

# Inelastic collision
plt.subplot(1, 2, 2)
plt.plot(time, pos1_elastic, label="Object 1 (Before)", linestyle="--")
plt.plot(time, pos2_elastic, label="Object 2 (Before)", linestyle="--")
plt.plot(time, pos_inelastic_after, label="Combined (After)")
plt.title("Inelastic Collision")
plt.xlabel("Time (s)")
plt.ylabel("Position (m)")
plt.legend()
plt.grid()

plt.tight_layout()
plt.show()

Code Explanation

1. Elastic Collision:

   • The final velocities $(v_{1f})$ and $(v_{2f})$ are derived using the conservation equations.
   • These velocities are used to calculate the positions over time for visualization.

2. Inelastic Collision:

   • A single final velocity $(v_f)$ is computed since the objects stick together.
   • This velocity is used to show the combined motion.

3. Visualization:

   • Two subplots illustrate the motion before and after both types of collisions.

Results Visualization

1. Elastic Collision:

   • The two objects bounce off each other, maintaining kinetic energy.
   • Their positions diverge after the collision.

2. Inelastic Collision:

   • The two objects stick together and move with the same final velocity.
   • Their positions are represented by a single line post-collision.

This simulation provides a clear comparison of elastic and inelastic collisions.

Problem Description

Simulate the trajectory of a projectile launched at an initial speed $( v_0 )$ and angle $( \theta )$ from the horizontal.

Assume:

1. The motion occurs in a vacuum (no air resistance).
2. The acceleration due to gravity is $( g = 9.8 , \text{m/s}^2 )$.

The equations of motion are:
$$
x(t) = v_0 \cos(\theta) t
$$
$$
y(t) = v_0 \sin(\theta) t - \frac{1}{2} g t^2
$$

• $ x(t) $: Horizontal position.
• $ y(t) $: Vertical position.
• $ t $: Time.

Python Implementation

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28

import numpy as np
import matplotlib.pyplot as plt

# Constants
g = 9.8  # Acceleration due to gravity (m/s^2)
v0 = 50  # Initial speed (m/s)
theta_deg = 45  # Launch angle in degrees
theta_rad = np.radians(theta_deg)  # Convert angle to radians

# Time of flight
t_flight = 2 * v0 * np.sin(theta_rad) / g

# Time array
t = np.linspace(0, t_flight, 500)

# Equations of motion
x = v0 * np.cos(theta_rad) * t
y = v0 * np.sin(theta_rad) * t - 0.5 * g * t**2

# Plot the trajectory
plt.figure(figsize=(10, 6))
plt.plot(x, y, label=f"Launch Angle: {theta_deg}°", color="blue")
plt.title("Projectile Motion Trajectory")
plt.xlabel("Horizontal Distance (m)")
plt.ylabel("Vertical Distance (m)")
plt.grid(True)
plt.legend()
plt.show()

Code Explanation

1. Constants:
   Define gravity $( g )$, initial speed $( v_0 )$, and launch angle $( \theta )$.
2. Time of Flight:
   Calculate the total time the projectile remains in the air using $( t_{\text{flight}} = \frac{2 v_0 \sin(\theta)}{g} )$.
3. Equations of Motion:
   Compute horizontal $( x )$ and vertical $( y )$ positions for a range of time values $( t )$.
4. Plotting:
   Plot the $( x )$ - $( y )$ trajectory to visualize the motion.

Results and Visualization

The graph shows the parabolic trajectory of the projectile:

1. The projectile starts at the origin $( x = 0, y = 0 )$.
2. It reaches a maximum height midway through its flight.
3. The projectile lands back at ground level $( y = 0 )$ after covering a horizontal distance determined by the initial velocity and angle.

This visualization helps understand the fundamental aspects of projectile motion, including its symmetric nature and dependency on the launch angle and speed.

Problem Description

According to Einstein’s theory of special relativity, time dilation occurs for objects moving at significant fractions of the speed of light.

The formula for time dilation is:

$$
t’ = \frac{t}{\sqrt{1 - \frac{v^2}{c^2}}}
$$

• $( t’ )$: Time experienced by the moving observer (proper time).
• $( t )$: Time experienced by the stationary observer.
• $( v )$: Velocity of the moving object.
• $( c )$: Speed of light $(( \approx 3 \times 10^8 , \text{m/s} ))$.

Goal

Simulate how time dilation depends on the velocity $( v )$ of the moving object as a fraction of the speed of light $( c )$.

Python Implementation

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23

import numpy as np
import matplotlib.pyplot as plt

# Constants
speed_of_light = 3e8  # Speed of light (m/s)
stationary_time = 1.0  # Time experienced by the stationary observer (s)

# Velocity fractions of the speed of light
velocity_fractions = np.linspace(0, 0.99, 500)  # From 0 to 99% of c
velocities = velocity_fractions * speed_of_light

# Calculate time dilation
time_dilation = stationary_time / np.sqrt(1 - (velocities**2 / speed_of_light**2))

# Plotting
plt.figure(figsize=(10, 6))
plt.plot(velocity_fractions, time_dilation, color='red', label="Time Dilation")
plt.title("Time Dilation vs. Velocity (Special Relativity)")
plt.xlabel("Velocity as a Fraction of the Speed of Light (v/c)")
plt.ylabel("Dilated Time (t')")
plt.grid(True)
plt.legend()
plt.show()

Code Explanation

1. Constants:
   The speed of light and stationary observer’s time are defined.
2. Velocity Fractions:
   The object’s velocity is expressed as fractions of $( c )$, ranging from $0$ to $0.99$ ($99$% of $( c )$).
3. Time Dilation Formula:
   The time dilation is computed using the relativistic formula.
4. Plotting:
   The relationship between velocity (as a fraction of $( c )$) and time dilation is plotted.

Results and Visualization

The graph illustrates the effect of velocity on time dilation:

1. At low velocities $(( v/c \approx 0 ))$, time dilation is minimal $(( t’ \approx t ))$.
2. As velocity approaches the speed of light $(( v/c \to 1 ))$, time dilation increases drastically $(( t’ \to \infty ))$.
3. This demonstrates how relativistic effects become significant only at very high velocities, a key feature of special relativity.

The double-slit experiment is a classic demonstration of the wave-particle duality of light and matter.

Below, I’ll present an example simulating the interference pattern observed in the experiment, solved using $Python$.

Problem Description

Simulate the intensity distribution of light on a screen in the double-slit experiment.

Assume:

1. The wavelength of light $( \lambda )$ is $500$ nm.
2. The distance between the slits $( d )$ is $0.1$ mm.
3. The screen is $1$ m away from the slits.
4. The slits emit coherent light.

Python Implementation

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34

import numpy as np
import matplotlib.pyplot as plt

# Constants
wavelength = 500e-9  # Wavelength of light (m)
slit_distance = 0.1e-3  # Distance between slits (m)
screen_distance = 1.0  # Distance to the screen (m)
screen_width = 0.02  # Width of the screen (m)
num_points = 1000  # Number of points on the screen

# Create a range of screen positions
x = np.linspace(-screen_width / 2, screen_width / 2, num_points)

# Calculate the path difference
path_difference = (slit_distance * x) / screen_distance

# Calculate the phase difference
phase_difference = (2 * np.pi / wavelength) * path_difference

# Calculate the intensity
intensity = (np.cos(phase_difference / 2)) ** 2

# Normalize intensity
intensity /= np.max(intensity)

# Plot the results
plt.figure(figsize=(10, 6))
plt.plot(x * 1e3, intensity, color='blue', label='Interference Pattern')
plt.title("Double-Slit Experiment Simulation")
plt.xlabel("Position on Screen (mm)")
plt.ylabel("Normalized Intensity")
plt.grid(True)
plt.legend()
plt.show()

Code Explanation

1. Constants: Define physical constants such as the wavelength of light, the distance between the slits, and the screen distance.
2. Screen Positions: Create an array of screen positions (x) to simulate the continuous nature of the screen.
3. Path Difference: Calculate the path difference for light reaching each point on the screen.
4. Phase Difference: Compute the corresponding phase difference.
5. Intensity Calculation: Use the formula $( I = \cos^2(\Delta \phi / 2) )$, normalized to a maximum of $1$.
6. Plotting: Plot the intensity distribution across the screen.

Results and Visualization

The graph shows the intensity distribution as bright and dark fringes.

Bright fringes occur where the waves constructively interfere, while dark fringes correspond to destructive interference.
The central maximum is the brightest fringe, gradually diminishing for higher-order fringes.

The quantum harmonic oscillator is a foundational problem in quantum mechanics.

The time-independent Schrödinger equation for this system is:

$$
-\frac{\hbar^2}{2m} \frac{d^2\psi(x)}{dx^2} + \frac{1}{2} m \omega^2 x^2 \psi(x) = E \psi(x)
$$

• $( \hbar )$: Reduced Planck’s constant
• $( m )$: Mass of the particle
• $( \omega )$: Angular frequency
• $( \psi(x) )$: Wavefunction
• $( E )$: Energy eigenvalue

We’ll solve this equation numerically using the finite difference method and plot the wavefunctions $ \psi(x) $ for the first few energy levels.

Python Code

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42

import numpy as np
import matplotlib.pyplot as plt
from scipy.linalg import eigh

# Constants
hbar = 1.0545718e-34  # Reduced Planck's constant (J·s)
m = 9.10938356e-31  # Mass of electron (kg)
omega = 1e15  # Angular frequency (rad/s)

# Grid setup
x_max = 1e-9  # Maximum x value (meters)
n_points = 1000  # Number of grid points
x = np.linspace(-x_max, x_max, n_points)
dx = x[1] - x[0]  # Spacing

# Kinetic energy operator (finite difference approximation)
kinetic = -(hbar**2 / (2 * m * dx**2)) * (
    np.diag(np.ones(n_points - 1), -1) - 2 * np.diag(np.ones(n_points), 0) + np.diag(np.ones(n_points - 1), 1)
)

# Potential energy operator
potential = 0.5 * m * omega**2 * np.diag(x**2)

# Hamiltonian
hamiltonian = kinetic + potential

# Solving for eigenvalues and eigenvectors
eigenvalues, eigenvectors = eigh(hamiltonian)

# Convert eigenvalues to energy levels (in eV)
eigenvalues_eV = eigenvalues / 1.60218e-19

# Plotting the wavefunctions
plt.figure(figsize=(10, 6))
for n in range(5):  # Plot first 5 wavefunctions
    plt.plot(x * 1e9, eigenvectors[:, n]**2, label=f"n={n}, E={eigenvalues_eV[n]:.2f} eV")
plt.title("Wavefunctions of Quantum Harmonic Oscillator", fontsize=14)
plt.xlabel("x (nm)", fontsize=12)
plt.ylabel("Probability Density", fontsize=12)
plt.legend(fontsize=10)
plt.grid(alpha=0.4)
plt.show()

Explanation of the Code

1. Grid Setup:

   • A finite grid is used to represent the position $( x )$ over a symmetric range.
   • The number of points $( n_{points} )$ determines the resolution.

2. Hamiltonian:

   • The kinetic energy is approximated using a finite difference method.
   • The potential energy is represented as $( \frac{1}{2} m \omega^2 x^2 )$.

3. Solving Eigenvalue Problem:

   • The eigenvalues correspond to energy levels.
   • The eigenvectors are the wavefunctions $( \psi(x) )$.

4. Wavefunction Plot:

   • The first $5$ energy levels’ wavefunctions are plotted as probability densities $( |\psi(x)|^2 )$.

Let’s execute the code to see the results.

Graph Explanation

1. Wavefunctions:

   • The plotted curves represent the probability densities ($(|\psi(x)|^2)$) of the quantum harmonic oscillator for the first $5$ energy levels ($n = 0$ to $n = 4$).
   • Each wavefunction corresponds to a distinct energy eigenvalue $(E_n)$, labeled in electron volts (eV).

2. Nodes:

   • Higher energy levels $(n > 0)$ have more nodes (points where $(|\psi(x)|^2 = 0)$) in their wavefunctions, a characteristic feature of quantum systems.

3. Bound States:

   • The wavefunctions are localized, indicating bound states of the particle within the potential well.

This result aligns with theoretical expectations and provides a clear visualization of quantum mechanical states in a harmonic potential.