Phase Transitions

January 9, 2025

Here’s an example of a Phase Transition in physics, demonstrated using the Ising Model in $Python$.

The Ising Model is a mathematical model used to describe ferromagnetism in statistical mechanics.

Problem

We want to simulate the 2D Ising Model to understand phase transitions.
Specifically, we will:

Simulate a 2D lattice of spins that can either be $+1$ or $-1$.
Study how the magnetization of the system changes with temperature.
Observe the critical temperature at which a phase transition occurs.

Python Code Implementation

import numpy as np
import matplotlib.pyplot as plt

def initialize_lattice(size):
    """Initialize a lattice with random spins (+1 or -1)."""
    return np.random.choice([-1, 1], size=(size, size))

def calculate_energy(lattice):
    """Calculate the total energy of the lattice."""
    energy = 0
    size = lattice.shape[0]
    for i in range(size):
        for j in range(size):
            spin = lattice[i, j]
            neighbors = lattice[(i + 1) % size, j] + lattice[i, (j + 1) % size] + \
                        lattice[(i - 1) % size, j] + lattice[i, (j - 1) % size]
            energy -= spin * neighbors
    return energy / 2  # Each pair is counted twice

def metropolis_step(lattice, beta):
    """Perform a single Metropolis step."""
    size = lattice.shape[0]
    for _ in range(size ** 2):
        i, j = np.random.randint(0, size, size=2)
        spin = lattice[i, j]
        neighbors = lattice[(i + 1) % size, j] + lattice[i, (j + 1) % size] + \
                    lattice[(i - 1) % size, j] + lattice[i, (j - 1) % size]
        delta_energy = 2 * spin * neighbors
        if delta_energy < 0 or np.random.rand() < np.exp(-beta * delta_energy):
            lattice[i, j] *= -1
    return lattice

def calculate_magnetization(lattice):
    """Calculate the magnetization of the lattice."""
    return np.abs(np.sum(lattice)) / lattice.size

def simulate_ising(size, temperatures, steps):
    """Simulate the Ising model for a range of temperatures."""
    lattice = initialize_lattice(size)
    magnetizations = []
    
    for temp in temperatures:
        beta = 1 / temp
        for _ in range(steps):  # Equilibrate
            lattice = metropolis_step(lattice, beta)
        magnetization = calculate_magnetization(lattice)
        magnetizations.append(magnetization)
    
    return magnetizations

# Parameters
size = 20
temperatures = np.linspace(1.5, 3.5, 20)
steps = 1000

# Simulation
magnetizations = simulate_ising(size, temperatures, steps)

# Plotting Results
plt.figure(figsize=(8, 6))
plt.plot(temperatures, magnetizations, marker='o', color='b', label="Magnetization")
plt.axvline(x=2.27, color='r', linestyle='--', label="Critical Temperature (Tc)")
plt.title("Phase Transition in 2D Ising Model")
plt.xlabel("Temperature (T)")
plt.ylabel("Magnetization")
plt.legend()
plt.grid()
plt.show()

Explanation

Initialization:
The lattice is a grid of spins, initialized randomly.
Metropolis Algorithm:
For each spin, we calculate the energy change (ΔE) if the spin is flipped.
The flip is accepted based on the Metropolis criterion, which ensures proper statistical weighting.
Magnetization:
The total magnetization of the system is recorded as the absolute sum of spins normalized by the number of spins.
Temperature Sweep:
By sweeping through a range of temperatures, we observe how the magnetization changes, showing a sharp drop near the critical temperature $(( T_c \approx 2.27 )$ for the 2D Ising Model).

Result

The graph shows magnetization as a function of temperature.

Below the critical temperature $(T_c)$, the system exhibits spontaneous magnetization, indicating a ferromagnetic phase.

Above $(T_c)$, the magnetization drops to zero, showing a transition to the paramagnetic phase.

The red vertical line indicates $(T_c)$.

Entropy Changes

January 8, 2025

Below is an example problem related to Entropy Changes, its solution in $Python$, and a detailed explanation.

Problem: Entropy Changes in Mixing Two Substances

Two substances, A and B, with masses $( m_A = 1 \text{kg} )$ and $( m_B = 2 \text{kg} )$, are mixed.

Both substances are water with a specific heat capacity $( c = 4184 \text{J/(kg·K)} )$.

Initial temperatures: $( T_A = 80^\circ \text{C} )$ and $( T_B = 20^\circ \text{C} )$.
Final equilibrium temperature $( T_f )$ is reached by energy conservation.

Calculate the total entropy change $( \Delta S )$ of the system during this process.

Solution Approach

Use energy conservation to calculate the final equilibrium temperature $( T_f )$:
$$
m_A c (T_A - T_f) = m_B c (T_f - T_B)
$$
Calculate the entropy change for each substance:
$$
\Delta S = \int_{T_i}^{T_f} \frac{dQ}{T}
$$
Since $( dQ = mc , dT )$:
$$
\Delta S = mc \ln\left(\frac{T_f}{T_i}\right)
$$
Sum up $( \Delta S_A )$ and $( \Delta S_B )$ for the total entropy change.

Python Implementation

import numpy as np
import matplotlib.pyplot as plt

# Constants
c = 4184  # Specific heat capacity (J/(kg·K))
m_A = 1.0  # Mass of substance A (kg)
m_B = 2.0  # Mass of substance B (kg)
T_A = 80 + 273.15  # Initial temperature of A (K)
T_B = 20 + 273.15  # Initial temperature of B (K)

