在 Python 上计算并绘制 txt 文件中粒子轨迹的均方位移

Question

我想计算并绘制粒子轨迹的 MSD，以确定轨迹是否是由于介质中的纯扩散造成的，或者它们是否受到约束或受到外力的影响。我尝试使用 Python 进行此练习，但我总是遇到同样的问题：即使我的数据集不同，但我的图形仍然相同，或者我的 MSD 是一条水平线。

这是我试图解决问题的代码：

import numpy as np
import matplotlib.pyplot as plt

# Loading data from the txt file
data = np.loadtxt('/Users/arthurcoet/Documents/PhD COËT/CiNAM/Python/Exercices/Dtxt1.txt')

# Retrieving data columns
temps = data[:, 0]
indices = data[:, 1]
positions_x = data[:, 2]
positions_y = data[:, 3]

# Calculation of the MSD for each particle
particle_indices = np.unique(indices)
msd_particles = []

for particle_index in particle_indices:
    indices_particle = np.where(indices == particle_index)[0]
    x_particle = positions_x[indices_particle]
    y_particle = positions_y[indices_particle]
    displacements = np.sqrt((x_particle - x_particle[0])**2 + (y_particle - y_particle[0])**2)
    msd_particle = np.mean(displacements ** 2)
    msd_particles.append(msd_particle)

# Calculation of average MSD
msd_mean = np.mean(msd_particles)

# Plot of average MSD against time
plt.plot(temps, msd_mean * np.ones_like(temps), '-')
plt.xlabel('Temps')
plt.ylabel('MSD')
plt.title('MSD moyen en fonction du temps')
plt.show()

我还会向您展示我的txt文件的摘录，格式如下：时间，粒子索引，X，Y：

2.000000000000000042e-02 1.000000000000000000e+00 1.240807597782861649e+00 -3.785534313700077980e-01 2.000000000000000042e-02 2.000000000000000000e+00 -7.309480801080029400e-01 -7.170078639286830979e-01 2.000000000000000042e-02 3.000000000000000000e+00 8.098986703909789586e-01 5.815296326847829711e-01 2.000000000000000042e-02 4.000000000000000000e+00 2.750517661172900349e-03 5.522082763915119319e-01 2.000000000000000042e-02 5.000000000000000000e+00 9.127289184940016176e-01 5.289936832099206843e-01 2.000000000000000042e-02 6.000000000000000000e+00 1.340115717383320693e+00 -2.206167086146532119e-01 2.000000000000000042e-02 7.000000000000000000e+00 -9.347632084034850353e-01 -1.031270972401856945e+00 2.000000000000000042e-02 8.000000000000000000e+00 -1.564283425436366892e+00 9.906037061246882880e-01 2.000000000000000042e-02 9.000000000000000000e+00 -1.092736507472809315e-01 -7.019693105410058642e-01 2.000000000000000042e-02 1.000000000000000000e+01 1.079635002230163246e-02 1.936345833496776025e+00

以及粒子轨迹图的图像：

您是否也遇到过这个问题？如果有，您知道如何解决吗？

Answer 1

您可以尝试正确绘制 MSD 值与时间的关系图。

...
# Create a time array for plotting
unique_times = np.unique(temps)

# Initialize an array to store the average MSD values at each time point
msd_mean = np.zeros(len(unique_times))

# Calculate the average MSD for each time point
for i, time in enumerate(unique_times):
    time_indices = np.where(temps == time)[0]
    msd_mean[i] = np.mean([msd_particles[j] for j in time_indices])

# Plot average MSD against time
plt.plot(unique_times, msd_mean, '-')
...