Quantum Harmonic Oscillator

March 23, 2025

Numerical Simulation and Visualization with Python

I’ll solve a practical quantum mechanics problem using $Python$ in Google Colab, including clear explanations and visualization.

import numpy as np
import matplotlib.pyplot as plt
from scipy import sparse
from scipy.linalg import expm
from scipy.special import hermite, factorial

# Set plotting parameters for better visualization
plt.rcParams['figure.figsize'] = [12, 8]
plt.rcParams['font.size'] = 14

# Physical constants
hbar = 1.0  # Reduced Planck constant (using natural units)
m = 1.0     # Mass of the particle
omega = 1.0 # Angular frequency of the oscillator

# Simulation parameters
n_points = 1000  # Number of position points for discretization
x_range = 10.0   # Position range for simulation
dx = 2.0 * x_range / n_points  # Position step size
x = np.linspace(-x_range, x_range, n_points)  # Position array

# Create potential energy function (harmonic oscillator: V(x) = 0.5 * m * omega^2 * x^2)
V = 0.5 * m * omega**2 * x**2

# Create second derivative operator (kinetic energy part)
diag = np.ones(n_points) * (2.0 / dx**2)
off_diag = np.ones(n_points-1) * (-1.0 / dx**2)
D2 = sparse.diags([off_diag, diag, off_diag], [-1, 0, 1])

# Create Hamiltonian operator: H = -hbar^2/(2m) * d^2/dx^2 + V(x)
H = sparse.diags(V) - (hbar**2 / (2.0 * m)) * D2

# Find eigenvalues and eigenvectors of the Hamiltonian
eigenvalues, eigenvectors = sparse.linalg.eigs(H, k=10, which='SM')  # Find lowest 10 energy states
idx = eigenvalues.argsort()  # Sort by energy
eigenvalues = eigenvalues[idx].real
eigenvectors = eigenvectors[:, idx]

# Normalize eigenvectors
for i in range(len(eigenvalues)):
    eigenvectors[:, i] = eigenvectors[:, i] / np.sqrt(np.sum(np.abs(eigenvectors[:, i])**2) * dx)

# Analytical solutions for comparison
def analytical_wavefunction(n, x):
    """
    Calculate the analytical wavefunction for quantum harmonic oscillator
    for quantum number n at position x
    """
    # Constants for the harmonic oscillator
    alpha = m * omega / hbar
    
    # Normalization constant
    normalization = (alpha / np.pi)**0.25 * (1.0 / np.sqrt(2**n * factorial(n)))
    
    # Hermite polynomial
    herm = hermite(n)(np.sqrt(alpha) * x)
    
    # Full wavefunction
    psi = normalization * herm * np.exp(-alpha * x**2 / 2.0)
    
    return psi

# Calculate analytical energy levels
def analytical_energy(n):
    """Calculate analytical energy for quantum harmonic oscillator with quantum number n"""
    return hbar * omega * (n + 0.5)

# Compute analytical solutions
analytical_energies = [analytical_energy(n) for n in range(10)]
analytical_wavefunctions = [analytical_wavefunction(n, x) for n in range(10)]

# Plot and compare numerical and analytical solutions
fig, axs = plt.subplots(5, 2, figsize=(15, 20))
axs = axs.flatten()

for i in range(10):
    axs[i].plot(x, eigenvectors[:, i].real, 'b-', label=f'Numerical ψ{i}')
    axs[i].plot(x, analytical_wavefunctions[i], 'r--', label=f'Analytical ψ{i}')
    axs[i].axhline(y=0, color='k', linestyle='-', alpha=0.3)
    axs[i].set_title(f'Energy Level {i}: E_numerical = {eigenvalues[i]:.5f}, E_analytical = {analytical_energies[i]:.5f}')
    axs[i].set_xlabel('Position (x)')
    axs[i].set_ylabel('Wave Function (ψ)')
    axs[i].legend()
    axs[i].grid(True, alpha=0.3)

plt.tight_layout()
plt.show()

# Time evolution of a wavepacket
# Create a Gaussian wavepacket as initial state (superposition of ground and first excited states)
psi_0 = 0.7 * eigenvectors[:, 0] + 0.7 * eigenvectors[:, 1]
psi_0 = psi_0 / np.sqrt(np.sum(np.abs(psi_0)**2) * dx)  # Normalize

# Time parameters
t_total = 10.0  # Total simulation time
n_steps = 100   # Number of time steps
dt = t_total / n_steps  # Time step

# Store wavefunction at different times
psi_t = np.zeros((n_steps, n_points), dtype=complex)
psi_t[0] = psi_0

# Time evolution operator
def time_evolution(psi, H, dt):
    """Time evolution using the operator exp(-i*H*dt/hbar)"""
    # We use the fact that eigenvectors of H provide a basis for expansion
    coeffs = np.zeros(10, dtype=complex)
    for i in range(10):
        coeffs[i] = np.sum(np.conj(eigenvectors[:, i]) * psi) * dx
    
    # Evolve each component and recombine
    psi_evolved = np.zeros_like(psi, dtype=complex)
    for i in range(10):
        psi_evolved += coeffs[i] * np.exp(-1j * eigenvalues[i] * dt / hbar) * eigenvectors[:, i]
    
    return psi_evolved

# Perform time evolution
for t in range(1, n_steps):
    psi_t[t] = time_evolution(psi_t[t-1], H, dt)

# Probability densities for animation
prob_densities = np.abs(psi_t)**2

# Plot time evolution
fig, ax = plt.subplots(figsize=(12, 8))

# Plot potential energy curve (scaled for visibility)
V_scaled = V / np.max(V) * np.max(prob_densities) * 0.8
ax.plot(x, V_scaled, 'k-', alpha=0.3, label='Potential (scaled)')

# Select specific timesteps to plot
time_points = [0, 25, 50, 75, 99]
colors = ['b', 'g', 'r', 'c', 'm']

for i, t in enumerate(time_points):
    ax.plot(x, prob_densities[t], colors[i], label=f't = {t*dt:.2f}')

ax.set_xlabel('Position (x)')
ax.set_ylabel('Probability Density |ψ(x,t)|²')
ax.set_title('Time Evolution of Quantum Harmonic Oscillator Wavepacket')
ax.legend()
ax.grid(True, alpha=0.3)

plt.tight_layout()
plt.show()

# Create an energy level diagram
n_levels = 10
fig, ax = plt.subplots(figsize=(8, 10))

# Plot energy levels
for n in range(n_levels):
    E = analytical_energy(n)
    ax.plot([-0.5, 0.5], [E, E], 'b-', linewidth=2)
    ax.text(0.6, E, f'n = {n}, E = {E:.2f}ħω', va='center')

# Add potential well visualization
x_pot = np.linspace(-0.4, 0.4, 100)
V_pot = 0.5 * m * omega**2 * x_pot**2
ax.plot(x_pot, V_pot, 'r--', label='V(x) = ½mω²x²')

ax.set_xlim(-1, 1.5)
ax.set_ylim(-0.5, analytical_energy(n_levels-1) + 1)
ax.set_ylabel('Energy (in units of ħω)')
ax.set_title('Energy Levels of Quantum Harmonic Oscillator')
ax.set_xticks([])
ax.legend(loc='upper right')
ax.grid(True, alpha=0.3)

plt.tight_layout()
plt.show()

Quantum Harmonic Oscillator Simulation

In this example, I’ll solve the quantum harmonic oscillator problem, which is one of the most important systems in quantum mechanics.

It models a particle moving in a parabolic potential well, like a mass on a spring at the quantum scale.

Theoretical Background

The quantum harmonic oscillator is described by the Hamiltonian:

$$\hat{H} = -\frac{\hbar^2}{2m}\frac{d^2}{dx^2} + \frac{1}{2}m\omega^2 x^2$$

Where:

$\hbar$ is the reduced Planck constant
$m$ is the mass of the particle
$\omega$ is the angular frequency of the oscillator
$x$ is the position

The energy eigenvalues are given by:

$$E_n = \hbar\omega(n + \frac{1}{2})$$

And the wavefunctions are:

$$\psi_n(x) = \left(\frac{\alpha}{\pi}\right)^{1/4}\frac{1}{\sqrt{2^n n!}}H_n(\sqrt{\alpha}x)e^{-\alpha x^2/2}$$

Where:

$\alpha = \frac{m\omega}{\hbar}$
$H_n$ are the Hermite polynomials
$n$ is the quantum number $(0, 1, 2, …)$

Code Explanation

The code solves this problem using two methods:

Numerical solution: Using a discretized position grid and sparse matrix techniques
Analytical solution: Calculating the exact expressions for comparison

Key Components:

Set up the position grid and potential:
- Create a discretized position space
- Define the harmonic potential $V(x) = \frac{1}{2}m\omega^2 x^2$
Create the Hamiltonian matrix:
- The kinetic energy term uses a finite difference approximation for the second derivative
- The potential energy is added as a diagonal matrix
Find eigenvalues and eigenvectors:
- Compute the lowest $10$ energy states
- Sort them by energy
Compare with analytical solutions:
- Calculate the exact energy levels $E_n = \hbar\omega(n + \frac{1}{2})$
- Calculate the exact wavefunctions using Hermite polynomials
Simulate time evolution:
- Create a wavepacket as a superposition of ground and first excited states
- Evolve it through time and visualize the probability density

Results and Visualization

The code produces three main visualizations:

Wavefunction Comparison:
The first plot shows the first $10$ energy eigenstates, comparing the numerical solutions (blue solid lines) with the analytical solutions (red dashed lines).
The close match verifies our numerical method is working correctly.

Time Evolution:
The second plot shows the time evolution of a wavepacket (initially a superposition of ground and first excited states).
The probability density is shown at five different times, demonstrating the oscillatory behavior of the wavepacket.

Energy Level Diagram:
The final plot shows the quantized energy levels of the harmonic oscillator.
Notice how they are equally spaced by $\hbar\omega$, which is a key characteristic of the quantum harmonic oscillator.

Physical Interpretation

Quantization of Energy: Unlike a classical oscillator that can have any energy, the quantum harmonic oscillator can only have discrete energy values.
Zero-Point Energy: Even in the ground state (n=0), the energy is $E_0 = \frac{1}{2}\hbar\omega$, not zero. This represents the quantum mechanical uncertainty principle at work.
Probability Distribution: The wavefunction squared gives the probability density of finding the particle at a particular position.
For the ground state, this is concentrated around the center, but for higher energy states, there are nodes where the probability is zero.
Time Evolution: The time-dependent wavepacket oscillates back and forth, similar to a classical oscillator, but maintains its quantum nature with the probability spread according to the uncertainty principle.

