企业官网建设流程全解析

各位大佬看看出了什么问题，为何模拟结果差那么多

我的代码，源代码没找到，不会用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*4000=1600万粒子防止内存溢出 # 星体配色调色板，最多支持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] ], dtype=np.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, endpoint=False) 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(parallel=True, cache=True, fastmath=True) 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, dtype=np.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, size=total_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, dtype=np.int16) # -1代表无归属 alive_full = np.ones(total_pix, dtype=np.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), dtype=np.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], dtype=np.float32) # 画布绘制 fig, ax = plt.subplots(figsize=(10, 10)) ax.imshow( rgb_img, extent=[-SPACE_HALF, SPACE_HALF, -SPACE_HALF, SPACE_HALF], origin="lower", interpolation="bilinear" ) # 绘制引力星体标记点，置于图层最上方 ax.scatter( STAR_POS[:, 0], STAR_POS[:, 1], c=STAR_COLORS, s=150, edgecolors="white", linewidth=1.2, zorder=10 ) ax.set_aspect("equal") ax.set_title( "正" + str(N_SIDE) + "边形引力吸引域 | 半隐式欧拉积分 | " + str(RESOLUTION) + "×" + str(RESOLUTION) + " | 速度上限Vmax=" + str(V_MAX) ) # 保存高清图片 plt.savefig("ngon_attractor_4000px.png", dpi=150, bbox_inches="tight") print("图片已保存至 ngon_attractor_4000px.png") plt.show()

大佬的模拟成果

我的结果