Wave Propagation

I’ll help you create a simulation of wave propagation, including both traveling waves and standing waves.

We’ll simulate and visualize the wave equation.

Wave Equation:
$$
\frac{\partial^2 y}{\partial t^2} = c^2\frac{\partial^2 y}{\partial x^2}
$$

  • $y(x,t)$ is the displacement
  • $c$ is the wave speed
  • $x$ is position
  • $t$ is time
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
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation
from mpl_toolkits.mplot3d import Axes3D

class WaveSimulator:
    def __init__(self, length=10, nx=1000, dt=0.01, c=1):
        """
        Initialize wave simulation parameters
        length: length of medium
        nx: number of spatial points
        dt: time step
        c: wave speed
        """
        self.length = length
        self.nx = nx
        self.dt = dt
        self.c = c
        
        self.x = np.linspace(0, length, nx)
        self.dx = self.x[1] - self.x[0]
        
        # Check stability condition (CFL condition)
        if c * dt / self.dx > 1:
            raise ValueError("Stability condition not met. Decrease dt or increase dx.")
            
    def traveling_wave(self, t, wavelength, amplitude=1):
        """Generate traveling wave"""
        k = 2 * np.pi / wavelength  # wave number
        omega = self.c * k          # angular frequency
        return amplitude * np.sin(k * self.x - omega * t)
    
    def standing_wave(self, t, mode=1, amplitude=1):
        """Generate standing wave"""
        k = mode * np.pi / self.length
        omega = self.c * k
        return amplitude * np.sin(k * self.x) * np.cos(omega * t)
    
    def wave_packet(self, t, x0=0, sigma=1, k0=5):
        """Generate Gaussian wave packet"""
        x = self.x - (x0 + self.c * t)
        return np.exp(-(x**2)/(2*sigma**2)) * np.cos(k0 * x)

# Create simulation instance
sim = WaveSimulator(length=10, nx=1000, dt=0.01, c=1)

# Time points for simulation
t = np.linspace(0, 10, 100)
X, T = np.meshgrid(sim.x, t)

# Generate different types of waves
traveling_wave = np.array([sim.traveling_wave(ti, wavelength=2) for ti in t])
standing_wave = np.array([sim.standing_wave(ti, mode=2) for ti in t])
wave_packet = np.array([sim.wave_packet(ti, x0=0, sigma=1) for ti in t])

# Create visualizations
fig = plt.figure(figsize=(15, 12))

# 3D plot of traveling wave
ax1 = fig.add_subplot(311, projection='3d')
ax1.plot_surface(X, T, traveling_wave, cmap='viridis')
ax1.set_xlabel('Position (x)')
ax1.set_ylabel('Time (t)')
ax1.set_zlabel('Amplitude')
ax1.set_title('Traveling Wave')

# 3D plot of standing wave
ax2 = fig.add_subplot(312, projection='3d')
ax2.plot_surface(X, T, standing_wave, cmap='viridis')
ax2.set_xlabel('Position (x)')
ax2.set_ylabel('Time (t)')
ax2.set_zlabel('Amplitude')
ax2.set_title('Standing Wave (Mode 2)')

# 3D plot of wave packet
ax3 = fig.add_subplot(313, projection='3d')
ax3.plot_surface(X, T, wave_packet, cmap='viridis')
ax3.set_xlabel('Position (x)')
ax3.set_ylabel('Time (t)')
ax3.set_zlabel('Amplitude')
ax3.set_title('Gaussian Wave Packet')

plt.tight_layout()
plt.show()

# Create animation of wave propagation
fig, (ax1, ax2, ax3) = plt.subplots(3, 1, figsize=(10, 12))

# Initialize lines
line1, = ax1.plot([], [], 'b-', label='Traveling Wave')
line2, = ax2.plot([], [], 'r-', label='Standing Wave')
line3, = ax3.plot([], [], 'g-', label='Wave Packet')

# Set axis limits
for ax in [ax1, ax2, ax3]:
    ax.set_xlim(0, sim.length)
    ax.set_ylim(-1.5, 1.5)
    ax.grid(True)
    ax.legend()

ax1.set_title('Traveling Wave')
ax2.set_title('Standing Wave')
ax3.set_title('Wave Packet')

