找到一组微分方程的稳态答案

【问题标题】：Find steady-state of a set of differential equations找到一组微分方程的稳态
【发布时间】：2015-05-18 08:05:04
【问题描述】：

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

你有更有效的方法吗？

【问题讨论】：

这真的取决于稳态的定义！ IE。考虑一个具有恒定频率和恒定幅度的振荡电路。虽然它不符合标准定义 dx/dt = 0，但它在拉普拉斯域中仍然是一个常数……
是的，你完全正确。让我们假设一个非常保守的情况，x=const 直到小幅波动。
更加吹毛求疵（抱歉），但高效是什么意思？ - 快速还是误报很少？让我参考这篇讨论检测稳态问题的文章：prism.groups.et.byu.net/uploads/Members/kelly_jpc2013.pdf
Henrik：在我的特殊情况下，可以证明只有一个具有特定高度的稳定状态。因此，我对速度和精度方面的效率相当感兴趣！

标签： python scipy differential-equations

【解决方案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

【讨论】：

【解决方案2】：

参考你的cmets，我没有看到更好的方法：

在这种情况下，您大致知道系统将在哪里稳定。（这在某些系统中相当容易预测，例如钟摆或充电电容器），并且它会稳定下来，我知道的最快方法是检查是否

(p[0] * x[i] + p[1] * x[i-1] ... + p[n] * x[i-n]  - mean(x[i-0:n]) )< epsilon)

困难在于确定 epsilon 的大小，参数 p[0:n] 来调整此检测。

显然你已经在使用这个窗口方法了：

epsilon = varmax
p[0:n] = 1
n = 22

并通过删除过滤器的参数并简单地使用方差来优化它。

对于许多系统（其中的稳定看起来像一阶微分方程），滤波器参数看起来像这样：

p[0] = n/2
p[n] = n/2
p[1:n-1] = 0

意思是如果你用这个简单的测试代替状态方差的计算，你可以让事情变得更快：

if( abs(x[-22] - x[0]) > epsilon)

如果仍然存在小的干扰，这将无法正确检测到，因为我们不了解您的系统，这很难说...

【讨论】：