This simulation provides a clear visualization of these fundamental quantum mechanical concepts using a system that has direct analogies to classical mechanics.

Geometric Brownian Motion

March 22, 2025

Stock Price Simulation and Probability Analysis

I’ll solve a practical probability theory problem using $Python$, present the mathematical formulas, and visualize the results with graphs.

import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
from scipy import stats

# Set random seed for reproducibility
np.random.seed(42)

# Parameters for our stock price model
initial_price = 100
drift = 0.0002  # Daily drift (expected return)
volatility = 0.015  # Daily volatility
trading_days = 252  # Number of trading days in a year
num_simulations = 1000  # Number of simulations to run
time_horizon = 252  # Simulate for 1 year

# Create an array to store all simulation results
all_simulations = np.zeros((num_simulations, time_horizon + 1))
all_simulations[:, 0] = initial_price  # Set initial price

# Perform the simulations using Geometric Brownian Motion
for i in range(num_simulations):
    # Generate random shocks
    random_shocks = np.random.normal(0, 1, time_horizon)
    
    # Calculate daily returns
    daily_returns = np.exp(drift - 0.5 * volatility**2 + volatility * random_shocks)
    
    # Calculate price path
    price_path = initial_price * np.cumprod(daily_returns)
    
    # Store the simulation
    all_simulations[i, 1:] = price_path

# Calculate statistics
final_prices = all_simulations[:, -1]
mean_final_price = np.mean(final_prices)
median_final_price = np.median(final_prices)
std_final_price = np.std(final_prices)

# Calculate probability of ending above certain thresholds
prob_above_initial = np.mean(final_prices > initial_price)
prob_above_110 = np.mean(final_prices > 110)
prob_above_120 = np.mean(final_prices > 120)
prob_below_90 = np.mean(final_prices < 90)

# Calculate confidence intervals
confidence_level = 0.95
lower_percentile = (1 - confidence_level) / 2
upper_percentile = 1 - lower_percentile
confidence_interval = np.percentile(final_prices, [lower_percentile * 100, upper_percentile * 100])

# Plotting
plt.figure(figsize=(15, 10))

# Plot 1: Sample paths
plt.subplot(2, 2, 1)
for i in range(50):  # Plot just 50 sample paths for clarity
    plt.plot(all_simulations[i], linewidth=0.8, alpha=0.6)
plt.title('Sample Stock Price Paths')
plt.xlabel('Trading Days')
plt.ylabel('Stock Price')
plt.grid(True)

# Plot 2: Histogram of final prices
plt.subplot(2, 2, 2)
sns.histplot(final_prices, kde=True, bins=30)
plt.axvline(initial_price, color='r', linestyle='--', label='Initial Price')
plt.axvline(mean_final_price, color='g', linestyle='-', label='Mean Final Price')
plt.axvline(confidence_interval[0], color='orange', linestyle=':', label=f'{confidence_level*100}% CI Lower')
plt.axvline(confidence_interval[1], color='orange', linestyle=':', label=f'{confidence_level*100}% CI Upper')
plt.title('Distribution of Final Stock Prices')
plt.xlabel('Stock Price')
plt.ylabel('Frequency')
plt.legend()
plt.grid(True)

# Plot 3: Empirical CDF
plt.subplot(2, 2, 3)
x = np.sort(final_prices)
y = np.arange(1, len(x) + 1) / len(x)
plt.plot(x, y)
plt.axvline(initial_price, color='r', linestyle='--', label='Initial Price')
plt.axhline(prob_above_initial, color='b', linestyle='--')
plt.axvline(110, color='g', linestyle='--', label='$110')
plt.axhline(prob_above_110, color='b', linestyle=':')
plt.title('Empirical Cumulative Distribution Function')
plt.xlabel('Stock Price')
plt.ylabel('Probability')
plt.grid(True)
plt.legend()

# Plot 4: Q-Q plot to check for normality of returns
plt.subplot(2, 2, 4)
log_returns = np.log(all_simulations[:, 1:] / all_simulations[:, :-1])
log_returns = log_returns.flatten()
stats.probplot(log_returns, dist="norm", plot=plt)
plt.title('Q-Q Plot of Log Returns')
plt.grid(True)

plt.tight_layout()
plt.savefig('stock_price_simulation.png', dpi=300)
plt.show()

# Display the statistical results
print(f"Initial stock price: ${initial_price:.2f}")
print(f"Mean final price: ${mean_final_price:.2f}")
print(f"Median final price: ${median_final_price:.2f}")
print(f"Standard deviation of final price: ${std_final_price:.2f}")
print(f"95% confidence interval: [${confidence_interval[0]:.2f}, ${confidence_interval[1]:.2f}]")
print(f"Probability of ending above initial price (${initial_price:.2f}): {prob_above_initial:.2%}")
print(f"Probability of ending above $110: {prob_above_110:.2%}")
print(f"Probability of ending above $120: {prob_above_120:.2%}")
print(f"Probability of ending below $90: {prob_below_90:.2%}")

Stock Price Movement Simulation Using Geometric Brownian Motion

In this example, I’ll simulate stock price movements using a probability model called Geometric Brownian Motion (GBM), which is widely used in financial mathematics and forms the foundation of the Black-Scholes option pricing model.

Mathematical Model

The stock price movement follows Geometric Brownian Motion, which can be expressed as:

$$dS_t = \mu S_t dt + \sigma S_t dW_t$$

Where:

$S_t$ is the stock price at time $t$
$\mu$ is the drift (expected return)
$\sigma$ is the volatility
$W_t$ is a Wiener process (standard Brownian motion)

The solution to this stochastic differential equation gives us the stock price at time $t$:

$$S_t = S_0 \exp\left(\left(\mu - \frac{\sigma^2}{2}\right)t + \sigma W_t\right)$$

For discrete time simulations, we can use:

$$S_{t+\Delta t} = S_t \exp\left(\left(\mu - \frac{\sigma^2}{2}\right)\Delta t + \sigma \sqrt{\Delta t} Z_t\right)$$

Where $Z_t \sim N(0,1)$ is a standard normal random variable.

Code Explanation

Setup and Parameters:
- We set initial stock price to $100
- Drift (µ) of $0.0002$ per day (approximately $5$% annual return)
- Volatility (σ) of $0.015$ per day (approximately $24$% annualized)
- Simulating $1,000$ possible price paths over $252$ trading days ($1$ year)
Simulation:
- We generate random shocks from a normal distribution
- Calculate daily returns using the GBM formula
- Compute cumulative price paths starting from the initial price
- Store all simulations in an array
Statistical Analysis:
- Calculate mean, median, and standard deviation of final prices
- Determine probabilities of the price exceeding certain thresholds
- Calculate 95% confidence intervals
Visualization:
- Sample price paths to visualize potential trajectories
- Histogram showing the distribution of final prices
- Empirical cumulative distribution function (ECDF)
- Q-Q plot to check if returns follow a normal distribution

Initial stock price: $100.00
Mean final price: $105.10
Median final price: $103.21
Standard deviation of final price: $24.78
95% confidence interval: [$63.76, $159.43]
Probability of ending above initial price ($100.00): 56.10%
Probability of ending above $110: 37.20%
Probability of ending above $120: 24.60%
Probability of ending below $90: 29.10%

Key Insights from the Visualization

Sample Paths:
The first plot shows $50$ possible price paths, illustrating the random nature of stock movements while maintaining the overall drift.
Final Price Distribution:
The histogram shows that final prices follow a log-normal distribution (which is expected from GBM).
The mean final price is typically higher than the median due to the right-skewed nature of the distribution.
Probabilities:
The ECDF graph shows cumulative probabilities, allowing us to read off values like:
- Probability of ending above $100 (initial price)
- Probability of ending above $110
- Probability of ending below $90
Log Returns Normality:
The Q-Q plot confirms that log returns closely follow a normal distribution, which is a key assumption of the GBM model.

Practical Applications

This simulation helps answer questions like:

What is the probability of a stock gaining more than $10$% in a year?
What is the expected range of stock prices after one year (confidence interval)?
What is the risk of the stock falling below a certain threshold?

These probability calculations are essential for risk management, option pricing, portfolio optimization, and investment decision-making.

The simulation demonstrates how seemingly complex financial phenomena can be modeled using probability theory and simulated using $Python$, providing actionable insights for investors and risk managers.

Option Pricing Analysis

March 21, 2025

Implementing the Black-Scholes Model with Python

I’ll provide you with a financial mathematics example and solve it using $Python$, including proper mathematical notation and visualization.

import numpy as np
import matplotlib.pyplot as plt
import scipy.stats as stats
from matplotlib import cm
import pandas as pd
from mpl_toolkits.mplot3d import Axes3D

# Black-Scholes Option Pricing Model
def black_scholes(S, K, T, r, sigma, option='call'):
    """
    Calculate Black-Scholes option price for a call or put
    
    Parameters:
    S: Current stock price
    K: Strike price
    T: Time to maturity (in years)
    r: Risk-free interest rate
    sigma: Volatility
    option: 'call' or 'put'
    
    Returns:
    Option price
    """
    d1 = (np.log(S / K) + (r + 0.5 * sigma**2) * T) / (sigma * np.sqrt(T))
    d2 = d1 - sigma * np.sqrt(T)
    
    if option == 'call':
        price = S * stats.norm.cdf(d1) - K * np.exp(-r * T) * stats.norm.cdf(d2)
    else:
        price = K * np.exp(-r * T) * stats.norm.cdf(-d2) - S * stats.norm.cdf(-d1)
        
    return price

# Example 1: Visualize call option price as a function of stock price and time to maturity
# Parameters
S_range = np.linspace(80, 120, 50)  # Stock price range
T_range = np.linspace(0.1, 1, 50)   # Time to maturity range
K = 100                             # Strike price
r = 0.05                            # Risk-free rate
sigma = 0.2                         # Volatility

# Create meshgrid
S_mesh, T_mesh = np.meshgrid(S_range, T_range)
option_price = np.zeros_like(S_mesh)

# Calculate option price for each combination
for i in range(len(T_range)):
    for j in range(len(S_range)):
        option_price[i, j] = black_scholes(S_mesh[i, j], K, T_mesh[i, j], r, sigma, 'call')

