让我们假设我有一组微分方程要与 scipy odeint 集成。现在我的目标是找到稳态(我选择了初始条件以使该状态存在)。目前我已经实现了类似的东西
cond = True
while cond:
x = integrate(interval = [0,t], steps = 200)
if var(x[-22::]) < maxvar:
cond = False
return mean(x)
else:
t*= 2
您有更有效的方法吗？
请您参考如下方法:
如果您使用 odeint ，那么您已经将微分方程编写为函数 f(x, t)(或者可能是 f(x, t, *args))。如果您的系统是自治的(即 f 实际上并不依赖于 t)，您可以通过求解 f(x, 0) == 0< 来找到平衡 为 x。例如，您可以使用 scipy.optimize.fsolve求解平衡。
下面是一个例子。它使用 "Coupled Spring Mass System" example from the scipy cookbook . scipy.optimize.fsolve 用于求平衡解 x1 = 0.5, y1 = 0, x2 = 1.5, y2 = 0。
from scipy.optimize import fsolve
def vectorfield(w, t, p):
"""
Defines the differential equations for the coupled spring-mass system.
Arguments:
w : vector of the state variables:
w = [x1, y1, x2, y2]
t : time
p : vector of the parameters:
p = [m1, m2, k1, k2, L1, L2, b1, b2]
"""
x1, y1, x2, y2 = w
m1, m2, k1, k2, L1, L2, b1, b2 = p
# Create f = (x1', y1', x2', y2'):
f = [y1,
(-b1 * y1 - k1 * (x1 - L1) + k2 * (x2 - x1 - L2)) / m1,
y2,
(-b2 * y2 - k2 * (x2 - x1 - L2)) / m2]
return f
if __name__ == "__main__":
# Parameter values
# Masses:
m1 = 1.0
m2 = 1.5
# Spring constants
k1 = 8.0
k2 = 40.0
# Natural lengths
L1 = 0.5
L2 = 1.0
# Friction coefficients
b1 = 0.8
b2 = 0.5
# Pack up the parameters and initial conditions:
p = [m1, m2, k1, k2, L1, L2, b1, b2]
# Initial guess to pass to fsolve. The second and fourth components
# are the velocities of the masses, and we know they will be 0 at
# equilibrium. For the positions x1 and x2, we'll try 1 for both.
# A better guess could be obtained by solving the ODEs for some time
# interval, and using the last point of that solution.
w0 = [1.0, 0, 1.0, 0]
# Find the equilibrium
eq = fsolve(vectorfield, w0, args=(0, p))
print "Equilibrium: x1 = {0:.1f} y1 = {1:.1f} x2 = {2:.1f} y2 = {3:.1f}".format(*eq)
输出是:
Equilibrium: x1 = 0.5 y1 = 0.0 x2 = 1.5 y2 = 0.0
