如何在R语言粒子碰撞模拟中高效实现多粒子扩展？

问题描述 #

我目前完成了单粒子模拟代码，原本计划为每个新粒子创建独立的数据框，再在ggplot中添加新的geom_point图层，但担心添加5个、10个甚至更多粒子时会导致代码过于冗长。请问是否存在无需编写大量重复代码即可为该模拟添加多粒子的方法？

原单粒子模拟代码：

library(tidyverse)
library(magick)
library(gganimate)
library(stringi)
library(ggtext)

particle_mass <- 4.65 #(x 10^-26) 
## A constant that dictates the number of decimal places that the data will be rounded to. Having this set as a constant avoids having to manually adjust a number which appears multiples times throughout the code if you wanted to change the number of decimal places for any reason.
rnd <- 4 
## initial particle coordinates 
x_i <- round(runif(n = 1, min = 0,max = 1), rnd) 
y_i <- round(runif(n = 1, min = 0,max = 1), rnd)

x_velocity <- 400 ## A constant for the x velocity in m/s
y_velocity <- 200 ##A  constant for the y velocity in m/s
time_step <- 0.001 ## A constant for the duration of a single time step in seconds
frame_total <- 200 ## A constant for the number of frames that the animation will be split into (two additional frames are added on at the end)
rate_factor <- 35 ## A factor that controls the playback rate of the animation. Increase this number and the animation will be slowed 

## The x and y displacement components per time step with a factor to slow down the animation to a more easily perceptible speed.
x_jump <- (x_velocity/rate_factor)/frame_total 
y_jump <- (y_velocity/rate_factor)/frame_total 

## Vectors of equal length for x position, y position, and corresponding time steps are initialized.
x_pos <- round(seq(x_i, (x_i + (frame_total*x_jump))-x_jump, x_jump), rnd)
y_pos <- round(seq(y_i, (y_i + (frame_total*y_jump))-y_jump, y_jump), rnd) 
times <- seq(0, (frame_total*time_step)-time_step, time_step)

## A loop that adjusts the direction of the x velocity vector when x = 1 or x = 0 to simulate wall contact
for (i in 2:length(x_pos)) {
  if (x_pos[i-1] > 1 | x_pos[i-1] < 0 ){x_jump = x_jump*(-1)}
  x_pos[[i]] <- round(x_pos[[i-1]] + x_jump, rnd)
}

##  A loop that adjusts the direction of the y velocity vector when y = 1 or y = 0 to simulate wall contact
for (i in 2:length(y_pos)) {
  if (y_pos[i-1] > 1 | y_pos[i-1] < 0 ){y_jump = y_jump*(-1)}
  y_pos[[i]] <- round(y_pos[[i-1]] + y_jump, rnd)
}


bounce <- data.frame(times, x_pos, y_pos) ##time and coordinate position information
bounce$x_dif <- c(0, diff(bounce$x_pos))## change in x position between time steps
bounce$x_turn <- c(0, diff(sign(bounce$x_dif))) ## helps detect changes in direction

## impulse magnitude experienced by the right wall at each time step
bounce <- bounce %>% mutate(impulse_right = case_when((bounce$x_turn) == -2 ~ 2*(x_velocity*particle_mass),
                                                      (bounce$x_turn) != -2 ~ 0))

bounce$cume_impulse_right <- cumsum(bounce$impulse_right) #cumulative impulse on the right wall 
## static plots of the moving particle
bounce_plot <- ggplot(bounce) +
  geom_point(aes(x_pos, y_pos), color = "azure4", size = 1) +
  labs(title = "Single gas particle model") + 
  ylab("") + xlab("") +
  theme_classic() +
  ## axis lines and tick marks
  theme(axis.line = element_blank(),
        axis.ticks.x = element_blank(),
        axis.text.x = element_blank(), axis.ticks.y = element_blank(),
        axis.text.y = element_blank()) +
  ## plot title aesthetics
  theme(plot.title = element_text(family = "Arial", face = "bold", hjust = 0.5, 
                                  size = 9, color = "gray38")) +
  ## The vertical and horizontal lines create a border around the 1 x 1 box that contains the particle
  geom_vline(xintercept=0, color="snow3", size=2) + 
  geom_hline(yintercept=0, color="snow3", size=2) + 
  geom_hline(yintercept=1, color="snow3", size=2) +
  geom_vline(xintercept=1, color="#CC6666", size=1) 

bounce_anim <- bounce_plot + transition_time(times)

## The desired number of frames and duration of the animation are set within the 'animate' function.
bounce_gif <- animate(bounce_anim, width = 3, height = 3, units = "in", res = 200,
                      nframes = frame_total, renderer = magick_renderer())

优化方案 #

核心思路是用函数封装单粒子轨迹生成逻辑，批量生成多粒子数据并统一管理，避免重复代码，同时用ggplot一次绘制所有粒子。