# Plot the 3D surface
fig = plt.figure(figsize=(12, 8))
ax = fig.add_subplot(111, projection='3d')
surf = ax.plot_surface(S_mesh, T_mesh, option_price, cmap=cm.coolwarm, alpha=0.8)

# Add labels
ax.set_xlabel('Stock Price ($)')
ax.set_ylabel('Time to Maturity (years)')
ax.set_zlabel('Call Option Price ($)')
ax.set_title('Black-Scholes Call Option Price')

# Add a color bar
fig.colorbar(surf, ax=ax, shrink=0.5, aspect=5)

plt.savefig('black_scholes_3d.png')
plt.show()

# Example 2: Visualize the Greeks (Delta, Gamma, Theta)
def bs_delta(S, K, T, r, sigma, option='call'):
    """Calculate option Delta"""
    d1 = (np.log(S / K) + (r + 0.5 * sigma**2) * T) / (sigma * np.sqrt(T))
    
    if option == 'call':
        return stats.norm.cdf(d1)
    else:
        return stats.norm.cdf(d1) - 1

def bs_gamma(S, K, T, r, sigma):
    """Calculate option Gamma (same for calls and puts)"""
    d1 = (np.log(S / K) + (r + 0.5 * sigma**2) * T) / (sigma * np.sqrt(T))
    return stats.norm.pdf(d1) / (S * sigma * np.sqrt(T))

def bs_theta(S, K, T, r, sigma, option='call'):
    """Calculate option Theta"""
    d1 = (np.log(S / K) + (r + 0.5 * sigma**2) * T) / (sigma * np.sqrt(T))
    d2 = d1 - sigma * np.sqrt(T)
    
    if option == 'call':
        theta = -S * stats.norm.pdf(d1) * sigma / (2 * np.sqrt(T)) - r * K * np.exp(-r * T) * stats.norm.cdf(d2)
    else:
        theta = -S * stats.norm.pdf(d1) * sigma / (2 * np.sqrt(T)) + r * K * np.exp(-r * T) * stats.norm.cdf(-d2)
        
    return theta / 365  # Convert to daily

# Calculate Greeks for various stock prices
S_range = np.linspace(80, 120, 100)
T = 0.5  # 6 months

delta_call = [bs_delta(S, K, T, r, sigma, 'call') for S in S_range]
delta_put = [bs_delta(S, K, T, r, sigma, 'put') for S in S_range]
gamma = [bs_gamma(S, K, T, r, sigma) for S in S_range]
theta_call = [bs_theta(S, K, T, r, sigma, 'call') for S in S_range]
theta_put = [bs_theta(S, K, T, r, sigma, 'put') for S in S_range]

# Plot the Greeks
fig, axs = plt.subplots(3, 1, figsize=(12, 15))

# Delta
axs[0].plot(S_range, delta_call, 'b-', label='Call Delta')
axs[0].plot(S_range, delta_put, 'r-', label='Put Delta')
axs[0].axhline(y=0, color='k', linestyle='-', alpha=0.3)
axs[0].axvline(x=K, color='k', linestyle='--', alpha=0.3)
axs[0].set_xlabel('Stock Price ($)')
axs[0].set_ylabel('Delta')
axs[0].set_title('Option Delta vs Stock Price')
axs[0].grid(True)
axs[0].legend()

# Gamma
axs[1].plot(S_range, gamma, 'g-')
axs[1].axhline(y=0, color='k', linestyle='-', alpha=0.3)
axs[1].axvline(x=K, color='k', linestyle='--', alpha=0.3)
axs[1].set_xlabel('Stock Price ($)')
axs[1].set_ylabel('Gamma')
axs[1].set_title('Option Gamma vs Stock Price')
axs[1].grid(True)

# Theta
axs[2].plot(S_range, theta_call, 'b-', label='Call Theta')
axs[2].plot(S_range, theta_put, 'r-', label='Put Theta')
axs[2].axhline(y=0, color='k', linestyle='-', alpha=0.3)
axs[2].axvline(x=K, color='k', linestyle='--', alpha=0.3)
axs[2].set_xlabel('Stock Price ($)')
axs[2].set_ylabel('Theta ($ per day)')
axs[2].set_title('Option Theta vs Stock Price')
axs[2].grid(True)
axs[2].legend()

plt.tight_layout()
plt.savefig('option_greeks.png')
plt.show()

# Example 3: Monte Carlo simulation for option pricing
def monte_carlo_option_pricing(S0, K, T, r, sigma, num_simulations=10000, option='call'):
    """
    Price options using Monte Carlo simulation
    
    Parameters:
    S0: Initial stock price
    K: Strike price
    T: Time to maturity (in years)
    r: Risk-free interest rate
    sigma: Volatility
    num_simulations: Number of simulation paths
    option: 'call' or 'put'
    
    Returns:
    Option price
    """
    # Generate random simulations of stock price at expiration
    z = np.random.standard_normal(num_simulations)
    ST = S0 * np.exp((r - 0.5 * sigma**2) * T + sigma * np.sqrt(T) * z)
    
    # Calculate payoff at expiration
    if option == 'call':
        payoff = np.maximum(ST - K, 0)
    else:
        payoff = np.maximum(K - ST, 0)
    
    # Calculate present value
    option_price = np.exp(-r * T) * np.mean(payoff)
    
    return option_price, ST

# Run Monte Carlo simulation
np.random.seed(42)  # For reproducibility
S0 = 100
num_simulations = 10000

mc_price, ST = monte_carlo_option_pricing(S0, K, T, r, sigma, num_simulations, 'call')
bs_price = black_scholes(S0, K, T, r, sigma, 'call')

# Create histogram of simulated stock prices at expiration
plt.figure(figsize=(12, 6))
plt.hist(ST, bins=50, alpha=0.5, color='blue', density=True)

# Add a normal distribution curve for comparison
x = np.linspace(min(ST), max(ST), 100)
mu = S0 * np.exp(r * T)  # Expected value under risk-neutral measure
std = S0 * np.exp(r * T) * np.sqrt(np.exp(sigma**2 * T) - 1)
plt.plot(x, stats.norm.pdf(x, mu, std), 'r-', linewidth=2)

plt.axvline(x=K, color='k', linestyle='--', label=f'Strike Price (K={K})')
plt.axvline(x=S0, color='g', linestyle='-', label=f'Initial Price (S0={S0})')

plt.title(f'Monte Carlo Simulation of Stock Prices at Expiration (T={T} years)')
plt.xlabel('Stock Price at Expiration ($)')
plt.ylabel('Probability Density')
plt.legend()
plt.grid(True)
plt.savefig('monte_carlo_simulation.png')
plt.show()

# Compare Monte Carlo with Black-Scholes
print(f"Monte Carlo Call Option Price: ${mc_price:.4f}")
print(f"Black-Scholes Call Option Price: ${bs_price:.4f}")
print(f"Difference: ${abs(mc_price - bs_price):.4f}")

Black-Scholes Option Pricing Model in Financial Mathematics

In this example, I’ll solve and analyze a classic financial mathematics problem: option pricing using the Black-Scholes model.

The Black-Scholes model is fundamental in financial mathematics for pricing European-style options.

Mathematical Background

The Black-Scholes formula for European call and put options is given by:

For a call option:
$$C(S,t) = S \cdot N(d_1) - K e^{-r(T-t)} \cdot N(d_2)$$

For a put option:
$$P(S,t) = K e^{-r(T-t)} \cdot N(-d_2) - S \cdot N(-d_1)$$

Where:
$$d_1 = \frac{\ln(S/K) + (r + \sigma^2/2)(T-t)}{\sigma\sqrt{T-t}}$$
$$d_2 = d_1 - \sigma\sqrt{T-t}$$

And:

$S$ is the current stock price
$K$ is the strike price
$r$ is the risk-free interest rate
$T-t$ is the time to maturity
$\sigma$ is the volatility of the stock
$N(·)$ is the cumulative distribution function of the standard normal distribution

Code Explanation

The code implements the Black-Scholes model and provides three detailed examples:

3D Visualization of Option Prices: Shows how call option prices change with stock price and time to maturity
Option Greeks Analysis: Calculates and visualizes Delta, Gamma, and Theta
Monte Carlo Simulation: Prices options using simulation and compares with the analytical solution

Key Functions:

black_scholes(): Implements the analytical Black-Scholes formula
bs_delta(), bs_gamma(), bs_theta(): Calculate option Greeks
monte_carlo_option_pricing(): Simulates stock price paths using geometric Brownian motion

Visualization Analysis

1. 3D Surface Plot of Call Option Prices

The 3D plot shows how the call option price (z-axis) varies with:

Stock price (x-axis): As stock price increases, call option value increases
Time to maturity (y-axis): The time value of the option is visible

This visualization helps understand the non-linear relationship between these variables.

2. Option Greeks Visualization

The Greeks represent sensitivities of option prices to various factors:

Delta: Measures the rate of change of option price with respect to changes in the underlying asset’s price
- Call delta ranges from $0$ to $1$
- Put delta ranges from $-1$ to $0$
- At-the-money options have delta near $0.5$ (calls) or $-0.5$ (puts)
Gamma: Measures the rate of change of delta with respect to changes in the underlying price
- Highest at-the-money and decreases as the option moves in or out of the money
- Same for both calls and puts
Theta: Measures the rate of change of option price with respect to the passage of time (time decay)
- Generally negative for both calls and puts (options lose value as time passes)
- Most significant for at-the-money options

3. Monte Carlo Simulation

Monte Carlo Call Option Price: $6.8874
Black-Scholes Call Option Price: $6.8887
Difference: $0.0013

The histogram shows the distribution of simulated stock prices at expiration:

The red curve represents the theoretical lognormal distribution
The vertical lines mark the initial stock price and strike price
The simulation validates the Black-Scholes model by producing a similar price

When we compare the Monte Carlo estimate with the analytical Black-Scholes solution, they should be very close, with differences attributable to simulation error.

Conclusion

This example demonstrates several fundamental concepts in financial mathematics:

Analytical solution of the Black-Scholes equation
Option Greeks for risk management
Monte Carlo methods for pricing financial instruments

The visualizations help develop intuition about option pricing behavior and the relationships between different variables in the model.

Visualizing the Prime Number Theorem with Python：A Comparison of $ \pi(x) $ and $ \frac{x}{\log x} $

March 20, 2025

Here’s an example from Analytic Number Theory with a complete solution and visualization in $Python$.

📚 Problem: Prime Number Theorem Approximation

The Prime Number Theorem (PNT) states that:

$$
\pi(x) \sim \frac{x}{\log x}
$$

Where:

