python - 找到一组微分方程的稳态

Question

让我们假设我有一组微分方程要与 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

您有更有效的方法吗？

Answer 1

如果您使用 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