Дифференциальные уравнения - ODEINT

Question

Мне нужно решить два дифференциальных уравнения с помощью ODEINT в Python, уравнения:

y''(t) = (l*q)/a * (1/y(p) * [1 - z'(p)*u]
z''(t) = a * (1/y(p) * y'(p)*u 

Поэтому мне сказали сделать:

y1=y
y2=y'
z1=z
z2=z'

и

y1' = y2
y2' = y'' = (l*q)/a * (1/y(p) * [1 - z'(p)*u]
z1' = z2
z2' = z''(t) = a * (1/y(p) * y'(p)*u

и теперь я должен решить эти 4 уравнения.l, q, a, u известны.

Я пробовал что-то вроде этого:

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

def rownanie(y, t, l, q, a, u):
    y1, y2, z1, z2 = y
    dydt = [y2, ((l*q)/a)*(1/y1)*(1-z2*u), z2, (a*y2*u)/y1]
    return dydt

l = 1
q = 1
a = 10
u = 0.25

y0 = 0
z0 = 0
t = np.linspace(0, 10, 101)
sol = odeint(rownanie, y0, z0, t, args=(l,q,a,u))
print(sol)

Нужна помощь с этим

Answer 1

Если вы прочитаете документы , вы увидите odeint

Решает проблему начальных значений для жестких или не жестких систем од-с первого порядка:

dy / dt = func (y, t, ...) [или func (t, y, ...)]

, где y может быть вектором

Это преобразование является стандартным математическим способом преобразования ODE второго порядка в векторное ODE первого порядка.

Поэтому вы создаете новую векторную переменную (я буду называть ее Y, чтобы избежать путаницы), состоящую извектора Y = [y, y_prime, z, z_prime]: Ваша реализация функции верна.

Также обратите внимание, что для численного решения необходимо указать начальные условия для всех векторов, в данном случае y0, z0, y'0 и z'0.Как указал Томас, вам нужно указать эти значения в качестве начального значения вектора при вызове odeint.

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

def rownanie(Y, t, l, q, a, u):
    y1, y2, z1, z2 = Y
    dydt = [y2, ((l*q)/a)*(1/y1)*(1-z2*u), z2, (a*y2*u)/y1]
    return dydt

l = 1
q = 1
a = 10
u = 0.25

y0 = 0
z0 = 0
y0_prime, z0_prime = 0, 0  # you need to specify a value for these too
t = np.linspace(0, 10, 101)
sol = odeint(rownanie, [y0, y0_prime, z0, z0_prime], t, args=(l,q,a,u))
print(sol)