在相空间中的轨迹 - Overflowerror：（34，“结果太大”）

Question

我虽然有C++和MATLAB的使用经验，但对Python还是新手。目前，我正在编写一个程序，用于绘制涉及dy/dt和dx/dt的非线性系统在相空间中的轨迹。然而，对于更复杂的函数形式，我遇到了Overflow Error错误。有没有办法可以解决这个问题呢？提前感谢！

以下是我的代码：

fig = plt.figure(figsize=(18,6))
dt = 0.01
def trajectories():
    #initial conditions
    x = y = 0.1
    xresult = [x]
    yresult = [y]
    for t in xrange(10000):
        # functional form: dx/dt = y, dy/dt = -r(x**2-1)*y-x
        nextx = x + (r*x-y+x*y**2) * dt
        nexty = y + (x + r*y + y**3) * dt
        x, y = nextx, nexty
        xresult.append(x)
        yresult.append(y)
    plt.plot(xresult, yresult)
    plt.axis('image')
    plt.axis([-3, 3, -3, 3])
    plt.title('r = ' + str(r))


rs = [-1, -0.1, 0, .1, 1]
for i in range(len(rs)):
    fig.add_subplot(1, len(rs), i + 1)
    r = rs[i]
    trajectories()
plt.show()

EDIT: this is the full traceback
Traceback (most recent call last):
  File "/Users/Griffin/Atom/NumInt.py", line 33, in <module>
    trajectories()
  File "/Users/Griffin/Atom/NumInt.py", line 18, in trajectories
    nextx = x + (r*x-y+x*y**2) * dt
OverflowError: (34, 'Result too large')

- griffinleow

---

请发布完整的回溯(traceback)信息。 - roganjosh

刚刚添加了完整的回溯信息。谢谢！ - griffinleow

你的邮政编码对我来说运行良好。我在考虑你的数据的小数精度...你的dt固定为0.01吗？对于更小的值会发生什么？ - eduardosufan

Answer 1

你遇到的直接错误与使用欧拉算法进行积分有关，因为在你使用的步长下，欧拉算法变得不稳定。实际上，根本问题是使用欧拉算法。下面的代码使用scipy.integrate.odeint来处理积分，并且由于能够使用可变步长而表现更好。一些积分仍然不完美，但至少我们得到了一些结果。

import numpy
import scipy.integrate
import matplotlib.pyplot as plt

def derivatives(states, t, r):
    x, y = states
    return [r*x - y + x*y**2,
            x + r*y + y**3]

def trajectories(r):
    initial_conditions = [0.1, 0.1]
    times = numpy.linspace(0, 100, 1000)
    result = scipy.integrate.odeint(derivatives, initial_conditions, times, args=(r,))
    xresult, yresult = result.T

    plt.plot(xresult, yresult)
    plt.axis('image')
    plt.axis([-3, 3, -3, 3])
    plt.title('r = ' + str(r))

fig = plt.figure(figsize=(18,6))

rs = [-1, -0.1, 0, .1, 1]

for i, r in enumerate(rs, 1):  # Avoid for i in range(len(rs))
    fig.add_subplot(1, len(rs), i)
    trajectories(r)
plt.show()

结果:

这段内容涉及IT技术，展示了一个模拟结果的截图。

- chthonicdaemon