1. 封装单粒子轨迹生成函数 #

把单个粒子的运动模拟逻辑写成可复用函数，输入粒子ID、初始参数，输出包含完整轨迹的数据框：

# 生成单个粒子的轨迹数据
generate_particle_trajectory <- function(particle_id, x_i, y_i, x_velocity, y_velocity, 
                                         time_step, frame_total, rate_factor, rnd, particle_mass) {
  # 计算每帧位移
  x_jump <- (x_velocity/rate_factor)/frame_total 
  y_jump <- (y_velocity/rate_factor)/frame_total 
  
  # 初始化位置向量
  x_pos <- numeric(frame_total)
  y_pos <- numeric(frame_total)
  x_pos[1] <- x_i
  y_pos[1] <- y_i
  
  # 模拟x方向碰壁反弹
  current_x_jump <- x_jump
  for (i in 2:frame_total) {
    if (x_pos[i-1] > 1 | x_pos[i-1] < 0) {
      current_x_jump <- current_x_jump * (-1)
    }
    x_pos[i] <- round(x_pos[i-1] + current_x_jump, rnd)
  }
  
  # 模拟y方向碰壁反弹
  current_y_jump <- y_jump
  for (i in 2:frame_total) {
    if (y_pos[i-1] > 1 | y_pos[i-1] < 0) {
      current_y_jump <- current_y_jump * (-1)
    }
    y_pos[i] <- round(y_pos[i-1] + current_y_jump, rnd)
  }
  
  # 生成时间序列
  times <- seq(0, (frame_total*time_step)-time_step, time_step)
  
  # 整理数据并计算冲量（保留原逻辑）
  trajectory <- data.frame(
    particle_id = particle_id,
    times = times,
    x_pos = x_pos,
    y_pos = y_pos
  ) %>%
    group_by(particle_id) %>%
    mutate(
      x_dif = c(0, diff(x_pos)),
      x_turn = c(0, diff(sign(x_dif))),
      impulse_right = case_when(
        x_turn == -2 ~ 2*(x_velocity*particle_mass),
        TRUE ~ 0
      ),
      cume_impulse_right = cumsum(impulse_right)
    ) %>%
    ungroup()
  
  return(trajectory)
}

2. 批量生成多粒子数据 #

设定粒子数量，用map_df批量调用函数，生成所有粒子的合并数据框：

# 模拟参数（保留你的原始设置）
particle_mass <- 4.65
rnd <- 4 
time_step <- 0.001 
frame_total <- 200 
rate_factor <- 35 

# 设定需要模拟的粒子数量
num_particles <- 10

# 批量生成所有粒子的轨迹数据
all_particles <- map_df(1:num_particles, function(id) {
  # 随机生成每个粒子的初始位置和速度（可根据需求调整）
  x_i <- round(runif(1, 0, 1), rnd)
  y_i <- round(runif(1, 0, 1), rnd)
  x_velocity <- sample(c(-400, 400), 1) # 随机左右运动方向
  y_velocity <- sample(c(-200, 200), 1) # 随机上下运动方向
  
  generate_particle_trajectory(
    particle_id = id,
    x_i = x_i,
    y_i = y_i,
    x_velocity = x_velocity,
    y_velocity = y_velocity,
    time_step = time_step,
    frame_total = frame_total,
    rate_factor = rate_factor,
    rnd = rnd,
    particle_mass = particle_mass
  )
})

3. 绘制多粒子动画 #

只需一次调用geom_point，用粒子ID区分颜色，无需重复添加图层：

# 绘制多粒子动画
bounce_plot <- ggplot(all_particles) +
  geom_point(aes(x_pos, y_pos, color = factor(particle_id)), size = 1) +
  labs(title = "多气体粒子模型", color = "粒子ID") + 
  ylab("") + xlab("") +
  theme_classic() +
  theme(
    axis.line = element_blank(),
    axis.ticks = element_blank(),
    axis.text = element_blank(),
    plot.title = element_text(family = "Arial", face = "bold", hjust = 0.5, 
                              size = 9, color = "gray38")
  ) +
  geom_vline(xintercept=0, color="snow3", size=2) + 
  geom_hline(yintercept=0, color="snow3", size=2) + 
  geom_hline(yintercept=1, color="snow3", size=2) +
  geom_vline(xintercept=1, color="#CC6666", size=1) 

bounce_anim <- bounce_plot + transition_time(times)

bounce_gif <- animate(bounce_anim, width = 3, height = 3, units = "in", res = 200,
                      nframes = frame_total, renderer = magick_renderer())

方案优势 #

• 扩展性强：只需修改num_particles数值，即可快速增加粒子数量
• 代码简洁：用函数封装重复逻辑，避免冗余代码
• 可维护性高：调整运动规则时，仅需修改generate_particle_trajectory函数

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