$( \pi(x) )$ is the number of primes less than or equal to $( x )$.
$( \frac{x}{\log x} )$ approximates $( \pi(x) )$ for large values of $( x )$.

🎯 Goal:

Calculate $( \pi(x) )$ for values up to $( x = 100,000 )$.
Compare it with the approximation $( \frac{x}{\log x} )$.
Visualize the results using a graph.

📝 Step 1: Define and Calculate Functions

pi(x): Count primes less than or equal to $( x )$.
Approximation: Use $( \frac{x}{\log x} )$.

🧩 Step 2: Python Code Implementation

# Import required libraries
import numpy as np
import matplotlib.pyplot as plt
from sympy import primepi
from math import log

# Define the range and step size
x_values = np.linspace(10, 100000, 200)  # Range from 10 to 100,000 with 200 points

# Calculate prime counting function pi(x)
pi_values = [primepi(int(x)) for x in x_values]

# Calculate the Prime Number Theorem (PNT) approximation
pnt_approx = [x / log(x) for x in x_values]

# Plot the results
plt.figure(figsize=(10, 6))
plt.plot(x_values, pi_values, label=r'$\pi(x)$ (Prime Count)', color='blue')
plt.plot(x_values, pnt_approx, label=r'$\frac{x}{\log x}$ (PNT Approximation)', linestyle='--', color='red')

# Add title and labels
plt.title("Prime Number Theorem: $\pi(x)$ vs. $\dfrac{x}{\log x}$", fontsize=14)
plt.xlabel("x", fontsize=12)
plt.ylabel("Value", fontsize=12)
plt.legend()
plt.grid(True)

# Show the plot
plt.show()

📊 Step 3: Explanation and Visualization

✅ Explanation:

pi(x): Using primepi() from sympy to compute the exact number of primes up to $( x )$.
PNT Approximation: Using the formula $( \frac{x}{\log x} )$.
Graph:
- Blue line: Exact prime counting function.
- Red dashed line: PNT approximation.

📈 Step 4: Interpretation of Results

The graph shows that the approximation $( \frac{x}{\log x} )$ becomes increasingly accurate as $( x )$ increases.
For smaller values of $( x )$, there is a noticeable gap, but for large $( x )$, the approximation closely tracks $( \pi(x) )$.

📢 Insights:

The Prime Number Theorem provides a remarkably simple yet powerful approximation for the distribution of prime numbers.
This visualization helps highlight how the error margin decreases with larger $( x )$.

Visualization and Analysis of Matrix Operations

March 19, 2025

Linear Transformations in Python

I’ll provide an example problem about linear transformations, solve it using $Python$, and visualize the results with graphs.

Here’s a comprehensive example:

Linear Transformation Example and Solution with Python

Let’s consider a linear transformation $T: \mathbb{R}^2 \rightarrow \mathbb{R}^2$ defined by the matrix:

$$T = \begin{pmatrix} 2 & 1 \ 1 & 3 \end{pmatrix}$$

We’ll solve the following problems:

Find the transformed coordinates of several points
Visualize how this transformation affects a unit square
Compute the eigenvalues and eigenvectors of the transformation
Visualize the effect of the transformation on the eigenvectors

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.patches import Polygon
from mpl_toolkits.mplot3d import Axes3D

# Define the linear transformation matrix
T = np.array([[2, 1], 
              [1, 3]])

print("Linear Transformation Matrix T:")
print(T)

# Part 1: Transform some example points
points = np.array([[1, 0], [0, 1], [1, 1], [-1, 2]])
transformed_points = np.dot(points, T)

print("\nOriginal Points:")
for p in points:
    print(f"[{p[0]}, {p[1]}]")
    
print("\nTransformed Points:")
for p in transformed_points:
    print(f"[{p[0]}, {p[1]}]")

# Part 2: Visualize the effect on a unit square
def plot_transformation(T):
    # Create a figure with two subplots
    fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12, 5))
    
    # Define the unit square
    square = np.array([[0, 0], [1, 0], [1, 1], [0, 1], [0, 0]])
    
    # Transform the square
    transformed_square = np.dot(square, T)
    
    # Plot the original square
    ax1.plot(square[:, 0], square[:, 1], 'b-')
    ax1.fill(square[:, 0], square[:, 1], 'lightblue', alpha=0.5)
    ax1.set_xlim(-0.5, 3.5)
    ax1.set_ylim(-0.5, 3.5)
    ax1.grid(True)
    ax1.set_title('Original Unit Square')
    ax1.set_aspect('equal')
    
    # Plot the transformed square
    ax2.plot(transformed_square[:, 0], transformed_square[:, 1], 'r-')
    ax2.fill(transformed_square[:, 0], transformed_square[:, 1], 'lightcoral', alpha=0.5)
    ax2.set_xlim(-0.5, 3.5)
    ax2.set_ylim(-0.5, 3.5)
    ax2.grid(True)
    ax2.set_title('Transformed Square')
    ax2.set_aspect('equal')
    
    plt.tight_layout()
    plt.show()

# Visualize the transformation
plot_transformation(T)

# Part 3: Find eigenvalues and eigenvectors
eigenvalues, eigenvectors = np.linalg.eig(T)

print("\nEigenvalues:")
for i, val in enumerate(eigenvalues):
    print(f"λ_{i+1} = {val:.4f}")

print("\nEigenvectors (as columns):")
print(eigenvectors)

# Part 4: Visualize eigenvectors and their transformations
def plot_eigenvectors(T, eigenvalues, eigenvectors):
    # Create a figure
    fig, ax = plt.subplots(figsize=(8, 8))
    
    # Set limits and grid
    ax.set_xlim(-3, 3)
    ax.set_ylim(-3, 3)
    ax.grid(True)
    ax.axhline(y=0, color='k', linestyle='-', alpha=0.3)
    ax.axvline(x=0, color='k', linestyle='-', alpha=0.3)
    
    # Colors for the eigenvectors
    colors = ['blue', 'green']
    
    # Plot eigenvectors and their transformations
    for i in range(len(eigenvalues)):
        # Get the i-th eigenvector
        v = eigenvectors[:, i]
        
        # Scale eigenvector for better visualization
        v_scaled = v * 2
        
        # Plot the original eigenvector
        ax.arrow(0, 0, v_scaled[0], v_scaled[1], head_width=0.1, head_length=0.1, 
                 fc=colors[i], ec=colors[i], label=f"Eigenvector {i+1}")
        
        # Transform the eigenvector
        tv = np.dot(T, v)
        tv_scaled = tv * 2
        
        # Plot the transformed eigenvector
        ax.arrow(0, 0, tv_scaled[0], tv_scaled[1], head_width=0.1, head_length=0.1, 
                 fc='lightcoral', ec='lightcoral', linestyle='--')
        
        # Add text for eigenvalue
        ax.text(v_scaled[0]*1.1, v_scaled[1]*1.1, 
                f"λ_{i+1}={eigenvalues[i]:.2f}", 
                color=colors[i])
    
    # Add a unit circle and its transformation
    theta = np.linspace(0, 2*np.pi, 100)
    circle_x = np.cos(theta)
    circle_y = np.sin(theta)
    
    # Stack coordinates to create points
    circle_points = np.column_stack((circle_x, circle_y))
    
    # Transform the circle
    transformed_circle = np.dot(circle_points, T)
    
    # Plot the unit circle
    ax.plot(circle_x, circle_y, 'k-', alpha=0.3)
    
    # Plot the transformed circle
    ax.plot(transformed_circle[:, 0], transformed_circle[:, 1], 'r-', alpha=0.3)
    
    ax.set_title('Eigenvectors and Their Transformations')
    ax.legend()
    plt.axis('equal')
    plt.tight_layout()
    plt.show()

# Visualize eigenvectors
plot_eigenvectors(T, eigenvalues, eigenvectors)

# Part 5: Visualize the transformation in 3D
def plot_transformation_3d():
    fig = plt.figure(figsize=(10, 8))
    ax = fig.add_subplot(111, projection='3d')
    
    # Create a grid of points
    x = np.linspace(-2, 2, 10)
    y = np.linspace(-2, 2, 10)
    X, Y = np.meshgrid(x, y)
    
    # Flatten the grid points
    points = np.column_stack((X.flatten(), Y.flatten()))
    
    # Transform the points
    transformed_points = np.dot(points, T)
    
    # Reshape the transformed coordinates
    X_transformed = transformed_points[:, 0].reshape(X.shape)
    Y_transformed = transformed_points[:, 1].reshape(Y.shape)
    
    # Create Z coordinates (all zeros for 2D transformation)
    Z = np.zeros(X.shape)
    Z_transformed = np.ones(X.shape)  # Offset for visualization
    
    # Plot the original grid
    ax.plot_surface(X, Y, Z, alpha=0.5, color='blue', label='Original')
    
    # Plot the transformed grid
    ax.plot_surface(X_transformed, Y_transformed, Z_transformed, alpha=0.5, color='red', label='Transformed')
    
    # Add lines connecting original and transformed points
    for i in range(len(points)):
        ax.plot([points[i, 0], transformed_points[i, 0]],
                [points[i, 1], transformed_points[i, 1]],
                [0, 1], 'k-', alpha=0.2)
    
    ax.set_xlabel('X')
    ax.set_ylabel('Y')
    ax.set_zlabel('Z (Visualization Only)')
    ax.set_title('3D Visualization of 2D Linear Transformation')
    
    plt.tight_layout()
    plt.show()

# Visualize the transformation in 3D
plot_transformation_3d()

Explanation of the Solution

1. Defining the Linear Transformation

We define a linear transformation $T$ using a $2 \times 2$ matrix:

$$T = \begin{pmatrix} 2 & 1 \ 1 & 3 \end{pmatrix}$$

This matrix represents the transformation, and we can apply it to any vector $\vec{v} \in \mathbb{R}^2$ by matrix multiplication: $T\vec{v}$.

2. Transforming Points

When we apply the transformation to specific points, we get:

The point $(1, 0)$ transforms to $(2, 1)$
The point $(0, 1)$ transforms to $(1, 3)$
The point $(1, 1)$ transforms to $(3, 4)$
The point $(-1, 2)$ transforms to $(1, 5)$

This demonstrates how the linear transformation maps points from the input space to the output space.

3. Visualization of the Transformation

Linear Transformation Matrix T:
[[2 1]
 [1 3]]

Original Points:
[1, 0]
[0, 1]
[1, 1]
[-1, 2]

Transformed Points:
[2, 1]
[1, 3]
[3, 4]
[0, 5]

