Python中测试粒子模拟计算长宽比的问题

Question

我不断收到此错误，但无法解决此问题。我想要每个时间步的长度和宽度（所有粒子的平均值）。因此，长度和宽度数组应该具有 (19,) 形状，而不是当前的 (1000,) 形状，这意味着它正在计算每个粒子。

错误：

‘plt.plot(t_输出，比率)

ValueError：x 和 y 必须具有相同的第一维度，但具有形状 (19,) 和 (1000,)’

我已经尽了我最好的“编码”知识，这对于经验丰富的专家编码员来说应该非常简单，但遗憾的是，我不是其中之一。在脚本的最后几行中，我从“#提取所有粒子和所有时间步长的速度”开始评论了我应该做什么以及我尝试了什么

请帮助我修复此错误并实现我想要做的事情！

import numpy as np
import matplotlib.pyplot as plt
from scipy.spatial.transform import Rotation
from mpl_toolkits.mplot3d import Axes3D


# Constants
M_Sun = 1.989e33  # Solar Mass in g
M_bh = 4.3e6 *M_Sun # Mass of BH
G = 6.67430e-8  # cm^3 g^(-1) s^(-2)
yr = 365 * 24 * 60 * 60  # 1 year in seconds
pc = 3.086e18 # Parsec in cm
R = .225 * 0.04 * pc / 2 # Radius of cloud

# Number of particles
num_particles = 1000

# Uniformly distributed particles in Sphere of radius R
phi = np.random.uniform(0, 2*np.pi, size=num_particles)
costheta = np.random.uniform(0, 1, size=num_particles)
u = np.random.uniform(0, 1, size=num_particles)

theta = np.arccos(costheta)
r = R * (u**(1/3))

x = r * np.sin(theta) * np.cos(phi)
y = r * np.sin(theta) * np.sin(phi)
z = r * np.cos(theta)

# Initial Conditions

RA = .425 # Right Ascension
Decl = .523 # Declination

vx = 550e5 # x-component of velocity
vy = 730e5 # y-component of velocity
vz = 0 # z-component of velocity

# Transform RA and Decl to x,y,z components (1 arcsec = 0.04 pc in the galactic center)
x_0 = RA * 0.04 * pc
y_0 = Decl * 0.04 * pc
z_0 = 0

# Empty lists for position and velocity
initial_pos = []
initial_vel = []

for i in range(num_particles):
    initial_position = (x_0 + x[i], y_0 + y[i], z_0 + z[i]) # x_0 + x (x is already multiplied by R), same for y and z axes
    initial_velocity = (vx, vy, vz) # Fixed velocity for all particles
    initial_pos.append(initial_position)
    initial_vel.append(initial_velocity)

initial_pos = np.array(initial_pos) #Shape (num_particles, coordinates)
initial_vel = np.array(initial_vel) #Shape (num_particles, coordinates)

# Time
t_end = 100*yr # Total time of integration
dt_constant = 0.1 # Magic number for time stepping scheme
intervals = 1000000 # seconds after which to store pos & vel
next_interval_time = intervals

# Initial conditions
time = np.zeros(1)
_pos_t = initial_pos.copy()
_vel_t = initial_vel.copy()

# Lists to store outputs
pos_output = []
vel_output = []
t_output = []

# Number of Stored outputs
output = 20

while time[-1] <= t_end: #We end up doing one more timestep after t_end
    r = np.linalg.norm(_pos_t, axis=1)
    acc = -G * M_bh / r[:, np.newaxis]**3 * _pos_t #np.newaxis for broadcasting with _pos_t

    # Calculate the time step for the current particle
    current_dt = dt_constant * np.sqrt(np.linalg.norm(_pos_t, axis=1)**3 / (G * M_bh))
    min_dt = np.min(current_dt) # Use the minimum time step for all particles

    # leap frog integration

    # calculate velocity at half timestep
    half_vel = _vel_t + 0.5 * acc * min_dt
    # current position only used in this timestep
    _pos_t = _pos_t + half_vel * min_dt # shape: (#particles, #coordinates)

    # Recalculate acceleration with the new position
    r = np.linalg.norm(_pos_t, axis=1) 
    # Acceleration at timestep
    acc = -G * M_bh / r[:, np.newaxis]**3 * _pos_t #np.newaxis for broadcasting with _pos_t
    # current velocity only used in this timestep
    _vel_t = half_vel + 0.5 * acc * min_dt # shape: (#particles, #coordinates)

    # time at timestep t
    _time_t = time[-1] + min_dt

    # Check if the elapsed time has surpassed the next interval
    if _time_t >= next_interval_time: # >= because time steps aren't constant and it may not be an exact multiple
        # Store data at intervals
        pos_output.append(_pos_t.copy())
        vel_output.append(_vel_t.copy())
        t_output.append(_time_t.copy())
        next_interval_time += intervals

    time = np.append(time, _time_t)

    # show current status by printing timestep number (-1 because initial conditions)
    print(f'timestep: {time.size - 1} [progress: {_time_t/t_end*100:.3f}%]')


