Как визуализировать коллективное поведение самоходки в двух измерениях? ⇐ Python

[Anonymous]

несколько номеров. Самодвижущиеся стержни (смоделированные с использованием нечетного числа соединенных твердых сфер) с фиксированной самодвижущейся скоростью движутся в среде (2D) с тремя различными константами диффузии для вращательной и поступательной диффузии (x и y), каждый стержень имеет разную ориентацию. когда две сферы двух разных стержней оказываются в пределах определенного радиуса взаимодействия, между ними действует сила взаимодействия, которая возникает из-за потенциала Леннарда-Джонса, который является функцией расстояния между двумя взаимодействующими сферами. составляющая силы взаимодействия складывается со скоростью перемещения, и возникающий крутящий момент медленно вращает оба стержня, и они начинают выравниваться.
уравнение движения стержней и форма потенциала приложены.введите описание изображения здесь введите описание изображения здесь
это код, который я написал,

Код: Выделить всё

import numpy as np
import matplotlib.pyplot as plt
import matplotlib.animation as animation
import matplotlib.patches as patches

# Constants and parameters
WIDTH, HEIGHT = 100, 100
NUM_SPHERES_PER_ROD = 9  # Odd number of spheres per rod
SPHERE_RADIUS = 0.5
NUM_RODS = 10
VELOCITY = 2
FRAMES = 10
NEIGHBOR_DISTANCE_THRESHOLD = 2 * SPHERE_RADIUS  # Diameter of a sphere
EPSILON = 0.000001  # Depth of the Lennard-Jones potential well
SIGMA = 2 * SPHERE_RADIUS  # Distance at which the potential is zero

def initialize_rods(num_rods):
positions = np.random.rand(num_rods, 2) * [WIDTH, HEIGHT]
angles = np.random.uniform(0, 2 * np.pi, num_rods)
return positions, angles

def update_positions(positions, angles, velocity):
theta = angles
dx = velocity * np.cos(theta)
dy = velocity * np.sin(theta)
positions[:, 0] += dx
positions[:, 1] += dy

# Periodic Boundary
positions[:, 0] = positions[:, 0] % WIDTH
positions[:, 1] = positions[:, 1] % HEIGHT
return positions

def get_sphere_positions(start_pos, angle, num_spheres):
sphere_positions = []
for i in range(int((num_spheres + 1) / 2)):
x = start_pos[0] - 2 * i * SPHERE_RADIUS * np.cos(angle)
y = start_pos[1] - 2 * i * SPHERE_RADIUS * np.sin(angle)
x1 = start_pos[0] + 2 * i * SPHERE_RADIUS * np.cos(angle)
y1 = start_pos[1] + 2 * i * SPHERE_RADIUS * np.sin(angle)
sphere_positions.append((x, y))
sphere_positions.append((x1, y1))
return sphere_positions

def lennard_jones_force(r):
if r == 0:
return 0  # Avoid division by zero
r6 = (SIGMA / r) ** 6
r12 = r6 * r6
force_magnitude = 24 * EPSILON * (2 * r12 - r6) / r
return force_magnitude

def compute_forces_and_neighbors(positions, angles):
num_rods = len(positions)
neighbor_sphere_pairs = []
forces = np.zeros_like(positions)  # Forces to be applied to each rod

for i in range(num_rods):
sphere_positions_i = get_sphere_positions(positions[i], angles[i], NUM_SPHERES_PER_ROD)
for j in range(i + 1, num_rods):
sphere_positions_j = get_sphere_positions(positions[j], angles[j], NUM_SPHERES_PER_ROD)
for sphere_pos_i in sphere_positions_i:
for sphere_pos_j in sphere_positions_j:
dist_vector = np.array(sphere_pos_j) - np.array(sphere_pos_i)
dist = np.linalg.norm(dist_vector)
if dist < NEIGHBOR_DISTANCE_THRESHOLD:
# Record the neighbor sphere pair
neighbor_sphere_pairs.append((sphere_pos_i, sphere_pos_j))
# Calculate the Lennard-Jones force
force_mag = lennard_jones_force(dist)
force_vec = force_mag * (dist_vector / dist)  # Normalize to get direction
# Apply forces to the rods
forces[i] += force_vec
forces[j] -= force_vec  # Equal and opposite force

return forces, neighbor_sphere_pairs

def update_angles(angles, forces):
for i in range(len(angles)):
# Compute resultant force
resultant_force = forces[i]
if np.linalg.norm(resultant_force) > 0:
# Adjust angle based on resultant direction
angles[i] = np.arctan2(resultant_force[1], resultant_force[0])
return angles

def draw_sphere(ax, center, radius, COLOUR):
circle = patches.Circle(center, radius, color=COLOUR, fill=True)
ax.add_patch(circle)

def draw_rod(ax, start_pos, angle, num_spheres):
sphere_positions = get_sphere_positions(start_pos, angle, num_spheres)
for idx, (x, y) in enumerate(sphere_positions):
COLOUR = "blue"  if idx % 2 == 0 else "red"
draw_sphere(ax, (x, y), SPHERE_RADIUS, COLOUR)

def animate(frame):
global positions, angles
positions = update_positions(positions, angles, VELOCITY)
forces, neighbor_sphere_pairs = compute_forces_and_neighbors(positions, angles)  # Find neighbors and compute forces
angles = update_angles(angles, forces)  # Update rod angles based on forces

ax.clear()
for idx, (pos, angle) in enumerate(zip(positions, angles)):
draw_rod(ax, pos, angle, NUM_SPHERES_PER_ROD)

ax.set_xlim(0, WIDTH)
ax.set_ylim(0, HEIGHT)
ax.set_aspect('equal')

# Initialize
positions, angles = initialize_rods(NUM_RODS)

# Set up the figure and axis
fig, ax = plt.subplots(figsize=(10, 8))

# Create the animation
anim = animation.FuncAnimation(fig, animate, frames=FRAMES, interval=10, repeat=True)
plt.show()

но вращение стержней после взаимодействия не является точным. можете ли вы что-нибудь предложить?

Подробнее здесь: https://stackoverflow.com/questions/790 ... -dimension