The code visualizes how the unit square (with vertices at $(0,0)$, $(1,0)$, $(1,1)$, and $(0,1)$) is transformed by $T$.

As you can see in the first visualization, the square is stretched, rotated, and skewed by the transformation.

4. Eigenvalues and Eigenvectors

Eigenvalues:
λ_1 = 1.3820
λ_2 = 3.6180

Eigenvectors (as columns):
[[-0.85065081 -0.52573111]
 [ 0.52573111 -0.85065081]]

The eigenvalues and eigenvectors of the transformation matrix tell us about the fundamental behavior of the transformation:

The eigenvalues are approximately $\lambda_1 \approx 1.38$ and $\lambda_2 \approx 3.62$
The corresponding eigenvectors are shown in the output

An eigenvector is a special vector that, when transformed by $T$, only gets scaled by its corresponding eigenvalue without changing direction.

The second visualization shows these eigenvectors and how they’re affected by the transformation.

5. 3D Visualization

The final visualization provides a 3D perspective of the transformation, showing how a grid of points in the original space maps to points in the transformed space.

Key Insights

Direction of Maximum Stretching:
The eigenvector corresponding to the larger eigenvalue ($\lambda_2 \approx 3.62$) indicates the direction in which the transformation stretches vectors the most.
Shape Distortion:
The transformation turns the square into a parallelogram, demonstrating how linear transformations preserve straight lines but can alter angles and distances.
Area Scaling:
The determinant of the transformation matrix ($\det(T) = 5$) tells us that the transformation scales areas by a factor of 5, which you can observe in the increased size of the transformed square.
Eigenvector Behavior:
Note how the eigenvectors are only scaled (not rotated) when the transformation is applied.
This property makes eigenvectors particularly useful in understanding linear transformations.

This example demonstrates the fundamental concepts of linear transformations and how we can visualize and analyze them using Python’s numerical and visualization libraries.

Exploring Nonlinear Dynamics

March 18, 2025

The Logistic Map Bifurcation Diagram in Python

Let’s dive into an example problem in nonlinear dynamics and solve it using Python in Google Colaboratory.

We’ll consider the logistic map, a classic example of a nonlinear dynamical system that exhibits rich behavior, including chaos, depending on its parameter values.

The logistic map is defined by the recursive equation:

$$
x_{n+1} = r x_n (1 - x_n)
$$

where:

$(x_n)$ is the state variable at iteration $(n)$, constrained between $0$ and $1$ (e.g., representing a normalized population),
$(r)$ is a growth rate parameter, typically between $0$ and $4$,
The nonlinearity arises from the quadratic term $(x_n (1 - x_n))$.

Problem Statement

We’ll simulate the logistic map for different values of $(r)$ and visualize how the system transitions from stable fixed points to periodic orbits and eventually to chaotic behavior. Specifically, we’ll:

Iterate the map for a range of $(r)$ values (e.g., $2.5$ to $4.0$),
Plot the long-term behavior (steady states or oscillations) as a bifurcation diagram,
Use Python to compute and graph the results.

Python Code

Here’s the complete code to simulate and visualize the logistic map’s bifurcation diagram:

import numpy as np
import matplotlib.pyplot as plt

# Define the logistic map function
def logistic_map(r, x):
    return r * x * (1 - x)

# Parameters
n_iterations = 1000  # Total iterations
n_transient = 800    # Discard initial transient iterations
n_points = 200       # Number of points to plot per r value
r_values = np.linspace(2.5, 4.0, 1000)  # Range of r values
x0 = 0.1  # Initial condition (must be between 0 and 1)

# Arrays to store results
r_plot = []
x_plot = []

# Iterate over r values
for r in r_values:
    x = x0
    # Run iterations, discard transients
    for i in range(n_iterations):
        x = logistic_map(r, x)
        # Store the last n_points after transient
        if i >= n_transient:
            r_plot.append(r)
            x_plot.append(x)

# Create the bifurcation diagram
plt.figure(figsize=(10, 6))
plt.plot(r_plot, x_plot, ',k', alpha=0.25)  # Plot as small black dots
plt.xlabel('$r$')
plt.ylabel('$x$')
plt.title('Bifurcation Diagram of the Logistic Map')
plt.grid(True)
plt.show()

Code Explanation

Imports:
- numpy for numerical computations (e.g., creating arrays of $(r)$ values).
- matplotlib.pyplot for plotting the bifurcation diagram.
Logistic Map Function:
- logistic_map(r, x) implements the equation $(x_{n+1} = r x_n (1 - x_n))$.
Parameters:
- n_iterations: Total number of iterations per $(r)$ value.
- n_transient: Number of initial iterations to discard (to focus on steady-state behavior).
- n_points: Number of points to plot after transients.
- r_values: Array of $(r)$ from $2.5$ to $4.0$ with $1000$ steps.
- x0: Initial condition set to $0.1$ (any value between $0$ and $1$ works).
Simulation:
- For each $(r)$, start with $(x = x0)$.
- Iterate the map n_iterations times.
- After n_transient iterations, store the subsequent $(x)$ values and corresponding $(r)$ values.
Plotting:
- Use plt.plot with a comma marker (,) to plot points as tiny dots, creating a dense bifurcation diagram.
- Add labels and a title with TEX-formatted variables (e.g., $r$ ).

Results and Graph Interpretation

When you run this code in Google Colaboratory, it generates a bifurcation diagram. Here’s what you’ll see:

For $(r < 3)$: The system converges to a single fixed point (e.g., at $(r = 2.5)$, $(x)$ stabilizes around $0.6$). This is a stable equilibrium.
For $(3 < r < 3.57)$: Period-doubling occurs. At $(r = 3.2)$, you’ll see two values (a period-2 cycle), then four (period-4), and so on.
For $(r > 3.57)$: The system enters chaotic regimes, interspersed with periodic windows (e.g., a clear period-3 window around $(r = 3.83)$).
At $(r = 4)$: Fully chaotic behavior, with $(x)$ values filling the interval [0, 1] densely.

The graph shows $(x)$ versus $(r)$.
Initially, you see a single line (fixed point), which splits into two, four, etc., before becoming a scattered cloud of points in the chaotic region.
This beautifully illustrates how a simple nonlinear equation can produce complex dynamics.

Conclusion

This example demonstrates key features of nonlinear dynamics: stability, bifurcations, and chaos.

The Python code is efficient and leverages NumPy for fast computation and Matplotlib for clear visualization.

You can tweak parameters like r_values or n_iterations to explore further!

Geodesic Trajectories on a Torus

March 17, 2025

Computational Analysis and Visualization in Differential Geometry

I’ll provide an example problem in differential geometry, solve it using $Python$, and visualize the results with graphs.

I’ll include both the code and mathematical expressions.

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp
from mpl_toolkits.mplot3d import Axes3D
from matplotlib import cm

# Define the torus parametrization
def torus(u, v, R=2, r=1):
    """
    Parametrize a torus with major radius R and minor radius r
    u, v are the parameters in [0, 2π)
    """
    x = (R + r * np.cos(v)) * np.cos(u)
    y = (R + r * np.cos(v)) * np.sin(u)
    z = r * np.sin(v)
    return x, y, z

# Compute the metric coefficients for the torus
def torus_metric(u, v, R=2, r=1):
    """
    Compute the metric coefficients E, F, G for the torus at point (u, v)
    """
    E = (R + r * np.cos(v))**2  # g_11
    F = 0                       # g_12 = g_21
    G = r**2                    # g_22
    return E, F, G

# Compute the Christoffel symbols for the torus
def christoffel_symbols(u, v, R=2, r=1):
    """
    Compute the Christoffel symbols for the torus at point (u, v)
    Returns a 2x2x2 array where Gamma[i,j,k] = Gamma^i_{jk}
    """
    E, F, G = torus_metric(u, v, R, r)
    
    # Derivatives of metric coefficients
    dE_du = 0
    dE_dv = -2 * (R + r * np.cos(v)) * r * np.sin(v)
    dF_du = dF_dv = 0
    dG_du = 0
    dG_dv = 0
    
    # Initialize Christoffel symbols
    Gamma = np.zeros((2, 2, 2))
    
    # g^{11} = 1/E, g^{12} = g^{21} = 0, g^{22} = 1/G
    g_inv_11 = 1/E
    g_inv_22 = 1/G
    
    # Gamma^1_{11}
    Gamma[0, 0, 0] = 0.5 * g_inv_11 * dE_du
    
    # Gamma^1_{12} = Gamma^1_{21}
    Gamma[0, 0, 1] = Gamma[0, 1, 0] = 0.5 * g_inv_11 * dE_dv
    
    # Gamma^1_{22}
    Gamma[0, 1, 1] = -0.5 * g_inv_11 * dG_du
    
    # Gamma^2_{11}
    Gamma[1, 0, 0] = -0.5 * g_inv_22 * dE_dv
    
    # Gamma^2_{12} = Gamma^2_{21}
    Gamma[1, 0, 1] = Gamma[1, 1, 0] = 0.5 * g_inv_22 * dF_du
    
    # Gamma^2_{22}
    Gamma[1, 1, 1] = 0.5 * g_inv_22 * dG_dv
    
    return Gamma

# Geodesic equation as a first-order system
def geodesic_equation(t, y, R=2, r=1):
    """
    The geodesic equation as a first-order system
    y = [u, v, u_dot, v_dot]
    Returns [u_dot, v_dot, u_ddot, v_ddot]
    """
    u, v, u_dot, v_dot = y
    
    # Compute Christoffel symbols
    Gamma = christoffel_symbols(u, v, R, r)
    
    # Geodesic equations (second derivatives)
    u_ddot = -Gamma[0, 0, 0] * u_dot**2 - 2 * Gamma[0, 0, 1] * u_dot * v_dot - Gamma[0, 1, 1] * v_dot**2
    v_ddot = -Gamma[1, 0, 0] * u_dot**2 - 2 * Gamma[1, 0, 1] * u_dot * v_dot - Gamma[1, 1, 1] * v_dot**2
    
    return [u_dot, v_dot, u_ddot, v_ddot]

# Solve the geodesic equation with given initial conditions
def solve_geodesic(u0, v0, u_dot0, v_dot0, t_span, R=2, r=1):
    """
    Solve the geodesic equation with initial conditions
    u0, v0: initial position
    u_dot0, v_dot0: initial velocity
    t_span: time span [t_start, t_end]
    """
    y0 = [u0, v0, u_dot0, v_dot0]
    
    # Solve the ODE system
    sol = solve_ivp(
        lambda t, y: geodesic_equation(t, y, R, r),
        t_span,
        y0,
        method='RK45',
        rtol=1e-8,
        atol=1e-8,
        dense_output=True
    )
    
    return sol

# Plot the torus and geodesic
def plot_torus_and_geodesic(sol, R=2, r=1, num_points=1000):
    """
    Plot the torus surface and the geodesic curve
    """
    # Create figure
    fig = plt.figure(figsize=(12, 10))
    ax = fig.add_subplot(111, projection='3d')
    
    # Create a mesh grid for the torus
    u = np.linspace(0, 2*np.pi, 100)
    v = np.linspace(0, 2*np.pi, 100)
    u_grid, v_grid = np.meshgrid(u, v)
    
    # Compute the torus coordinates
    x, y, z = torus(u_grid, v_grid, R, r)
    
    # Plot the torus surface with translucent appearance
    ax.plot_surface(x, y, z, cmap=cm.viridis, alpha=0.6, edgecolor='none')
    
    # Extract the solution at evenly spaced points
    t = np.linspace(sol.t[0], sol.t[-1], num_points)
    y_interp = sol.sol(t)
    u_sol, v_sol = y_interp[0], y_interp[1]
    
    # Convert the geodesic curve to Cartesian coordinates
    x_geo, y_geo, z_geo = torus(u_sol, v_sol, R, r)
    
    # Plot the geodesic curve
    ax.plot(x_geo, y_geo, z_geo, 'r-', linewidth=2, label='Geodesic')
    
    # Mark the start and end points
    ax.plot([x_geo[0]], [y_geo[0]], [z_geo[0]], 'go', markersize=8, label='Start')
    ax.plot([x_geo[-1]], [y_geo[-1]], [z_geo[-1]], 'bo', markersize=8, label='End')
    
    # Set labels and title
    ax.set_xlabel('X')
    ax.set_ylabel('Y')
    ax.set_zlabel('Z')
    ax.set_title('Geodesic on a Torus (R={}, r={})'.format(R, r))
    ax.legend()
    
    plt.tight_layout()
    plt.show()
    
    return fig, ax

# Analyze the length and curvature of the geodesic
def analyze_geodesic(sol, R=2, r=1, num_points=1000):
    """
    Analyze properties of the geodesic such as length and curvature
    """
    # Extract the solution at evenly spaced points
    t = np.linspace(sol.t[0], sol.t[-1], num_points)
    y_interp = sol.sol(t)
    u_sol, v_sol = y_interp[0], y_interp[1]
    u_dot_sol, v_dot_sol = y_interp[2], y_interp[3]
    
    # Compute the metric coefficients along the geodesic
    E, F, G = [], [], []
    for i in range(len(t)):
        E_i, F_i, G_i = torus_metric(u_sol[i], v_sol[i], R, r)
        E.append(E_i)
        F.append(F_i)
        G.append(G_i)
    
    E = np.array(E)
    F = np.array(F)
    G = np.array(G)
    
    # Compute the speed along the geodesic
    speed = np.sqrt(E * u_dot_sol**2 + 2*F * u_dot_sol * v_dot_sol + G * v_dot_sol**2)
    
    # Compute the length of the geodesic
    dt = t[1] - t[0]
    length = np.sum(speed * dt)
    
    # Plot the speed along the geodesic
    plt.figure(figsize=(10, 6))
    plt.plot(t, speed)
    plt.title('Speed along the Geodesic')
    plt.xlabel('t')
    plt.ylabel('Speed')
    plt.grid(True)
    plt.show()
    
    print(f"Geodesic Length: {length:.6f}")
    
    return length

# Main execution code
if __name__ == "__main__":
    # Parameters for the torus
    R = 2.0  # Major radius
    r = 1.0  # Minor radius
    
    # Initial conditions for the geodesic
    u0 = 0.0
    v0 = 0.0
    u_dot0 = 1.0
    v_dot0 = 0.5
    
    # Time span for integration
    t_span = [0, 10]
    
    # Solve the geodesic equation
    sol = solve_geodesic(u0, v0, u_dot0, v_dot0, t_span, R, r)
    
    # Plot the results
    plot_torus_and_geodesic(sol, R, r)
    
    # Analyze the geodesic
    geodesic_length = analyze_geodesic(sol, R, r)
    
    # Plot a few different geodesics
    plt.figure(figsize=(12, 10))
    
    initial_conditions = [
        (0.0, 0.0, 1.0, 0.0, 'r', 'Horizontal'),
        (0.0, 0.0, 0.0, 1.0, 'g', 'Vertical'),
        (0.0, 0.0, 1.0, 0.5, 'b', 'Diagonal'),
        (0.0, 0.0, 1.0, 2.0, 'purple', 'Steep Diagonal'),
    ]
    
    ax = plt.axes(projection='3d')
    
    # Plot the torus wireframe
    u = np.linspace(0, 2*np.pi, 30)
    v = np.linspace(0, 2*np.pi, 30)
    u_grid, v_grid = np.meshgrid(u, v)
    x, y, z = torus(u_grid, v_grid, R, r)
    ax.plot_surface(x, y, z, color='lightgray', alpha=0.2)
    
    # Compute and plot different geodesics
    for u0, v0, u_dot0, v_dot0, color, label in initial_conditions:
        sol = solve_geodesic(u0, v0, u_dot0, v_dot0, t_span, R, r)
        
        # Extract the solution
        t = np.linspace(sol.t[0], sol.t[-1], 1000)
        y_interp = sol.sol(t)
        u_sol, v_sol = y_interp[0], y_interp[1]
        
        # Convert to Cartesian coordinates
        x_geo, y_geo, z_geo = torus(u_sol, v_sol, R, r)
        
        # Plot the geodesic
        ax.plot(x_geo, y_geo, z_geo, color=color, linewidth=2, label=label)
    
    ax.set_xlabel('X')
    ax.set_ylabel('Y')
    ax.set_zlabel('Z')
    ax.set_title('Different Geodesics on a Torus (R={}, r={})'.format(R, r))
    ax.legend()
    
    plt.tight_layout()
    plt.show()

Differential Geometry Example: Geodesics on a Torus

In this example, I’ll solve the geodesic equations on a torus surface and visualize the results.

Geodesics are curves that locally minimize distance between points on a manifold, representing the “straightest possible paths” on the surface.

Mathematical Background

Torus Parameterization

A torus with major radius $R$ and minor radius $r$ can be parameterized as:

$$x(u,v) = (R + r\cos v)\cos u$$
$$y(u,v) = (R + r\cos v)\sin u$$
$$z(u,v) = r\sin v$$

where $u, v \in [0, 2\pi)$.

Metric Tensor

The components of the first fundamental form (metric tensor) are:

$$E = g_{11} = (R + r\cos v)^2$$
$$F = g_{12} = g_{21} = 0$$
$$G = g_{22} = r^2$$

Geodesic Equations

The geodesic equations are given by:

$$\frac{d^2u}{dt^2} + \Gamma^1_{11}\left(\frac{du}{dt}\right)^2 + 2\Gamma^1_{12}\frac{du}{dt}\frac{dv}{dt} + \Gamma^1_{22}\left(\frac{dv}{dt}\right)^2 = 0$$

$$\frac{d^2v}{dt^2} + \Gamma^2_{11}\left(\frac{du}{dt}\right)^2 + 2\Gamma^2_{12}\frac{du}{dt}\frac{dv}{dt} + \Gamma^2_{22}\left(\frac{dv}{dt}\right)^2 = 0$$

where $\Gamma^i_{jk}$ are the Christoffel symbols.

Code Explanation

Surface Parameterization: I defined the torus using the standard parameterization with major radius $R$ and minor radius $r$.
Metric Calculation: I computed the components of the metric tensor ($E$, $F$, $G$).
Christoffel Symbols: I calculated the Christoffel symbols from the metric and its derivatives.
Geodesic Equations: I formulated the geodesic equations as a first-order ODE system.
Numerical Integration: I used solve_ivp from $SciPy$ to numerically integrate the geodesic equations. This converts our second-order differential equations into a system of first-order ODEs.
Visualization: I created 3D plots showing the torus surface and the geodesic curves on it.
Analysis: I calculated and plotted the speed along the geodesic and computed its total length.

Key Functions in the Code

torus(u, v, R, r): Computes the Cartesian coordinates (x, y, z) for given parameters (u, v)
torus_metric(u, v, R, r): Calculates the metric tensor components E, F, G
christoffel_symbols(u, v, R, r): Computes the Christoffel symbols for the torus
geodesic_equation(t, y, R, r): Implements the geodesic equations as a first-order system
solve_geodesic(u0, v0, u_dot0, v_dot0, t_span, R, r): Solves the geodesic initial value problem
plot_torus_and_geodesic(sol, R, r, num_points): Visualizes the torus and geodesic
analyze_geodesic(sol, R, r, num_points): Computes geodesic properties like length and speed

Results Interpretation

When you run the code, you’ll see several visualizations:

Main Geodesic Visualization:

Shows a single geodesic curve (in red) on the torus surface, with green and blue markers for the start and end points.
Speed Profile:

A 2D plot showing how the speed of a particle following the geodesic changes over time.
For a geodesic with proper parameterization, this speed should be constant, but numerical integration can introduce small variations.
Multiple Geodesics:

A comparison of different geodesics with various initial velocities:
- Horizontal geodesic ($u_dot0 = 1, v_dot0 = 0$)
- Vertical geodesic ($u_dot0 = 0, v_dot0 = 1$)
- Diagonal geodesic ($u_dot0 = 1, v_dot0 = 0.5$)
- Steep diagonal geodesic ($u_dot0 = 1, v_dot0 = 2.0$)

Mathematical Insights

Horizontal and Vertical Geodesics: Special cases where the geodesics follow circles on the torus.
When $u_dot0 = 1$ and $v_dot0 = 0$, the geodesic follows a circle of constant v.
When $u_dot0 = 0$ and $v_dot0 = 1$, it follows a circle of constant u.
Diagonal Geodesics: More complex paths that wind around the torus, potentially never closing back on themselves (creating dense coverage if extended long enough).
Geodesic Behavior: Notice how geodesics tend to avoid the inner part of the torus and prefer to travel along the outer part where the curvature is less extreme.
This is a direct consequence of the geodesic equations trying to minimize the path length.

Applications

This type of analysis has applications in:

Computer graphics (for path planning on surfaces)
General relativity (geodesics represent the paths of free-falling particles)
Robotics (finding optimal paths on curved surfaces)
Differential geometry (studying the intrinsic properties of surfaces)

You can modify the code to explore different surfaces by changing the parameterization, metric, and Christoffel symbols accordingly.

The Central Limit Theorem

March 16, 2025

A Visual Exploration with Python

I’ll solve a probability theory problem using $Python$, explain the mathematical concepts, and visualize the results with graphs.

Let’s look at a classic probability problem: the Central Limit Theorem demonstration.

import numpy as np
import matplotlib.pyplot as plt
from scipy import stats

# Set random seed for reproducibility
np.random.seed(42)

# Parameters
n_samples = 10000  # Number of samples to generate
sample_sizes = [1, 2, 5, 10, 30, 50]  # Different sample sizes to demonstrate CLT
distributions = {
    'uniform': stats.uniform(0, 1),
    'exponential': stats.expon(scale=1),
    'poisson': stats.poisson(5)
}

# Function to generate samples and calculate means
def generate_sample_means(distribution, sample_size, n_samples):
    # Generate n_samples sets of sample_size samples and calculate their means
    sample_means = np.array([
        np.mean(distribution.rvs(size=sample_size)) 
        for _ in range(n_samples)
    ])
    return sample_means

# Setup the plots
fig, axes = plt.subplots(len(distributions), len(sample_sizes), figsize=(20, 12))
fig.suptitle('Central Limit Theorem Demonstration', fontsize=16)

# Create plots for each distribution and sample size
for i, (dist_name, distribution) in enumerate(distributions.items()):
    # Calculate the theoretical mean and std of the distribution
    if dist_name == 'uniform':
        true_mean = 0.5
        true_std = 1/np.sqrt(12)
    elif dist_name == 'exponential':
        true_mean = 1
        true_std = 1
    else:  # poisson
        true_mean = 5
        true_std = np.sqrt(5)
    
    for j, sample_size in enumerate(sample_sizes):
        # Generate sample means
        means = generate_sample_means(distribution, sample_size, n_samples)
        
        # Calculate expected standard deviation of the sampling distribution
        expected_std = true_std / np.sqrt(sample_size)
        
        # Plot histogram
        ax = axes[i, j]
        ax.hist(means, bins=50, density=True, alpha=0.6, color='skyblue')
        
        # Plot the normal distribution curve for comparison
        x = np.linspace(min(means), max(means), 1000)
        normal_pdf = stats.norm.pdf(x, true_mean, expected_std)
        ax.plot(x, normal_pdf, 'r-', linewidth=2)
        
        # Add theoretical mean line
        ax.axvline(true_mean, color='green', linestyle='--', linewidth=1.5)
        
        # Set title and labels
        if i == 0:
            ax.set_title(f'Sample Size = {sample_size}')
        if j == 0:
            ax.set_ylabel(f'{dist_name.capitalize()}')
        
        # Add some stats
        sample_mean = np.mean(means)
        sample_std = np.std(means)
        ax.text(0.05, 0.95, f'Mean: {sample_mean:.4f}\nSTD: {sample_std:.4f}', 
                transform=ax.transAxes, fontsize=8, 
                verticalalignment='top', bbox=dict(boxstyle='round', facecolor='white', alpha=0.7))

plt.tight_layout()
plt.subplots_adjust(top=0.9)

# Create additional visualization for convergence of variance
plt.figure(figsize=(10, 6))
colors = ['blue', 'red', 'green']
markers = ['o', 's', '^']

for i, (dist_name, distribution) in enumerate(distributions.items()):
    if dist_name == 'uniform':
        true_std = 1/np.sqrt(12)
    elif dist_name == 'exponential':
        true_std = 1
    else:  # poisson
        true_std = np.sqrt(5)
    
    std_ratios = []
    for sample_size in range(1, 51):
        means = generate_sample_means(distribution, sample_size, n_samples)
        sample_std = np.std(means)
        expected_std = true_std / np.sqrt(sample_size)
        std_ratio = sample_std / expected_std
        std_ratios.append(std_ratio)
    
    plt.plot(range(1, 51), std_ratios, color=colors[i], marker=markers[i], 
             linestyle='-', label=f'{dist_name.capitalize()}', markersize=4)

plt.axhline(1, color='black', linestyle='--')
plt.xlabel('Sample Size')
plt.ylabel('Ratio of Sample STD to Expected STD')
plt.title('Convergence of Sampling Distribution Standard Deviation')
plt.legend()
plt.grid(True, alpha=0.3)
plt.show()

# Let's also visualize the original distributions
plt.figure(figsize=(15, 5))
x_values = {
    'uniform': np.linspace(0, 1, 1000),
    'exponential': np.linspace(0, 5, 1000),
    'poisson': np.arange(0, 20)
}

for i, (dist_name, distribution) in enumerate(distributions.items()):
    plt.subplot(1, 3, i+1)
    
    if dist_name == 'poisson':
        # For discrete distribution
        x = x_values[dist_name]
        plt.bar(x, distribution.pmf(x), alpha=0.6, color=colors[i])
    else:
        # For continuous distributions
        x = x_values[dist_name]
        plt.plot(x, distribution.pdf(x), color=colors[i], linewidth=2)
        plt.fill_between(x, distribution.pdf(x), alpha=0.3, color=colors[i])
    
    plt.title(f'{dist_name.capitalize()} Distribution')
    plt.xlabel('Value')
    plt.ylabel('Density' if dist_name != 'poisson' else 'Probability')
    plt.grid(True, alpha=0.3)

plt.tight_layout()
plt.show()

Central Limit Theorem Example

The Central Limit Theorem (CLT) is one of the most important results in probability theory.

It states that when independent random variables are added, their properly normalized sum tends toward a normal distribution even if the original variables are not normally distributed.

Mathematical Formulation

Let $X_1, X_2, \ldots, X_n$ be independent and identically distributed random variables with mean $\mu$ and finite variance $\sigma^2$.

The CLT states that as $n$ approaches infinity, the distribution of:

$$Z_n = \frac{\bar{X}_n - \mu}{\sigma/\sqrt{n}}$$

converges to the standard normal distribution $N(0,1)$,
where is the sample mean.

In other words, the sampling distribution of the mean approaches a normal distribution with mean $\mu$ and standard deviation $\frac{\sigma}{\sqrt{n}}$ as the sample size increases.

Code Explanation

My code demonstrates the Central Limit Theorem by:

Generating samples from three different distributions (uniform, exponential, and Poisson)
Computing the mean of each sample for various sample sizes
Plotting the distribution of these sample means alongside a normal distribution curve
Tracking how quickly the standard deviation of the sample means converges to the theoretical value

Key components:

generate_sample_means() creates multiple samples of a specified size and calculates their means
The first plot shows histograms of sample means for different distributions and sample sizes
The second plot shows how the standard deviation ratio converges to the expected value
The third plot visualizes the original distributions to highlight their different shapes

Theoretical Insights

For each distribution, I’ve calculated the true mean and standard deviation:

Uniform distribution on [0,1]:
- Mean: $\mu = \frac{0+1}{2} = 0.5$
- Variance: $\sigma^2 = \frac{(1-0)^2}{12} = \frac{1}{12}$
Exponential distribution with scale=1:
- Mean: $\mu = \lambda = 1$
- Variance: $\sigma^2 = \lambda^2 = 1$
Poisson distribution with λ=5:
- Mean: $\mu = \lambda = 5$
- Variance: $\sigma^2 = \lambda = 5$

The standard deviation of the sampling distribution should be $\frac{\sigma}{\sqrt{n}}$ where n is the sample size.

Results Interpretation

The graphs demonstrate several important aspects of the CLT:

As sample size increases, the distribution of sample means becomes more normal (bell-shaped), regardless of the original distribution shape
The larger the sample size, the smaller the standard deviation of the sampling distribution
The convergence happens relatively quickly - even with sample sizes of $30$, the approximation is quite good
The standard deviation of the sample means scales proportionally to $\frac{1}{\sqrt{n}}$

The bottom plot showing the ratio of observed to theoretical standard deviation confirms that our experimental results match the theory as it approaches $1.0$ for all distributions as sample size increases.

This demonstrates why the CLT is so powerful in statistics - it allows us to make inferences about population parameters regardless of the underlying distribution, as long as we have a sufficiently large sample size.

Solving the Lotka-Volterra Predator-Prey Model

March 15, 2025

Mathematical Analysis and Python Visualization

I’ll provide an example problem in biomathematics, solve it using $Python$, and include visualizations to illustrate the results.

Let me walk you through a classic example: the Lotka-Volterra predator-prey model.

Lotka-Volterra Predator-Prey Model

The Lotka-Volterra equations describe the dynamics of biological systems in which two species interact, one as a predator and one as prey.

This model is fundamental in mathematical biology.

The basic equations are:

$$\frac{dx}{dt} = \alpha x - \beta xy$$
$$\frac{dy}{dt} = \delta xy - \gamma y$$

Where:

$x$ is the prey population
$y$ is the predator population
$\alpha$ is the prey’s natural growth rate
$\beta$ is the predation rate
$\delta$ is the predator’s reproduction rate (proportional to food intake)
$\gamma$ is the predator’s natural mortality rate

Let’s solve this system of differential equations numerically using $Python$:

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

# Define the Lotka-Volterra model
def lotka_volterra(state, t, alpha, beta, delta, gamma):
    x, y = state
    dx_dt = alpha * x - beta * x * y
    dy_dt = delta * x * y - gamma * y
    return [dx_dt, dy_dt]

# Set parameter values
alpha = 1.0   # Prey growth rate
beta = 0.1    # Predation rate
delta = 0.02  # Conversion efficiency (predator reproduction)
gamma = 0.4   # Predator mortality rate

# Initial conditions
x0 = 40       # Initial prey population
y0 = 9        # Initial predator population
state0 = [x0, y0]

# Time points to solve for
t = np.linspace(0, 100, 1000)

# Solve the differential equations
solution = odeint(lotka_volterra, state0, t, args=(alpha, beta, delta, gamma))
x_solution = solution[:, 0]  # Prey population over time
y_solution = solution[:, 1]  # Predator population over time

# Create a figure with multiple plots
plt.figure(figsize=(16, 12))

# Plot 1: Population vs. Time
plt.subplot(2, 2, 1)
plt.plot(t, x_solution, 'g-', label='Prey (x)')
plt.plot(t, y_solution, 'r-', label='Predator (y)')
plt.grid(True)
plt.xlabel('Time')
plt.ylabel('Population')
plt.title('Population Dynamics Over Time')
plt.legend()

