0

我已经实现了以下方程组的解决方案

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 Forbidden 消息，请在 url 字段内按 enter）

后来我把同样的代码写进了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 Forbidden 消息，请在 url 字段中按 enter）

如果在较小的字段中再次绘制 Mathematica 解：

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

他将得到与 python 解决方案类似的结果。只有轴'会错位。

为什么剧情差别这么大？我究竟做错了什么？

我怀疑我对模型的 python 实现是错误的，尤其是在计算 f1 的地方。或者，在这种情况下， plot() 函数对于绘制参数方程可能根本不方便。

谢谢。

ps：很抱歉没有拍到文字里面的图片，让你的生活变得艰难；我还没有足够的声望。

4

1 回答 1

0

您t在输入向量中用作第三个参数，而不是作为单独的参数。tinf(z,t)从不使用；相反，您使用z[2]，这将不等于t您之前定义的范围 ( t=np.linspace(0,20.,1000))。for的lambda功能在g这里无济于事：您只使用一次来设置 a 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在评论中的内容，因为我觉得现在更容易看到事情是正确的。我得到的数字：

在此处输入图像描述

于 2014-06-03T09:51:29.770 回答