# Animation function
def animate(frame):
    line1.set_data(sim.x, traveling_wave[frame])
    line2.set_data(sim.x, standing_wave[frame])
    line3.set_data(sim.x, wave_packet[frame])
    return line1, line2, line3

# Create animation
anim = FuncAnimation(fig, animate, frames=len(t), interval=50, blit=True)
plt.tight_layout()
plt.show()

Let me explain the key components of this simulation:

  1. Wave Types Implemented:

    a) Traveling Wave:

    • Form: $y(x,t) = A\sin(kx - \omega t)$
    • Moves with constant speed in one direction
    • Maintains shape while propagating

    b) Standing Wave:

    • Form: $y(x,t) = A\sin(kx)\cos(\omega t)$
    • Fixed nodes and antinodes
    • Result of interference between waves traveling in opposite directions

    c) Gaussian Wave Packet:

    • Localized wave group
    • Combines carrier wave with Gaussian envelope
    • Shows dispersion effects

Pendulum Dynamics

I’ll help you create a $Python$ program to simulate and visualize the dynamics of a pendulum, including both simple and damped pendulum cases.

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
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
import numpy as np
from scipy.integrate import odeint
import matplotlib.pyplot as plt

class PendulumSimulator:
    def __init__(self, L=1.0, g=9.81, m=1.0, b=0.0):
        """
        Initialize pendulum parameters
        L: length of pendulum (m)
        g: gravitational acceleration (m/s^2)
        m: mass of bob (kg)
        b: damping coefficient (kg/s)
        """
        self.L = L
        self.g = g
        self.m = m
        self.b = b

    def derivatives(self, state, t):
        """
        Calculate derivatives for the pendulum system
        state[0] = theta (angle)
        state[1] = omega (angular velocity)
        """
        theta, omega = state
        
        # Calculate acceleration (second derivative of theta)
        # Include nonlinear pendulum dynamics and damping
        theta_dot = omega
        omega_dot = (-self.b * omega / self.m 
                    - (self.g / self.L) * np.sin(theta))
        
        return [theta_dot, omega_dot]

    def simulate(self, theta0, omega0, t_span):
        """
        Simulate pendulum motion
        theta0: initial angle (radians)
        omega0: initial angular velocity (rad/s)
        t_span: time points for simulation
        """
        initial_state = [theta0, omega0]
        solution = odeint(self.derivatives, initial_state, t_span)
        return solution

# Create simulation instances
simple_pendulum = PendulumSimulator(L=1.0, b=0.0)  # No damping
damped_pendulum = PendulumSimulator(L=1.0, b=0.5)  # With damping

# Time points for simulation
t = np.linspace(0, 10, 1000)

# Initial conditions
theta0 = np.pi/3  # 60 degrees
omega0 = 0.0      # starting from rest

# Run simulations
simple_solution = simple_pendulum.simulate(theta0, omega0, t)
damped_solution = damped_pendulum.simulate(theta0, omega0, t)

# Create visualization
fig = plt.figure(figsize=(15, 10))

# Plot angular displacement
ax1 = fig.add_subplot(221)
ax1.plot(t, simple_solution[:, 0], 'b-', label='Simple Pendulum')
ax1.plot(t, damped_solution[:, 0], 'r--', label='Damped Pendulum')
ax1.set_xlabel('Time (s)')
ax1.set_ylabel('Angle (radians)')
ax1.set_title('Angular Displacement vs Time')
ax1.grid(True)
ax1.legend()

# Plot angular velocity
ax2 = fig.add_subplot(222)
ax2.plot(t, simple_solution[:, 1], 'b-', label='Simple Pendulum')
ax2.plot(t, damped_solution[:, 1], 'r--', label='Damped Pendulum')
ax2.set_xlabel('Time (s)')
ax2.set_ylabel('Angular Velocity (rad/s)')
ax2.set_title('Angular Velocity vs Time')
ax2.grid(True)
ax2.legend()

# Phase space plot
ax3 = fig.add_subplot(223)
ax3.plot(simple_solution[:, 0], simple_solution[:, 1], 'b-', label='Simple Pendulum')
ax3.plot(damped_solution[:, 0], damped_solution[:, 1], 'r--', label='Damped Pendulum')
ax3.set_xlabel('Angle (radians)')
ax3.set_ylabel('Angular Velocity (rad/s)')
ax3.set_title('Phase Space')
ax3.grid(True)
ax3.legend()

