Python和Mathematica中2x2微分系统解的参数图

Question

我实现了一个方程组的解法。

dy/dt = -t*y(t) - x(t)
dx/dt = 2*x(t) - y(t)^3

y(0) = x(0) = 1.
0 <= t <= 20

首先是在Mathematica中实现的，之后又在Python中实现。

我在Mathematica中的代码是：

s = NDSolve[
{x'[t] == -t*y[t] - x[t], y'[t] == 2 x[t] - y[t]^3, x[0] == y[0] == 1},
{x, y}, {t, 20}]

ParametricPlot[Evaluate[{x[t], y[t]} /. s], {t, 0, 20}]

从这个代码中，我得到了以下的图像：Plot1（如果出现403禁止访问的提示，请在网址栏按回车键）

后来我把同样的内容用Python编写了：

import numpy as np
from scipy.integrate import odeint
import matplotlib.pyplot as plt

g = lambda t: t

def f(z,t):
    xi = z[0]
    yi = z[1]
    gi = z[2]

    f1 = -gi*yi-xi
    f2 = 2*xi-yi**3
    return [f1,f2]

# Initial Conditions
x0 = 1.
y0 = 1.
g0 = g(0)
z0 = [x0,y0,g0]
t= np.linspace(0,20.,1000)

# Solve the ODEs
soln = odeint(f,z0,t)
x = soln[:,0]
y = soln[:,1]

plt.plot(x,y)
plt.show()

这是我得到的图像： Plot2（如果出现403禁止访问的提示，请在网址栏按回车键）

如果再在一个更小的区域中绘制Mathematica的解：

ParametricPlot[Evaluate[{x[t], y[t]} /. s], {t, 0, 6}]

你会发现得到的结果和Python的解类似，只是坐标轴的位置会不对。

为什么这些图像之间会有这么大的差别呢？我哪里做错了？

我怀疑我在Python中实现的模型有问题，尤其是在计算f1的部分。或者可能是plot()函数在绘制参数方程时并不太好用。

谢谢。

附言：抱歉让你们看起来很麻烦，因为我没有把图片放在文本里；我还没有足够的声望。

Answer 1

你在输入向量中把 t 当作第三个参数用，而不是作为一个单独的参数。其实在 f(z,t) 中的 t 根本没有被用到；你用的是 z[2]，这和你之前定义的 t 的范围（t=np.linspace(0,20.,1000)）是不一样的。对于 g 的 lambda 函数也没什么帮助：你只用它来设置一次 t0，之后就没再用了。

简化你的代码，把输入向量中的第三个参数去掉（还有那个 lambda 函数）。比如：

import numpy as np
from scipy.integrate import odeint
import matplotlib.pyplot as plt

def f(z,t):
    xi = z[0]
    yi = z[1]

    f1 = -t*yi-xi
    f2 = 2*xi-yi**3
    return [f1,f2]

# Initial Conditions
x0 = 1.
y0 = 1.
#t= np.linspace(0,20.,1000)
t = np.linspace(0, 10., 100)

# Solve the ODEs
soln = odeint(f,[x0,y0],t)
x = soln[:,0]
y = soln[:,1]

ax = plt.axes()
#plt.plot(x,y)
plt.plot(t,x)
# Put those axes at their 0 value position
ax.spines['left'].set_position('zero')
ax.spines['bottom'].set_position('zero')
ax.spines['right'].set_color('none')
ax.spines['top'].set_color('none')
ax.xaxis.set_ticks_position('bottom')
ax.yaxis.set_ticks_position('left')
#plt.axis([-0.085, 0.085, -0.05, 0.07])
plt.show()

我把你想要的实际图像注释掉了，改成了画 x 和 t 的关系，就像你注释里写的那样，因为我觉得这样更容易看出结果是否正确。得到的图像是：