pos_output = np.array(pos_output).transpose((1, 0, 2))  # shape: (#objects, #stored timesteps, #coordinates)
vel_output = np.array(vel_output).transpose((1, 0, 2))  # shape: (#objects, #stored timesteps, #coordinates)
t_output = np.array(t_output)

# Plotting the ratio of length to breadth as a function of time for the entire cloud 
# find the velocity vector at each timestep 
# move to COM frame
# find the angle of x or y position to the vel vector #dot product
# rotate one axis so it aligns with vel vector # use rotation matrix
# then x and y are perpendicular and one is aligned with vel vector
# then do the ratio calculation

# Extract velocities for all particles and all timesteps
velocities = vel_output[:, :, :]

# Move to the Center of Mass (COM) frame
com_velocity = np.mean(velocities, axis=0)  # Calculate the COM velocity for all timesteps
velocities_rel_com = velocities - com_velocity  # Relative velocities to the COM frame

# Calculate the angle between x axis and the velocity vector
dot_products = np.sum(pos_output[:, :, 0] * velocities_rel_com[:, :, 0], axis=-1) # dot product for each particle and each timestep (x and velocity vector)

# Find the angles between the position and velocity vectors
magnitude_pos = np.linalg.norm(pos_output[:, :, 0]) # magnitude of the position vector
magnitudes_vel = np.linalg.norm(velocities_rel_com[:, :, 0]) # magnitude of the velocity vector
angles = np.arccos(dot_products / (magnitude_pos * magnitudes_vel)) # angle between the position and velocity vectors

# Rotate the axes so that one axis aligns with the velocity vector
rotation = Rotation.from_euler('y', angles)
rotation_matrix = rotation.as_matrix()

positions_ = pos_output.copy()
# Rotate the position vectors
pos_rotated = np.einsum('ijk,ikl->ijl', positions_, rotation_matrix)

# Calculate the length and width for each timestep and all particles
length = np.max(pos_rotated[:, :, 0], axis=-1) - np.min(pos_rotated[:, :, 0], axis=-1)
width = np.max(pos_rotated[:, :, 1], axis=-1) - np.min(pos_rotated[:, :, 1], axis=-1)

# Calculate the ratio for each timestep
ratio = length / width

# Plot the ratio as a function of time
plt.plot(t_output, ratio)
plt.xlabel('Time (s)')
plt.ylabel('Ratio')
plt.title('Ratio of Length to Width of Cloud')
plt.show()

Answer 1

错误发生在

plt.plot(t_output, ratio)

。向量的大小需要相同，但错误是

t_output

的长度为 19，而

ratio

的长度为 1000。我做了以下修改

length = np.max(pos_rotated[:, :, 0], axis=0) - np.min(pos_rotated[:, :, 0], axis=0)
width = np.max(pos_rotated[:, :, 1], axis=0) - np.min(pos_rotated[:, :, 1], axis=0)

之前是

axis=-1

。我把它改为

axis=0

。假设您采用

pos_rotated[:, :, 0]

- 该数组的形状为 (1000, 19) - not (1000, 19, 3)，这就是

pos_rotated

。鉴于其形状为 (1000, 19)，如果我们想要在 19 个时间步长中的每个时间步中找到所有 1000 个粒子的最大值，我们需要指定

axis=0

，即在粒子维度上执行

np.max

。您最终会得到一个形状为 (19) 的数组，对于每个时间点，它是超过 1000 个粒子的最大值。