# Energy plot
def calculate_energy(solution, pendulum):
    theta, omega = solution[:, 0], solution[:, 1]
    KE = 0.5 * pendulum.m * (pendulum.L * omega)**2
    PE = pendulum.m * pendulum.g * pendulum.L * (1 - np.cos(theta))
    return KE + PE

ax4 = fig.add_subplot(224)
E_simple = calculate_energy(simple_solution, simple_pendulum)
E_damped = calculate_energy(damped_solution, damped_pendulum)
ax4.plot(t, E_simple, 'b-', label='Simple Pendulum')
ax4.plot(t, E_damped, 'r--', label='Damped Pendulum')
ax4.set_xlabel('Time (s)')
ax4.set_ylabel('Total Energy (J)')
ax4.set_title('Total Energy vs Time')
ax4.grid(True)
ax4.legend()

plt.tight_layout()
plt.show()

Let me explain the code and the physics:

  1. Physical System:

    • The pendulum consists of a mass m at the end of a rigid, massless rod of length $L$
    • Motion is governed by gravity ($g$) and optional damping ($b$)
    • The equation of motion is: $θ̈ + (b/m)θ̇ + (g/L)sin(θ) = 0$

  2. Code Structure:

    • PendulumSimulator class handles the physics calculations
    • Uses scipy.integrate.odeint for numerical integration of the equations of motion
    • Simulates both simple (undamped) and damped pendulum cases

  3. Visualizations (4 plots):

a) Angular Displacement vs Time:
- Shows how the angle changes over time
- Simple pendulum oscillates continuously
- Damped pendulum’s oscillations decay

b) Angular Velocity vs Time:
- Shows how the angular velocity changes
- Simple pendulum maintains amplitude
- Damped pendulum shows decreasing velocity magnitude

c) Phase Space Plot:
- Shows relationship between angle and angular velocity
- Simple pendulum forms closed orbit (conservation of energy)
- Damped pendulum spirals toward center (energy dissipation)

d) Total Energy vs Time:
- Shows sum of kinetic and potential energy
- Simple pendulum maintains constant energy
- Damped pendulum shows energy decay

  1. Parameters Used:
    • Length ($L$) = $1.0$ m
    • Initial angle = $π/3$ radians ($60$ degrees)
    • Initial angular velocity = $0$
    • Damping coefficient (b) = $0.5$ kg/s for damped case

The simulation shows key physical principles:

  1. Conservation of energy in the simple pendulum
  2. Energy dissipation in the damped pendulum
  3. Nonlinear behavior at large angles
  4. Phase space representation of the dynamics

Magnetic Field of a Current Loop

I’ll help you create a $Python$ program to calculate and visualize the magnetic field of a current loop.

We’ll solve a specific example and create visualizations to better understand the field distribution.

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
83
84
85
86
87
88
89
90
91
92
93
94
95
96
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

def magnetic_field_loop(x, y, z, R, I):
    """
    Calculate the magnetic field components (Bx, By, Bz) at point (x, y, z)
    for a circular current loop of radius R carrying current I
    centered at the origin in the x-y plane.
    
    Parameters:
    x, y, z: coordinates of the point where B-field is calculated (meters)
    R: radius of the current loop (meters)
    I: current in the loop (amperes)
    
    Returns:
    tuple: (Bx, By, Bz) magnetic field components in Tesla
    """
    mu0 = 4 * np.pi * 1e-7  # permeability of free space
    
    # Calculate intermediate values
    r = np.sqrt(x**2 + y**2 + z**2)
    if r == 0:
        return (0, 0, 0)
    
    # Calculate components
    factor = mu0 * I * R**2 / (4 * np.pi * (R**2 + r**2)**(3/2))
    
    Bx = factor * 3 * x * z / r**2
    By = factor * 3 * y * z / r**2
    Bz = factor * (3 * z**2 / r**2 - 1)
    
    return (Bx, By, Bz)

# Example parameters
R = 0.1  # radius of loop = 10 cm
I = 1.0  # current = 1 A