# Plot 2: Phase Plot
plt.subplot(2, 2, 2)
plt.plot(x_solution, y_solution, 'b-')
plt.grid(True)
plt.xlabel('Prey Population (x)')
plt.ylabel('Predator Population (y)')
plt.title('Phase Plot: Predator vs. Prey')

# Mark the initial point
plt.plot(x0, y0, 'ro', label='Initial State')
plt.legend()

# Plot 3: 3D visualization of the dynamics
ax = plt.subplot(2, 2, 3, projection='3d')
ax.plot(x_solution, y_solution, t, color='purple')
ax.set_xlabel('Prey Population (x)')
ax.set_ylabel('Predator Population (y)')
ax.set_zlabel('Time (t)')
ax.set_title('3D Trajectory of the System')

# Plot 4: Vector field showing the dynamics
plt.subplot(2, 2, 4)
x_range = np.linspace(0.5 * min(x_solution), 1.5 * max(x_solution), 20)
y_range = np.linspace(0.5 * min(y_solution), 1.5 * max(y_solution), 20)
X, Y = np.meshgrid(x_range, y_range)

# Calculate derivatives
DX = alpha * X - beta * X * Y
DY = delta * X * Y - gamma * Y

# Normalize the vectors for better visualization
norm = np.sqrt(DX**2 + DY**2)
# Avoid division by zero
norm[norm == 0] = 1
DX = DX / norm
DY = DY / norm

plt.quiver(X, Y, DX, DY, norm, cmap=cm.viridis)
plt.colorbar(label='Vector Magnitude')
plt.xlabel('Prey Population (x)')
plt.ylabel('Predator Population (y)')
plt.title('Vector Field of the System')

# Overlay the solution trajectory
plt.plot(x_solution, y_solution, 'r-', alpha=0.3, label='Trajectory')
plt.legend()

plt.tight_layout()
plt.savefig('lotka_volterra_analysis.png', dpi=300)
plt.show()

# Calculate statistics about the cycles
from scipy.signal import find_peaks

# Find peaks in prey population to analyze cycles
peaks, _ = find_peaks(x_solution, height=np.mean(x_solution))
peak_times = t[peaks]
cycle_periods = np.diff(peak_times)

print(f"Average cycle period: {np.mean(cycle_periods):.2f} time units")
print(f"Maximum prey population: {np.max(x_solution):.2f}")
print(f"Maximum predator population: {np.max(y_solution):.2f}")
print(f"Minimum prey population: {np.min(x_solution):.2f}")
print(f"Minimum predator population: {np.min(y_solution):.2f}")

Code Explanation

Model Definition: I’ve defined the Lotka-Volterra differential equations in a function that returns the rate of change for both populations.
Parameters:
- $α (alpha) = 1.0$: Growth rate of prey in absence of predators
- $β (beta) = 0.1$: Rate at which predators kill prey
- $δ (delta) = 0.02$: Rate at which predators increase by consuming prey
- $γ (gamma) = 0.4$: Natural death rate of predators
Numerical Solution: Using $SciPy$’s odeint to solve the differential equations over time.
Visualization: The code creates four different plots:
- Time series of both populations
- Phase plot showing the predator-prey cycle
- 3D trajectory showing how the system evolves
- Vector field that illustrates the dynamics of the system
Analysis: The code calculates statistics about the cycles, such as average period and population extremes.

Results Interpretation

When you run this code, you’ll observe the classic predator-prey relationship:

Average cycle period: 10.30 time units
Maximum prey population: 40.52
Maximum predator population: 15.95
Minimum prey population: 7.95
Minimum predator population: 5.75

Cyclical Behavior: Both populations oscillate, but out of phase with each other.
When prey increase, predators follow with a delay. More predators reduce the prey population, which then causes predator decline.
Phase Plot: The closed loops in the phase plot indicate periodic behavior - the system returns to similar states over time rather than reaching a stable equilibrium.
Vector Field: The arrows show the direction of change at each point in the state space.
This helps understand how the system would evolve from any initial condition.
Ecological Insights:
- The system exhibits natural oscillations without external forcing
- The predator population peaks after the prey population
- Neither species goes extinct under these parameters
- The system is sensitive to initial conditions and parameter values

This model, despite its simplicity, captures fundamental dynamics of predator-prey systems observed in nature, such as the famous case of lynx and snowshoe hare populations in Canada, which show similar oscillatory patterns.

You can modify the parameters to explore different ecological scenarios, such as:

Higher predation rates (increase $β$)
Lower predator mortality (decrease $γ$)
Different initial population ratios

These changes would result in different dynamics, potentially including extinction of one or both species under extreme parameter values.

Integer Programming：Solving the Knapsack Problem in Python with PuLP

March 14, 2025

I’ll solve an integer programming example using $Python$, showing the implementation, solution, and visualization.

Let’s solve a classic integer programming problem - the Knapsack Problem.

Integer Programming Example: Knapsack Problem

The knapsack problem is a classic optimization problem: Given a set of items, each with a weight and a value, determine which items to include in a collection so that the total weight is less than or equal to a given limit and the total value is maximized.

Mathematically, we can formulate this as:

$$\begin{align}
\text{Maximize} & \sum_{i=1}^{n} v_i x_i \
\text{Subject to} & \sum_{i=1}^{n} w_i x_i \leq W \
& x_i \in {0,1} \quad \forall i \in {1,2,\ldots,n}
\end{align}$$

Where:

$v_i$ is the value of item $i$
$w_i$ is the weight of item $i$
$x_i$ is a binary variable indicating whether item $i$ is selected
$W$ is the maximum weight capacity

Let’s implement this using $PuLP$, a popular linear programming library in $Python$:

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
from pulp import *

# Problem data
items = ["Item_" + str(i) for i in range(1, 11)]
values = [10, 13, 18, 31, 7, 15, 28, 22, 11, 17]
weights = [5, 8, 12, 16, 3, 9, 14, 11, 5, 8]
capacity = 40

# Create a dataframe for better visualization
data = pd.DataFrame({
    'Item': items,
    'Value': values,
    'Weight': weights,
    'Value/Weight': [v/w for v, w in zip(values, weights)]
})

# Create the LP problem
prob = LpProblem("Knapsack_Problem", LpMaximize)

# Create binary variables for each item
x = LpVariable.dicts("Include", items, lowBound=0, upBound=1, cat='Integer')

# Objective function: maximize total value
prob += lpSum([values[i] * x[items[i]] for i in range(len(items))]), "Total Value"

# Constraint: total weight must not exceed capacity
prob += lpSum([weights[i] * x[items[i]] for i in range(len(items))]) <= capacity, "Weight_Constraint"

# Solve the problem
prob.solve()

# Extract results
selected_items = []
selected_values = []
selected_weights = []
total_value = 0
total_weight = 0

for i in range(len(items)):
    if value(x[items[i]]) == 1:
        selected_items.append(items[i])
        selected_values.append(values[i])
        selected_weights.append(weights[i])
        total_value += values[i]
        total_weight += weights[i]

# Create results dataframe
results = pd.DataFrame({
    'Item': selected_items,
    'Value': selected_values,
    'Weight': selected_weights
})

# Print solution information
print("Status:", LpStatus[prob.status])
print("\nSelected Items:")
print(results)
print("\nTotal Value:", total_value)
print("Total Weight:", total_weight)
print("Capacity:", capacity)

# Visualization of the solution
plt.figure(figsize=(12, 10))

# Plot 1: All items with Value vs Weight
plt.subplot(2, 2, 1)
sns.scatterplot(x='Weight', y='Value', data=data, s=100)
for i, row in data.iterrows():
    plt.annotate(row['Item'], (row['Weight'], row['Value']), fontsize=9)
plt.title('All Items: Value vs Weight')
plt.grid(True, alpha=0.3)

# Plot 2: Selected vs Not Selected
plt.subplot(2, 2, 2)
selected_status = ['Selected' if item in selected_items else 'Not Selected' for item in items]
plot_data = data.copy()
plot_data['Status'] = selected_status
sns.scatterplot(x='Weight', y='Value', hue='Status', data=plot_data, s=100)
plt.title('Selected vs Not Selected Items')
plt.grid(True, alpha=0.3)

# Plot 3: Value/Weight Efficiency
plt.subplot(2, 2, 3)
sns.barplot(x='Item', y='Value/Weight', data=data)
plt.xticks(rotation=45)
plt.title('Value/Weight Ratio of Items')
plt.grid(True, alpha=0.3)

# Plot 4: Contribution to Solution
plt.subplot(2, 2, 4)
if not results.empty:
    results['Percentage'] = results['Value'] / total_value * 100
    sns.barplot(x='Item', y='Percentage', data=results)
    plt.xticks(rotation=45)
    plt.title('Contribution to Total Value (%)')
    plt.ylabel('Percentage of Total Value')
    plt.grid(True, alpha=0.3)

plt.tight_layout()
plt.show()

Code Explanation

Setting up the problem:
- We defined $10$ items with corresponding values and weights
- The knapsack capacity is set to $40$ units
Formulating the integer program with PuLP:
- Created binary decision variables for each item (x)
- Set the objective function to maximize total value
- Added a constraint to ensure the total weight doesn’t exceed capacity
Solving the problem:
- $PuLP$ uses branch and bound algorithms to solve the integer program
- Results are extracted and organized into a pandas DataFrame
Visualization:
- Four plots are created to analyze different aspects of the solution:
  - All items plotted by weight vs. value
  - Selected vs. not selected items
  - Value/weight efficiency ratio for each item
  - Each selected item’s contribution to the total value

Solution Analysis

When executed, this code will select the optimal combination of items that maximizes value while staying within the weight limit.

The optimization model makes trade-offs based on both the absolute value of items and their weight efficiency (value per unit weight).

Items with high values and relatively low weights are preferred.

The visualizations help us understand:

The distribution of all items in the value-weight space
Which items were selected by the optimizer
The efficiency (value/weight) of each item
How much each selected item contributes to the total solution value

This example demonstrates how integer programming can be applied to combinatorial optimization problems, finding the optimal solution among a finite but potentially very large set of possibilities.

Result

Status: Optimal

Selected Items:
     Item  Value  Weight
0  Item_1     10       5
1  Item_4     31      16
2  Item_5      7       3
3  Item_8     22      11
4  Item_9     11       5

Total Value: 81
Total Weight: 40
Capacity: 40