Решить уравнение в питоне за заданный интервал времени с заданным начальным условием

Question

Я хочу решить уравнение в Python за промежуток времени I = [0,10] с начальным условием (x_0, y_0) = (1,0) и значениями параметров μ ∈ {-2, -1, 0, 1, 2} с помощью функции

scipy.integrate.odeint

Затем я хочу построить решения (x (t; x_0, y_0), y (t; x_0, y_0)) вплоскость ху.

Первоначально заданная линейная система:

dx / dt = y, x (0) = x_0

dy / dt = - x - μy, y (0) =y_0

Пожалуйста, смотрите мой код ниже:

import numpy as np

from scipy.integrate import odeint

sol = odeint(myode, y0, t , args=(mu,1)) #mu and 1 are the coefficients when set equation to 0

y0 = 0

myode(y, t, mu) = -x-mu*y

def t = np.linspace(0,10, 101) #time interval

dydt = [y[1], -y[0] - mu*y[1]]

return dydt

Может ли кто-нибудь проверить, правильно ли я определил вызываемую функцию myode ?Эта функция оценивает правую часть ODE.

Также для этой строки кода появилось сообщение об ошибке синтаксиса

def t = np.linspace(0,10, 101) #time interval

, сообщающее, что синтаксис неверен.Должен ли я как-то использовать

for * in ** 

, чтобы избавиться от сообщения об ошибке?Если да, то как именно?

Я очень плохо знаком с Python и ODE.Может ли кто-нибудь помочь мне с этим вопросом?Большое спасибо!

Answer 1

myode должно быть определением функции, таким образом

def myode(u, t, mu): x,y = u; return [ y, -x-mu*y]

Массив времени представляет собой простое объявление / присваивание переменной, там не должно быть def. Поскольку система двумерная, начальное значение также должно иметь размерность два

sol = odeint(myode, [x0,y0], t, args=(mu,) )

Таким образом, минимальная модификация вашего скрипта

def myode(u, t, mu): x,y = u; return [ y, -x-mu*y]
t = np.linspace(0,10, 101) #time interval
x0,y0 = 1,0  # initial conditions
for mu in [-2,-1,0,1,2]:
    sol = odeint(myode, [x0,y0], t, args=(mu,) )
    x,y = sol.T
    plt.plot(x,y)
a=5; plt.xlim(-a,a); plt.ylim(-a,a)
plt.grid(); plt.show()

дача участка

Answer 2

Попробуйте использовать метод solve_ivp.

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

i = 0
u = [-2,-1,0,1,2]
for x in u:
    def rhs2(t,y):
       return [y[1], -1*y[0] - u[x]*y[1]]
    value = u[i]
    res2 = solve_ivp(rhs2, [0,10],  [1,0] , t_eval=[0,1,2,3,4,5,6,7,8,9,10],  method = 'RK45') 

    t = np.array(res2.t[1:-1])
    x = np.array(res2.y[0][1:-1])
    y = np.array(res2.y[1][1:-1])
    fig = plt.figure()
    plt.plot(t, x, 'b-', label='X(t)')
    plt.plot(t, y, 'g-', label='Y(t)')
    plt.title("u = {}".format(value))
    plt.legend(loc='lower right')
    plt.show() 
    i = i + 1

Вот метод solve_ivp Документация

Здесь очень похожая проблема случшее объяснение.