基于欧拉方法的粒子扩散数值方程启停问题Python求解

Got it, let's tackle this 2D diffusion problem with absorbing boundaries using the explicit Euler method in Python. I'll break this down step by step so you can follow along easily.

Problem Breakdown #

First, let's align on the core details of the problem:

• We're simulating particle diffusion on a 24x24 grid representing a 2x2 mm surface. Each grid cell has a spatial step of dx = dy = 2/24 = 1/12 ≈ 0.0833 mm.
• Diffusion constant D = 0.025 (we'll assume units of mm²/s, standard for such problems).
• Absorbing boundaries: Any particle reaching the grid edge is immediately removed, so boundary cells always hold 0 particles.
• Initial condition: The 4 center grid cells start with a total of N moles of particles (we'll split this equally across the 4 cells for realism, but you can adjust this if needed).
• We need to run simulations for N = 10, 100, 10000 and stop when nearly all particles are absorbed.

Implementation Approach #

For the explicit Euler (finite difference) method, we must respect a stability condition to avoid numerical divergence. For 2D diffusion, this rule is:

dt ≤ (dx²) / (4*D)

We'll calculate dt to stay well within this limit to keep the simulation stable.

The update rule for each internal cell (non-boundary) comes from discretizing the 2D diffusion equation:

n_new[i,j] = n_old[i,j] + D*dt*( (n_old[i+1,j] - 2*n_old[i,j] + n_old[i-1,j])/dx² + (n_old[i,j+1] - 2*n_old[i,j] + n_old[i,j-1])/dy² )

Boundary cells remain 0 throughout the simulation to enforce the absorbing condition.

Full Python Code #

Here's a complete, configurable implementation with visualization:

import numpy as np
import matplotlib.pyplot as plt

def simulate_diffusion(N, grid_size=24, surface_size=2.0, D=0.025, termination_threshold=1e-6):
    # Grid and time step setup
    dx = surface_size / grid_size
    dy = dx
    # Use a factor of 5 to stay safely within the stability limit
    dt = (dx**2) / (5*D)
    
    # Initialize grid: all zeros, split N equally across 4 center cells
    grid = np.zeros((grid_size, grid_size))
    center = grid_size // 2
    grid[center-1:center+1, center-1:center+1] = N / 4
    
    # Track total particles over time to monitor absorption
    total_particles = [np.sum(grid)]
    
    # Simulation loop
    while total_particles[-1] > termination_threshold:
        new_grid = np.copy(grid)
        
        # Update only internal cells (skip boundaries)
        for i in range(1, grid_size-1):
            for j in range(1, grid_size-1):
                # Calculate Laplacian terms for x and y directions
                laplacian_x = (grid[i+1,j] - 2*grid[i,j] + grid[i-1,j]) / dx**2
                laplacian_y = (grid[i,j+1] - 2*grid[i,j] + grid[i,j-1]) / dy**2
                new_grid[i,j] = grid[i,j] + D * dt * (laplacian_x + laplacian_y)
        
        # Enforce absorbing boundaries (reset all edge cells to 0)
        new_grid[0, :] = 0
        new_grid[-1, :] = 0
        new_grid[:, 0] = 0
        new_grid[:, -1] = 0
        
        grid = new_grid
        total_particles.append(np.sum(grid))
    
    print(f"Simulation for N={N} finished after {len(total_particles)-1} time steps.")
    print(f"Final remaining particles: {total_particles[-1]:.6f}\n")
    
    return grid, total_particles

# Run simulations for all required N values
N_values = [10, 100, 10000]
simulation_results = {}
for N in N_values:
    final_grid, particle_history = simulate_diffusion(N)
    simulation_results[N] = (final_grid, particle_history)

# Plot particle loss over time for each N
plt.figure(figsize=(10,6))
for N in N_values:
    plt.plot(simulation_results[N][1], label=f'Initial N={N}')
plt.xlabel('Time Step')
plt.ylabel('Total Particles (mol)')
plt.title('Particle Absorption Over Time')
plt.legend()
plt.grid(True)
plt.show()

# Plot final particle distribution for N=10000 (most detailed pattern)
plt.figure(figsize=(8,8))
plt.imshow(simulation_results[10000][0], cmap='viridis', origin='lower')
plt.colorbar(label='Particles per Cell (mol)')
plt.title('Final Particle Distribution (N=10000, Absorbing Boundaries)')
plt.xlabel('Grid X Index')
plt.ylabel('Grid Y Index')
plt.show()

Key Notes & Customizations #

• Stability: The time step dt uses a safety factor of 5 to avoid divergence. You can adjust this to 4.5 for slightly faster simulations, but never exceed 4 (the theoretical limit).
• Termination Condition: The simulation stops when total particles drop below 1e-6 (nearly all absorbed). You can tweak this threshold (e.g., stop when 99% of particles are gone by checking total_particles[-1] < 0.01*N).
• Initial Condition: If you want each center cell to hold N mol instead of splitting the total, change the initialization line to grid[center-1:center+1, center-1:center+1] = N—just note this makes the total initial particles 4*N.
• Animation: To visualize the diffusion process in motion, you can add a FuncAnimation block inside the simulation function (let me know if you want help with that!).

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