# Create grid of points
x = np.linspace(-0.2, 0.2, 20)
y = np.linspace(-0.2, 0.2, 20)
z = np.linspace(-0.2, 0.2, 20)
X, Y, Z = np.meshgrid(x, y, z)

# Calculate B-field components
Bx = np.zeros_like(X)
By = np.zeros_like(Y)
Bz = np.zeros_like(Z)

for i in range(len(x)):
    for j in range(len(y)):
        for k in range(len(z)):
            Bx[i,j,k], By[i,j,k], Bz[i,j,k] = magnetic_field_loop(X[i,j,k], Y[i,j,k], Z[i,j,k], R, I)

# Create plots
fig = plt.figure(figsize=(15, 5))

# XZ-plane plot (y = 0)
ax1 = fig.add_subplot(131)
y_mid = len(y)//2
magnitude = np.sqrt(Bx[:,y_mid,:]**2 + Bz[:,y_mid,:]**2)
ax1.streamplot(x, z, Bx[:,y_mid,:], Bz[:,y_mid,:], 
              color=magnitude, cmap='viridis')
ax1.set_title('Magnetic field lines in XZ-plane (y = 0)')
ax1.set_xlabel('x (m)')
ax1.set_ylabel('z (m)')
ax1.add_patch(plt.Circle((0,0), R, fill=False, color='red'))

# Field strength along z-axis
ax2 = fig.add_subplot(132)
z_axis = np.linspace(-0.2, 0.2, 100)
Bz_axis = [magnetic_field_loop(0, 0, z_val, R, I)[2] for z_val in z_axis]
ax2.plot(z_axis, Bz_axis)
ax2.set_title('Bz along z-axis')
ax2.set_xlabel('z (m)')
ax2.set_ylabel('Bz (T)')
ax2.grid(True)

# 3D vector plot
ax3 = fig.add_subplot(133, projection='3d')
skip = 3
ax3.quiver(X[::skip,::skip,::skip], 
           Y[::skip,::skip,::skip], 
           Z[::skip,::skip,::skip],
           Bx[::skip,::skip,::skip], 
           By[::skip,::skip,::skip], 
           Bz[::skip,::skip,::skip],
           length=0.05, normalize=True)
ax3.set_title('3D magnetic field vectors')
ax3.set_xlabel('x (m)')
ax3.set_ylabel('y (m)')
ax3.set_zlabel('z (m)')

plt.tight_layout()
plt.show()

Let me explain the code and the physics behind it:

  1. Physics Background:

    • A circular current loop generates a magnetic field described by the $Biot$-$Savart$ $law$
    • The field has axial symmetry around the $z$-axis
    • The field is strongest near the loop and decreases with distance
    • Inside the loop, the field is mostly uniform and parallel to the $z$-axis

  2. Code Structure:

    • The magnetic_field_loop function calculates the magnetic field components (Bx, By, Bz) at any point in space
    • We use the analytical solution for the magnetic field of a current loop
    • The calculation includes the permeability of free space (μ₀)
    • We create a $3D$ grid of points to evaluate the field

  3. Visualizations:

    • Left plot: Magnetic field lines in the $xz$-plane ($y = 0$)

      • The red circle represents the current loop
      • The streamlines show the direction of the magnetic field
      • Colors indicate field strength

    • Middle plot: Bz component along the $z$-axis

      • Shows how the field strength varies along the axis of the loop
      • The field is strongest near the loop and decreases with distance

    • Right plot: $3D$ vector representation

      • Arrows show the direction of the magnetic field in $3D$ space
      • The length of the arrows is normalized for better visualization

  4. Parameters Used:

    • Loop radius ($R$) = $0.1$ meters ($10$ cm)
    • Current ($I$) = $1$ ampere
    • Calculated field is in Tesla units

The visualizations show that:

  1. The field is symmetric around the $z$-axis
  2. Field lines emerge from one side of the loop and curve around to enter the other side
  3. The field is strongest near the loop and decreases with distance
  4. Inside the loop, the field is relatively uniform and points along the $z$-axis

Spring-Mass System

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

Quantum Mechanics - Particle in a Box

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 $.

  1. 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.

Time Dilation

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.

Electric Field Due to a Continuous Line Charge

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.

Planetary Orbits Using Kepler's Laws

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.

Elastic and Inelastic Collisions in One Dimension

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.

Projectile Motion

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.