基于MCMC的一维铁磁伊辛模型Python实现技术求助

Hey David, let's get your 1D Ising model MCMC code working correctly! Your draft had a few key missteps (like accidentally using a 2D lattice and the wrong Metropolis proposal mechanism), so let's break this down step by step.

Key Issues in Your Original Code #

• You set up a 2D lattice (Shape=(20,20)) but the problem specifies a 1D system
• The energy calculation was written for 2D interactions, not the 1D Hamiltonian provided
• The Metropolis-Hastings loop was using continuous values instead of flipping discrete ±1 spins (the core of the Ising model)

Corrected Full Code #

import numpy as np
import matplotlib.pyplot as plt

# --------------------------
# Parameter Setup (as requested)
# --------------------------
L = 20  # Lattice size (1D)
T = 2   # Temperature
Beta = 1 / T  # β = 1/T
B = 0   # External magnetic field

# Initialize spin configuration: random ±1 for each site in 1D
spins = np.random.choice([-1, 1], size=L)

# --------------------------
# Helper Functions
# --------------------------
def calculate_energy(spins, B):
    """Compute the Hamiltonian for the 1D Ising model."""
    # Nearest-neighbor interaction term: sum(σ_i * σ_{i+1}) for i=1 to L-1
    neighbor_interaction = np.sum(spins[:-1] * spins[1:])
    # External field term: sum(σ_i)
    field_term = np.sum(spins)
    # Hamiltonian: H = -J*sum(σ_iσ_{i+1}) - B*sum(σ_i) (J=1 here, ferromagnetic)
    return -neighbor_interaction - B * field_term

def calculate_magnetization(spins):
    """Compute the average magnetization per site."""
    return np.sum(spins) / L

# --------------------------
# Metropolis-Hastings MCMC Loop
# --------------------------
n_steps = 10000  # Number of MCMC steps
energy_history = []
magnetization_history = []

# Initial energy and magnetization
current_energy = calculate_energy(spins, B)
current_magnetization = calculate_magnetization(spins)
energy_history.append(current_energy)
magnetization_history.append(current_magnetization)

for step in range(n_steps):
    # 1. Propose a new state: flip a single random spin
    site_idx = np.random.randint(0, L)
    spins[site_idx] *= -1  # Flip the spin
    
    # 2. Calculate new energy and energy difference
    new_energy = calculate_energy(spins, B)
    delta_E = new_energy - current_energy
    
    # 3. Metropolis acceptance rule
    if delta_E <= 0:
        # Accept the flip: update current state
        current_energy = new_energy
    else:
        # Calculate acceptance probability
        accept_prob = np.exp(-Beta * delta_E)
        # Randomly accept/reject
        if np.random.rand() < accept_prob:
            current_energy = new_energy
        else:
            # Reject the flip: revert the spin
            spins[site_idx] *= -1
    
    # 4. Update magnetization and record values
    current_magnetization = calculate_magnetization(spins)
    energy_history.append(current_energy)
    magnetization_history.append(current_magnetization)

# --------------------------
# Post-Processing & Plotting
# --------------------------
plt.figure(figsize=(12, 5))

# Plot 1: Magnetization Histogram (skip burn-in period)
plt.subplot(1, 2, 1)
plt.hist(magnetization_history[1000:], bins=20, edgecolor='black')
plt.xlabel('Magnetization per Site')
plt.ylabel('Frequency')
plt.title('Magnetization Histogram (T=2, B=0)')

# Plot 2: Energy vs. MCMC Steps
plt.subplot(1, 2, 2)
plt.plot(energy_history)
plt.xlabel('MCMC Step')
plt.ylabel('Total Energy')
plt.title('Energy Evolution Over Time')

plt.tight_layout()
plt.show()

Key Improvements & Explanations #

1. 1D Lattice Setup: We use a 1D numpy array for spins instead of 2D, matching your problem's requirements.
2. Exact Hamiltonian Implementation: The calculate_energy function directly follows the formula you provided, with nearest-neighbor interactions and external field terms.
3. Correct Metropolis Proposal: Instead of continuous values, we flip a single random spin—this is the standard, efficient proposal mechanism for Ising models.
4. Burn-In Period: We skip the first 1000 steps in the magnetization histogram to let the chain reach equilibrium before collecting meaningful statistics.
5. Optimization Tip: For larger lattices, you can avoid recalculating the entire energy each time. Instead, compute only the energy change from flipping a single spin:

   # Handle open boundary conditions (no wrap-around)
   left = spins[site_idx-1] if site_idx > 0 else 0
   right = spins[site_idx+1] if site_idx < L-1 else 0
   delta_E = 2 * spins[site_idx] * (left + right + B)
   

   This cuts down computation time drastically for large systems.

Expected Results #

With T=2 (above the 1D Ising model's critical temperature T_c≈2.269), you'll see a magnetization histogram centered near 0 (paramagnetic phase) and the energy will stabilize quickly after the burn-in period.

内容的提问来源于stack exchange，提问作者David