用Python绘制函数的n次迭代图

Question

我正在研究动态系统，特别是一个叫做逻辑斯蒂家族的函数 g(x) = cx(1-x)。我需要多次迭代这个函数，以便理解它的行为。对于给定的某个点 x_0，我可以顺利地进行迭代，但我想要绘制整个函数及其迭代的图，而不仅仅是一个点。为了绘制单个函数，我有以下这段代码：

import numpy as np
import scipy as sp
import matplotlib.pyplot as plt

def logplot(c, n = 10):
    dt = .001
    x = np.arange(0,1.001,dt)
    y = c*x*(1-x)
    plt.plot(x,y)
    plt.axis([0, 1, 0, c*.25 + (1/10)*c*.25])
    plt.show()

我想我可以通过一种比较繁琐的方法来解决这个问题，具体来说就是明确地创建每次迭代的范围列表，像这样：

def log(c,x0):
    return c*x0*(1-x)

def logiter(c,x0,n):
    i = 0
    y = []
    while i <= n:
        val = log(c,x0)
        y.append(val)
        x0 = val
        i += 1
    return y

但这样做似乎太麻烦了，我在想有没有更好的方法。谢谢！

Answer 1

我想要达到的目标是通过可视化的方式，间接理解函数 g(c, x) = cx(1-x) 的初始条件的表现。

def jam(c, n):

    x = np.linspace(0,1,100)
    y = c*x*(1-x)
    for i in range(n):
        plt.plot(x, y)
        y = c*y*(1-y)

    plt.show()

Answer 2

一些不同的选项

这其实是风格的问题。你的解决方案有效，而且不难理解。如果你想继续这个方向，我建议稍微调整一下：

def logiter(c, x0, n):
    y = []
    x = x0
    for i in range(n):
        x = c*x*(1-x)
        y.append(x)

    return np.array(y)

修改的地方：

• 使用for循环比while循环更容易阅读
• x0在迭代中没有被使用（这增加了一个变量，但从数学上讲更容易理解；x0是一个常量）
• 函数写得很简单，直接用一行代码表示（如果不是这样，它的名字应该改成别的，不要用log，因为容易和对数混淆）
• 结果被转换成了numpy数组。（这通常是我需要绘图时的做法）

在我看来，这个函数现在足够清晰了。

你也可以采用面向对象的方法，创建一个逻辑函数对象：

class Logistics():
    def __init__(self, c, x0):
        self.x = x0
        self.c = c

    def next_iter(self):
        self.x = self.c * self.x * (1 - self.x)
        return self.x

然后你可以这样使用：

 def logiter(c, x0, n):
    l = Logistics(c, x0)
    return np.array([ l.next_iter() for i in range(n) ])

或者你可以把它做成生成器：

def log_generator(c, x0):
    x = x0
    while True:
        x = c * x * (1-x)
        yield x

def logiter(c, x0, n):
    l = log_generator(c, x0)
    return np.array([ l.next() for i in range(n) ])

如果你需要性能，并且有大表格的话，我建议：

def logiter(c, x0, n):
    res = np.empty((n, len(x0)))
    res[0] = c * x0 * (1 - x0)
    for i in range(1,n):
        res[i] =  c * res[i-1] * (1 - res[i-1])
    return res    

这样可以避免慢慢地转换成np.array和一些数据的复制。内存只分配一次，避免了从列表转换成数组的高开销。

（顺便说一下，如果你返回一个数组，初始的x0作为第一行，最后的版本会看起来更整洁。现在如果想避免复制向量，第一行必须单独计算。）

---

哪种方法最好？我不知道。在我看来，所有的方法都能读懂，都是合理的，这只是风格问题。不过，我的Python水平很差，所以可能还有其他更好的理由，或者上面的方法中有些并不好！

性能

关于性能：我在我的机器上尝试了以下方法：

logiter(3.2, linspace(0,1,1000), 10000)

前面三种方法的时间基本相同，大约是1.5秒。最后一种方法（预分配数组）的运行时间是0.2秒。不过，如果去掉从列表转换成数组的步骤，第一种方法的运行时间是0.16秒，所以时间主要花在了转换上。

可视化

我能想到两种有用但完全不同的可视化方法。你提到你会有大约100或1000个不同的x0作为起始值。你没有提到想要多少次迭代，但我们可以先从100次开始。所以，让我们创建一个包含100个不同x0和100次迭代的数组，c = 3.2。

data = logiter(3.6, np.linspace(0,1,100), 100)

一种标准的方法是画100条线，每条线代表一个起始值。这很简单：

import matplotlib.pyplot as plt

plt.plot(data)
plt.show()

这样得到：

看起来所有的值都在某个地方震荡，但除此之外，我们只有一团颜色。这种方法可能在你使用更窄的x0值范围时更有用：

data = logiter(3.6, np.linspace(0.8,0.81,100), 100)

你可以通过例如：

color1 = np.array([1,0,0])
color2 = np.array([0,0,1])
for i,k in enumerate(np.linspace(0, 1, data.shape[1])):
    plt.plot(data[:,i], '.', color=(1-k)*color1 + k*color2)

来给起始值上色，第一列（对应x0 = 0.80）用红色，最后一列用蓝色，中间用渐变色。（请注意，越蓝的点表示绘制得越晚，因此蓝色会覆盖红色。）

不过，也可以采取完全不同的方法。

data = logiter(3.6, np.linspace(0,1,1000), 50)
plt.imshow(data.T, cmap=plt.cm.bwr, interpolation='nearest', origin='lower',extent=[1,21,0,1], vmin=0, vmax=1)
plt.axis('tight')
plt.colorbar()

得到：

这是我个人最喜欢的方法。我不想过多解释，以免破坏大家的乐趣，但在我看来，这种方式很容易展示出行为的许多特性。