各位大佬看看出了什么问题为何模拟结果差那么多我的代码源代码没找到不会用github原作者网址网页链接import numpy as np import matplotlib.pyplot as plt from numba import njit, prange # # 全局仿真参数区 # 适用场景N边形引力中心划分吸引域、半隐式欧拉粒子追踪、4000×4000高分辨率网格 # N_SIDE 6 # 引力中心构成正N边形边数 R_CIRCUM 1 # 正N边形外接圆半径 COLLIDE_RADIUS 0.05 # 引力星体捕获碰撞半径粒子进入即判定归属 SPACE_HALF 4.0 # 仿真空间半边长全域 [-4, 4] × [-4, 4] RESOLUTION 600 # 图像分辨率 4000×4000 网格粒子 V_MAX 20.0 # 【修改需求】速度上限替代原加速度限幅速度超过15会截断 V0 0.0 # 粒子初始速度全部置0 DAMP_ALPHA 0.012 # 速度阻尼系数模拟耗散摩擦加速粒子向星体收敛 DT 0.01 # 积分时间步长 MAX_STEP 10000 # 最大迭代步数粒子全部捕获可提前终止 G 1.0 # 万有引力常数 STAR_MASS 1.0 # 单个引力星体质量统一为1 SOFTEN 1e-3 # 引力软化系数避免星体附近引力无穷大奇点 CHUNK_SIZE 400000 # 分块粒子数量拆分4000*40001600万粒子防止内存溢出 # 星体配色调色板最多支持10个引力中心 COLOR_PALETTE np.array([ [0.92, 0.18, 0.18], [0.16, 0.74, 0.22], [0.18, 0.42, 0.90], [0.94, 0.82, 0.08], [0.84, 0.16, 0.82], [0.06, 0.88, 0.86], [0.96, 0.48, 0.06], [0.52, 0.30, 0.70], [0.26, 0.82, 0.58], [0.78, 0.22, 0.46] ], dtypenp.float32) # Matplotlib全局绘图配置解决中文、负号显示问题 plt.rcParams.update({ font.sans-serif: [SimHei, Microsoft YaHei], axes.unicode_minus: False, figure.dpi: 150 }) def generate_ngon_stars(n: int, r: float) - np.ndarray: 生成正N边形顶点坐标作为引力星体位置 :param n: 多边形边数 :param r: 外接圆半径 :return: shape(n,2) float32 星体xy坐标数组 theta np.linspace(0, 2 * np.pi, n, endpointFalse) xs r * np.cos(theta) ys r * np.sin(theta) return np.column_stack((xs, ys)).astype(np.float32) # 预计算星体坐标、数量、对应颜色 STAR_POS generate_ngon_stars(N_SIDE, R_CIRCUM) STAR_NUM STAR_POS.shape[0] STAR_COLORS COLOR_PALETTE[:STAR_NUM] # # Numba高性能核心半隐式欧拉积分求解粒子运动 # 并行分块计算nopython无Python对象模式禁用f-string兼容旧版numba # 核心更新移除加速度限幅增加速度上限V_MAX截断 # njit(parallelTrue, cacheTrue, fastmathTrue) def implicit_euler_particle_chunk( pos: np.ndarray, vel: np.ndarray, target_id: np.ndarray, alive_mask: np.ndarray, star_pos: np.ndarray, star_mass: float, g: float, damp_alpha: float, dt: float, collide_r: float, max_step: int, soften: float, v_max: float ): n_particles pos.shape[0] # 全局时间迭代循环 for step in range(max_step): # 无存活粒子直接提前退出节省计算 if not np.any(alive_mask): break # 提取当前存活粒子下标 alive_idx np.nonzero(alive_mask)[0] n_alive alive_idx.size pos_alive pos[alive_idx] vel_alive vel[alive_idx] acc_alive np.zeros_like(pos_alive) # 并行遍历每一个存活粒子计算引力加速度 for i in prange(n_alive): px, py pos_alive[i, 0], pos_alive[i, 1] ax, ay 0.0, 0.0 # 累加所有星体对当前粒子的引力 for s in range(STAR_NUM): sx, sy star_pos[s] dx sx - px dy sy - py r_sq dx * dx dy * dy soften r np.sqrt(r_sq) coeff g * star_mass / r_sq ax coeff * dx / r ay coeff * dy / r acc_alive[i, 0] ax acc_alive[i, 1] ay # --------------------------半隐式欧拉核心公式-------------------------- # 阻尼项隐式处理避免显式阻尼带来的数值发散 denom 1.0 damp_alpha * dt vel_new (vel_alive dt * acc_alive) / denom # 【需求修改】速度上限截断替代原加速度限幅 for i in range(n_alive): vx, vy vel_new[i, 0], vel_new[i, 1] v_mag np.sqrt(vx * vx vy * vy) if v_mag v_max: scale v_max / v_mag vel_new[i, 0] vx * scale vel_new[i, 1] vy * scale # 更新粒子位置 pos_new pos_alive dt * vel_new # 碰撞捕获判定粒子进入星体碰撞半径则标记死亡归属对应星体 for i in range(n_alive): px, py pos_new[i, 0], pos_new[i, 1] captured False for s in range(STAR_NUM): sx, sy star_pos[s] dist np.sqrt((sx - px) ** 2 (sy - py) ** 2) if dist collide_r: gid alive_idx[i] target_id[gid] s alive_mask[gid] False captured True break # 未被捕获则更新位置与速度 if not captured: pos_alive[i] pos_new[i] vel_alive[i] vel_new[i] # 写回分块数组 pos[alive_idx] pos_alive vel[alive_idx] vel_alive # 每1000步打印进度禁用f-string避免Numba报错 if step % 1000 0: remain np.sum(alive_mask) print(当前迭代步数:, step, 剩余未捕获粒子:, remain) return pos, vel, target_id # # 分块调度函数解决4000×4000网格一次性加载内存溢出问题 # 将全部粒子按CHUNK_SIZE切分逐块送入Numba核心计算结果回写全局数组 # def solve_full_domain(): # 生成均匀网格坐标 axis np.linspace(-SPACE_HALF, SPACE_HALF, RESOLUTION, dtypenp.float32) xx, yy np.meshgrid(axis, axis) total_pix RESOLUTION * RESOLUTION # 全局全量粒子数组 pos_full np.column_stack((xx.ravel(), yy.ravel())) vel_full np.zeros_like(pos_full) # rand_angle np.random.uniform(0, 2 * np.pi, sizetotal_pix).astype(np.float32) # vel_full[:, 0] V0 * np.cos(rand_angle) # vel_full[:, 1] V0 * np.sin(rand_angle) target_id_full np.full(total_pix, -1, dtypenp.int16) # -1代表无归属 alive_full np.ones(total_pix, dtypenp.bool_) # 分块循环处理 chunk_idx 0 for start in range(0, total_pix, CHUNK_SIZE): chunk_idx 1 end min(start CHUNK_SIZE, total_pix) print(\n开始处理第, chunk_idx, 块粒子区间 [, start, :, end, ]) # 切片拷贝当前分块数据 pos_chunk pos_full[start:end].copy() vel_chunk vel_full[start:end].copy() tid_chunk target_id_full[start:end].copy() alive_chunk alive_full[start:end].copy() # 调用Numba半隐式求解器传入速度上限V_MAX pos_chunk, vel_chunk, tid_chunk implicit_euler_particle_chunk( pos_chunk, vel_chunk, tid_chunk, alive_chunk, STAR_POS, STAR_MASS, G, DAMP_ALPHA, DT, COLLIDE_RADIUS, MAX_STEP, SOFTEN, V_MAX ) # 计算结果回写全局数组 pos_full[start:end] pos_chunk vel_full[start:end] vel_chunk target_id_full[start:end] tid_chunk alive_full[start:end] alive_chunk # 将一维归属标签重塑为二维图像矩阵 result_img target_id_full.reshape(RESOLUTION, RESOLUTION) return result_img # # 程序入口仿真计算 结果绘图输出 # if __name__ __main__: total_particles RESOLUTION * RESOLUTION print(仿真任务启动) print(分辨率, RESOLUTION, ×, RESOLUTION) print(总粒子数量, total_particles) print(引力中心正, N_SIDE, 边形布局) print(速度上限限制, V_MAX) # 执行分块粒子追踪仿真 result_label solve_full_domain() print(\n仿真计算完成开始渲染图像) # 构建RGB彩色图像 rgb_img np.zeros((RESOLUTION, RESOLUTION, 3), dtypenp.float32) # 按星体ID填充对应颜色 for sid in range(STAR_NUM): mask result_label sid rgb_img[mask] STAR_COLORS[sid] # 未被捕获粒子填充灰色背景 rgb_img[result_label -1] np.array([0.3, 0.3, 0.3], dtypenp.float32) # 画布绘制 fig, ax plt.subplots(figsize(10, 10)) ax.imshow( rgb_img, extent[-SPACE_HALF, SPACE_HALF, -SPACE_HALF, SPACE_HALF], originlower, interpolationbilinear ) # 绘制引力星体标记点置于图层最上方 ax.scatter( STAR_POS[:, 0], STAR_POS[:, 1], cSTAR_COLORS, s150, edgecolorswhite, linewidth1.2, zorder10 ) ax.set_aspect(equal) ax.set_title( 正 str(N_SIDE) 边形引力吸引域 | 半隐式欧拉积分 | str(RESOLUTION) × str(RESOLUTION) | 速度上限Vmax str(V_MAX) ) # 保存高清图片 plt.savefig(ngon_attractor_4000px.png, dpi150, bbox_inchestight) print(图片已保存至 ngon_attractor_4000px.png) plt.show()大佬的模拟成果我的结果