非线性动力学系统的平衡点求解与稳定性分析咨询

Hey there! Let's work through your nonlinear dynamical system to figure out the equilibrium points and their stability—since SymPy is spitting out overly long solutions, breaking this down step by step should make things easier.

First, Simplify the Independent Variable #

Notice the equation for (x_5) is completely decoupled from other variables:
[
\frac{dx_5}{dt} = a_5 - d_5x_5
]
Setting this to zero for equilibrium gives an immediate, simple solution:
[
x_5^* = \frac{a_5}{d_5}
]
This is a constant (assuming parameters (a_5, d_5) are fixed), so we can substitute this into all other equations right away to reduce the system's complexity.

Solve for Remaining Variables Step-by-Step #

Let's go through each equation with (x_5^* = \frac{a_5}{d_5}) substituted in:

1. Equation for (x_6)
   [
   \frac{dx_6}{dt} = a_6x_6 + a_7x_5^* - k_5x_4x_6 = 0
   ]
   Rearranging terms to isolate (x_6):
   [
   x_6(a_6 - k_5x_4) = -\frac{a_7a_5}{d_5}
   ]

• If (a_6 \neq k_5x_4), we can solve explicitly for (x_6):
  [
  x_6^* = \frac{a_7a_5}{d_5(k_5x_4^* - a_6)}
  ]
• If (a_6 = k_5x_4), the equation only holds if (-\frac{a_7a_5}{d_5} = 0) (which usually isn't the case for positive biological/physical parameters), so we can focus on the first scenario for most cases.

2. Equation for (x_3)
   [
   \frac{dx_3}{dt} = a_3x_6 - d_3x_3 - k_1x_1x_3 = 0
   ]
   Factor out (x_3) and solve:
   [
   x_3^(d_3 + k_1x_1^) = a_3x_6^* \implies x_3^* = \frac{a_3x_6^}{d_3 + k_1x_1^}
   ]

3. Equation for (x_2)
   [
   \frac{dx_2}{dt} = a_2 + k_4x_4 - d_2x_2 - k_8x_2x_3 = 0
   ]
   Factor out (x_2):
   [
   x_2^(d_2 + k_8x_3^) = a_2 + k_4x_4^* \implies x_2^* = \frac{a_2 + k_4x_4^}{d_2 + k_8x_3^}
   ]

4. Equation for (x_4)
   [
   \frac{dx_4}{dt} = a_4 + k_{10}x_3 - d_6x_2x_4 - d_4x_4 + k_6x_6 - k_3x_2 = 0
   ]
   Group terms with (x_4):
   [
   x_4^*(-d_6x_2* - d_4) = -a_4 - k_{10}x_3^* - k_6x_6^* + k_3x_2^*
   ]
   Solve for (x_4^):
   [
   x_4^ = \frac{a_4 + k_{10}x_3^* + k_6x_6^* - k_3x_2^*}{d_6x_2* + d_4}
   ]
   Note: The numerator must be non-negative (assuming (x_4) represents a concentration or positive quantity), which will constrain valid parameter sets.

5. Equation for (x_1)
   This is the one that's giving SymPy trouble, but let's rewrite it with (x_5^) substituted:
   [
   a_1 + k_7\left(\frac{a_5}{d_5}\right)\left(\frac{x_1}{K+x_1}\right) + k_2x_1x_2 - d_1x_1 - d_7x_1\left(\frac{a_5}{d_5}\right) = 0
   ]
   Multiply through by (K+x_1) to eliminate the fraction, which turns it into a quadratic equation in (x_1):
   [
   a_1(K+x_1) + \frac{k_7a_5}{d_5}x_1 + x_1\left(k_2x_2^ - d_1 - \frac{d_7a_5}{d_5}\right)(K+x_1) = 0
   ]
   You can use SymPy's solve function on this quadratic directly (instead of the original fractional equation) and filter for positive real roots (since (x_1) is likely a positive quantity). This should give a more manageable solution than the raw SymPy output for the original equation.

Stability Analysis #

To analyze stability, we need to compute the Jacobian matrix of the system at the equilibrium point, then check the eigenvalues of this matrix. If all eigenvalues have negative real parts, the equilibrium is asymptotically stable; if any has a positive real part, it's unstable.

Here's the Jacobian (J) for your system (with (x_5^* = \frac{a_5}{d_5}) already substituted):
[
J = \begin{bmatrix}
k_7x_5^{*\frac{K}{(K+x_1})^2} + k_2x_2^ - d_1 - d_7x_5^* & k_2x_1^* & 0 & 0 & k_7\frac{x_1^*}{K+x_1} - d_7x_1^ & 0 \
0 & -d_2 - k_8x_3^* & -k_8x_2^* & k_4 & 0 & 0 \
-k_1x_3^* & 0 & -d_3 - k_1x_1^* & 0 & 0 & a_3 \
0 & -d_6x_4^* - k_3 & k_{10} & -d_6x_2^* - d_4 & 0 & k_6 \
0 & 0 & 0 & 0 & -d_5 & 0 \
0 & 0 & 0 & -k_5x_6^* & a_7 & a_6 - k_5x_4^*
\end{bmatrix}
]

Key Observations for Stability:

• The (x_5) row has a single non-zero element: (-d_5). This means one of the eigenvalues is (-d_5), which is negative (assuming (d_5 > 0)), so this component is always stable.
• For the remaining 5x5 submatrix, calculating analytical eigenvalues is extremely complex due to the nested equilibrium expressions. Instead, you can:
  1. Plug in specific parameter values to compute numerical eigenvalues (use NumPy's linalg.eigvals function).
  2. Use numerical simulation (e.g., SciPy's solve_ivp) to track system trajectories starting near the equilibrium and see if they converge or diverge.
  3. For parameter ranges, use the Routh-Hurwitz criterion (though this is tedious for 5th-order systems) to derive stability conditions.

Tips for SymPy #

• Reduce the system first by substituting (x_5^*) before running solve—this cuts down on the complexity SymPy has to handle.
• When solving for (x_1), explicitly tell SymPy to return only positive real roots using solveset with domain restrictions, e.g.:

  from sympy import solveset, S, symbols
  x1 = symbols('x1', positive=True)
  solveset(your_quadratic_eq, x1, domain=S.Reals)
  

备注：内容来源于stack exchange，提问作者konyong