Параметрический график решения 2х2 разн. система на Python, Mathematica

Я реализовал решение следующей системы уравнений

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)

Позже я написал то же самое в 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)

Если снова построить график решения Mathematica в меньшем поле:

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

он получит результат, аналогичный решению на питоне. Только ось 'будет не на своем месте.

Почему такая большая разница в сюжетах? Что я делаю неправильно?

Я подозреваю, что моя реализация модели на Python неверна, особенно там, где вычисляется f1. Или, может быть, функция plot () совсем не удобна для построения параметрических уравнений, как в этом случае.

Спасибо.

ps: извините за то, что усложнил себе жизнь, не вставив изображения внутри текста; У меня пока недостаточно репутации.


python wolfram-mathematica differential-equations odeint parametric-equations
person milia    schedule 03.06.2014    source источник
comment
Возможно, это просто автоматическое масштабирование, которое Mathematica делает иначе, чем matplotlib. Последний попытается включить все точки данных внутри рисунка; легко увидеть, что Mathematica этого не делает (поскольку некоторые части строки обрезаны). Попробуйте масштабировать фигуру matplotlib до фигуры Mathematic и посмотрите, согласовано ли это, например plt.axis([-0.085, 0.085, -0.05, 0.07]) (что, по сути, является инверсией вашего упоминания о построении решения Mathematica в меньшем поле).   -  person    schedule 03.06.2014
comment
Спасибо за ответ. Я сделал это и получил этот результат (пожалуйста, нажмите Enter). Что меня беспокоит, так это возможность того, что расчеты в модели Python каким-то образом ошибочны. Мне просто кажется, что код на Python производит меньше оборотов спирали за тот же промежуток времени.   -  person milia    schedule 03.06.2014
comment
Опасаясь того, что x (t) и y (t) в моем коде Python вычисляются иначе, чем в коде Mathematica, я построил график x (t) для (t, 0,10) в Python и Mathematica. И я получил следующие два графика: Math, Pyth. Понятно, что функции очень разные.   -  person milia    schedule 03.06.2014
comment
Хорошо, я нашел ошибку! Функция f1 должна была быть f1=-g(t)*yi-xi. Также gi и g0 в этом случае не нужны. Теперь функция x и параметрические графики похожи на результаты Mathematica.   -  person milia    schedule 03.06.2014

Ответы (1)


arrow_upward
0
arrow_downward

Вы используете t в качестве третьего параметра во входном векторе, а не как отдельный параметр. t в f(z,t) никогда не используется; вместо этого вы используете z[2], который не будет равняться диапазону t, как вы определили ранее (t=np.linspace(0,20.,1000)). Функция lambda для g здесь не поможет: вы используете ее только один раз для установки t0, но никогда после.

Упростите свой код и удалите этот третий параметр из входного вектора (а также лямбда-функцию). Например:

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, что вы указали в комментариях, так как я чувствую, что теперь легче увидеть, что все правильно. Цифра у меня получается:

введите описание изображения здесь

person Community    schedule 03.06.2014