Numerical Methods in Physics

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

def monte_carlo_integrate(func, dim, n_samples=100000):
    """
    Monte Carlo integration in arbitrary dimensions
    ∫ f(x) dx over [0,1]^dim
    """
    # Generate random samples
    samples = np.random.random((n_samples, dim))
    
    # Evaluate function at sample points
    values = func(samples)
    
    # Calculate Monte Carlo estimate
    volume = 1.0  # Volume of [0,1]^dim
    integral_estimate = volume * np.mean(values)
    
    # Estimate error (standard deviation of mean)
    error_estimate = volume * np.std(values) / np.sqrt(n_samples)
    
    return integral_estimate, error_estimate

def example_function_6d(x):
    """
    6-dimensional example function: 
    f(x) = exp(-∑ (x_i - 0.5)^2)
    """
    r_squared = np.sum((x - 0.5)**2, axis=1)
    return np.exp(-10 * r_squared)

# Compare Monte Carlo with other methods
def compare_integration_methods():
    dim = 6
    n_samples = 100000
    
    print(f"Integrating 6-dimensional function")
    print(f"Monte Carlo samples: {n_samples:,}")
    print("-" * 50)
    
    # Monte Carlo integration
    start_time = time.time()
    mc_integral, mc_error = monte_carlo_integrate(
        example_function_6d, dim, n_samples
    )
    mc_time = time.time() - start_time
    
    # Quadrature (only works for low dimensions)
    if dim <= 3:
        start_time = time.time()
        quad_result = integrate.nquad(
            lambda *args: example_function_6d(np.array([args]).T),
            [[0,1]]*dim
        )[0]
        quad_time = time.time() - start_time
    
    print(f"Monte Carlo Result: {mc_integral:.6f} ± {mc_error:.6f}")
    print(f"Monte Carlo Time: {mc_time:.3f} seconds")
    print(f"Monte Carlo Efficiency: {n_samples/mc_time:.0f} samples/sec")
    
    if dim <= 3:
        print(f"\nQuadrature Result: {quad_result:.6f}")
        print(f"Quadrature Time: {quad_time:.3f} seconds")
        print(f"Speedup: {quad_time/mc_time:.1f}x faster")
    
    # Error convergence analysis
    print("\n" + "-" * 50)
    print("Monte Carlo Error Convergence (1/√N):")
    
    sample_sizes = [1000, 5000, 20000, 50000, 100000, 500000]
    errors = []
    
    for n in sample_sizes:
        _, error = monte_carlo_integrate(example_function_6d, dim, n)
        errors.append(error)
    
    # Plot convergence
    plt.figure(figsize=(10, 6))
    plt.loglog(sample_sizes, errors, 'bo-', linewidth=2, markersize=8)
    plt.loglog(sample_sizes, 1/np.sqrt(sample_sizes), 'r--', 
               linewidth=2, label='1/√N reference')
    plt.xlabel('Number of Samples (N)', fontsize=12)
    plt.ylabel('Error Estimate', fontsize=12)
    plt.title('Monte Carlo Error Convergence (6D Integration)', fontsize=14)
    plt.grid(True, alpha=0.3)
    plt.legend()
    plt.tight_layout()
    plt.show()

# Run comparison
if __name__ == "__main__":
    compare_integration_methods()