# Calculate final equilibrium temperature (energy conservation)
T_f = (m_A * c * T_A + m_B * c * T_B) / (m_A * c + m_B * c)

# Entropy change calculations
def entropy_change(m, c, T_i, T_f):
    return m * c * np.log(T_f / T_i)

delta_S_A = entropy_change(m_A, c, T_A, T_f)
delta_S_B = entropy_change(m_B, c, T_B, T_f)
total_entropy_change = delta_S_A + delta_S_B

# Display results
print(f"Final equilibrium temperature (T_f): {T_f - 273.15:.2f} °C")
print(f"Entropy change for substance A (ΔS_A): {delta_S_A:.2f} J/K")
print(f"Entropy change for substance B (ΔS_B): {delta_S_B:.2f} J/K")
print(f"Total entropy change (ΔS_total): {total_entropy_change:.2f} J/K")

# Visualization
temperatures = np.linspace(T_B, T_A, 100)  # Range of temperatures
entropy_A = m_A * c * np.log(T_f / temperatures)
entropy_B = m_B * c * np.log(T_f / temperatures[::-1])

plt.figure(figsize=(8, 6))
plt.plot(temperatures - 273.15, entropy_A, label="Entropy Change (A)")
plt.plot(temperatures - 273.15, entropy_B, label="Entropy Change (B)")
plt.axhline(y=total_entropy_change, color='r', linestyle='--', label="Total ΔS")
plt.title("Entropy Changes During Mixing")
plt.xlabel("Temperature (°C)")
plt.ylabel("Entropy Change (J/K)")
plt.legend()
plt.grid()
plt.show()

Explanation of the Code

Constants and Inputs:
- The specific heat capacity $( c )$ is set for water.
- Initial temperatures are converted to Kelvin for calculations.
Final Temperature Calculation:
- Using energy conservation, $( T_f )$ is calculated.
Entropy Change Function:
- $ \Delta S = mc \ln(T_f / T_i) $ is implemented for both substances.
Results:
- The entropy change for substances A and B and the total entropy change are displayed.
Graph:
- The graph shows how the entropy of A and B changes with temperature, along with the total entropy change.

Results Interpretation

Final equilibrium temperature (T_f): 40.00 °C
Entropy change for substance A (ΔS_A): -502.96 J/K
Entropy change for substance B (ΔS_B): 552.27 J/K
Total entropy change (ΔS_total): 49.31 J/K

Final Temperature: The equilibrium temperature is approximately $( 40^\circ \text{C} )$.
Entropy Changes: Both substances experience positive entropy changes, as the system moves toward equilibrium.
Total Entropy Change: The total entropy change is positive, consistent with the second law of thermodynamics.

Wave Propagation

January 7, 2025

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

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:

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

January 6, 2025

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

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:

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

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:

Conservation of energy in the simple pendulum
Energy dissipation in the damped pendulum
Nonlinear behavior at large angles
Phase space representation of the dynamics

Magnetic Field of a Current Loop

January 5, 2025

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.

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:

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
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
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
Parameters Used:
- Loop radius ($R$) = $0.1$ meters ($10$ cm)
- Current ($I$) = $1$ ampere
- Calculated field is in Tesla units

The visualizations show that:

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

Spring-Mass System

January 4, 2025

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

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:

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$
Three Cases Simulated:
- Simple Harmonic Motion (no damping, no forcing)
- Damped Oscillation (with damping, no forcing)
- Forced Oscillation (with damping and external force)
Mathematical Model:
- Based on Newton’s Second Law: F = ma
- Differential equation: $mẍ + cẋ + kx = F₀sin(ωt)$
- Solved using scipy’s odeint function
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
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

January 3, 2025

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:

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

Equations

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

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

Constants:
- $ L $, $ m $, and $ h $ are defined in $SI$ units.
Energy Calculation:
- The energy levels $ E_n $ are calculated for $( n = 1, 2, 3 )$ using the energy formula.
Wavefunction Calculation:
- Wavefunctions are computed for each $ n $ using the sinusoidal formula.
Visualization:
- The first graph shows the spatial wavefunctions $ \psi_n(x) $.
- The second graph represents the energy levels.

Graph Explanation

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

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

Results

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

Time Dilation

January 2, 2025

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

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

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.
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 )$.
Dilated Time Calculation:
- Using the time dilation formula, the dilated time is calculated for each velocity fraction.
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

At $ v/c = 0.1 $, the dilated time is approximately $ 1.005 \text{s} $.
At $ v/c = 0.9 $, the dilated time is approximately $ 2.294 \text{s} $.
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

January 1, 2025

Example Problem

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

Problem Description

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

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

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.
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.
Vector Normalization:
- The field vectors are normalized for plotting, and their magnitude is used as the color map.
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

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

December 31, 2024

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

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

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

Initial Conditions:
- The planet starts at a distance $( r_0 ) (1 AU)$ with an initial velocity perpendicular to the radius.
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 )$.
Scaling:
- The positions are scaled to astronomical units $(AU)$ for clear visualization.
Plot:
- The orbit is plotted, showing the elliptical path of the planet around the star.

Results Visualization

Orbit Shape:
- The planet’s trajectory forms an ellipse with the star at one focus.
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.