绘制周期轨迹

Question

我有一些数据，描述了一个粒子在一个有边界的走廊里移动的情况。把这个粒子的轨迹画出来后，发现它的轨迹是锯齿形的。

我想知道怎么让plot()这个函数在粒子回到起点时，不把这些点连起来。就像图片上面那样，但不想要"."这个点。

我最初的想法是找出numpy数组a[:-1]-a[1:]中，哪里变成正数的索引，然后从0画到那个索引。但是我该怎么找到a[:-1]-a[1:]中第一个正数的索引呢？也许还有其他的想法。

Answer 1

根据Thorsten Kranz的回答，这个版本在'y'值穿过某个区间时，会给原始数据添加一些点。这一点很重要，特别是当数据点的密度不高时，比如说np.linspace(0., 100., 100)和原始的np.linspace(0., 100., 1000)。在这个过程中，曲线的变化位置是通过线性插值来计算的。把这些内容封装成一个函数，代码如下：

import numpy as np
def periodic2plot(x, y, period=np.pi*2.):
    indexes = np.argwhere(np.abs(np.diff(y))>.5*period).flatten()
    index_shift = 0
    for i in indexes:
        i += index_shift
        index_shift += 3  # in every loop it adds 3 elements

        if y[i] > .5*period:
            x_transit = np.interp(period, np.unwrap(y[i:i+2], period=period), x[i:i+2])
            add = np.ma.array([ period, 0., 0.], mask=[0,1,0])
        else:
            # interpolate needs sorted xp = np.unwrap(y[i:i+2], period=period)
            x_transit = np.interp(0, np.unwrap(y[i:i+2], period=period)[::-1], x[i:i+2][::-1])
            add = np.ma.array([ 0., 0., period], mask=[0,1,0])
        x_add = np.ma.array([x_transit]*3, mask=[0,1,0])

        x = np.ma.hstack((x[:i+1], x_add, x[i+1:]))
        y = np.ma.hstack((y[:i+1], add, y[i+1:]))
    return x, y

这段代码是为了和Thorsten Kranz的原始回答进行比较，主要是因为数据点的密度较低。

import matplotlib.pyplot as plt

x = np.linspace(0., 100., 100)
y = (x*0.03 + np.sin(x) * 0.1) % 1

#Thorsten Kranz: Make a masked array with jump points masked
abs_d_data = np.abs(np.diff(y))
mask = np.hstack([np.abs(np.diff(y))>.5, [False]])
masked_y = np.ma.MaskedArray(y, mask)

# Plot
plt.figure()
plt.plot(*periodic2plot(x, y, period=1), label='This answer')
plt.plot(x, masked_y, label='Thorsten Kranz')

plt.autoscale(enable=True, axis='both', tight=True)
plt.legend(loc=1)
plt.tight_layout()

Answer 2

要找出粒子穿过上边界的地方，你可以这样做：

>>> import numpy as np
>>> a = np.linspace(0, 10, 50) % 5
>>> a = np.linspace(0, 10, 50) % 5 # some sample data
>>> np.nonzero(np.diff(a) < 0)[0] + 1
array([25, 49])
>>> a[24:27]
array([ 4.89795918,  0.10204082,  0.30612245])
>>> a[48:]
array([ 4.79591837,  0.        ])
>>> 

np.diff(a) 这个函数会计算数组 a 中相邻元素之间的差值，而 np.nonzero 则用来找出那些差值小于零的地方，也就是说，粒子是向下移动的。

Answer 3

我会采取不同的方法。首先，我不会通过查看导数的符号来确定跳跃点，因为运动可能会向上或向下，甚至可能有一些周期性的变化。我会关注那些导数值最大的点。

其次，一个优雅的方法来在图表中制造断点，就是在每个跳跃点上遮罩一个值。这样，matplotlib会自动生成线段。我的代码是：

import pylab as plt
import numpy as np

xs = np.linspace(0., 100., 1000.)
data = (xs*0.03 + np.sin(xs) * 0.1) % 1

plt.subplot(2,1,1)
plt.plot(xs, data, "r-")

#Make a masked array with jump points masked
abs_d_data = np.abs(np.diff(data))
mask = np.hstack([ abs_d_data > abs_d_data.mean()+3*abs_d_data.std(), [False]])
masked_data = np.ma.MaskedArray(data, mask)
plt.subplot(2,1,2)
plt.plot(xs, masked_data, "b-")

plt.show()

这样得到的结果是：

当然，这样的缺点是每个断点会损失一个数据点，但考虑到你似乎有的采样率，我想你可以用这个方法换取更简单的代码。