2x2 diff の解のパラメトリック プロット。Python、Mathematica のシステム

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}]

そこから次のプロットが得られます。 プロット1 (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()

そして、これが私が得たプロットです:プロット2 (403 Forbidden メッセージが表示される場合は、URL フィールド内で Enter キーを押してください)

Mathematica の解をより小さいフィールドに再度プロットすると:

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

Python ソリューションと同様の結果が得られます。軸だけがずれてしまいます。

なぜプロットにこれほど大きな違いがあるのでしょうか?私の何が間違っているのでしょうか？

モデルの Python 実装、特に f1 が計算される部分が間違っているのではないかと思います。あるいは、この場合のように、plot() 関数はパラメトリック方程式をプロットするのにまったく便利ではないかもしれません。

ありがとう。

ps:テキスト内の画像を平手打ちしなかったことにより、あなたの生活を困難にさせて申し訳ありません。まだ知名度が足りません。

Answer 1

あなたが使っているのは t 別個のパラメーターとしてではなく、入力ベクトルの 3 番目のパラメーターとして使用します。の t で f(z,t) 決して使用されません。代わりに、使用します z[2], の範囲と等しくありません。 t 前に定義したように (t=np.linspace(0,20.,1000)）。の lambda のための関数 g ここでは役に立ちません:設定するために一度だけ使用します。 t0, 、しかしその後は決して。

コードを簡素化し、入力ベクトル (ラムダ関数も同様) から 3 番目のパラメーターを削除します。例えば：

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, 、コメントに書いてあることは、物事が正しいことがより簡単にわかると思うので。私が得た図: