微分方程组

Question

我必须用Python中的ODEINT解两个微分方程，方程：

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

如果您阅读docs，您将看到odeint

Solves the initial value problem for stiff or non-stiff systems of first order ode-s:

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

where y can be a vector

这种转换是将二阶常微分方程转换为一阶矢量常微分方程的标准数学方法。你知道吗

因此，您创建了一个新的向量变量（我称之为Y以避免混淆），由向量Y = [y, y_prime, z, z_prime]组成：您的函数实现是正确的。你知道吗

还要注意的是，为了进行数值求解，需要指定所有向量的初始条件，在本例中为y0，z0，y'0和z'0。正如Thomas指出的，在调用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)