書記が数学やるだけ#877 リミットサイクルの存在条件
力学系について，リミットサイクルの性質を扱っていく。
問題
説明
リミットサイクルとは，力学系における相空間上での閉軌道であり，時間 t を無限大またはマイナス無限大にしたとき，その閉軌道に収束する軌道が少なくとも1つ存在するものである。物理的には自励振動の数理モデルとなり，リミットサイクルを持つ例としてファン・デル・ポール振動子がある。
閉軌道が存在するか否かについて，様々な判定法がある。
リアプノフ関数は，ロシアの数学者であるアレクサンドル・リアプノフにちなんで命名された関数であり，力学系や自励系を成す常微分方程式系における不動点の安定性を証明するために用いられる。
デュラックの判定法は，グリーンの定理に基づいた判定法である。
ポアンカレ・ベンディクソンの定理とは，平面上の連続力学系あるいは自励的常微分方程式系では，有界な軌道が時間経過後に最終的に落ち着く先は，平衡点を含まなければ周期軌道であることを述べる数学の定理である。
別の方法として，リエナールの定理を用いたものがある。
解答
まず「勾配系ならば閉軌道を持たない」例から。
リアプノフ関数について，リアプノフ関数が正定値関数で，その時間変化が負定値であるならば閉軌道を持たない。
デュラックの判定法による方法について。
ポアンカレ-ベンディクソンの定理を用いたリミットサイクルの存在を示す例，ここではトラッピング領域から示していく。
トラッピング領域とは一度入った軌道が外に出ない領域であり，境界でのベクトル場が内部に向かう，軌道が領域内に閉じ込められるといった条件で表せる。
Sel’kov modelは，解糖系の律速段階である「F6P + ATP → F1,6BP + ADP」のモデルであり，あるa,bにおいてリミットサイクルを持つ（基質濃度が振動する）ことが示せる。
Sel’kov modelにおいては，以下の領域がトラッピング領域であることが示せる。
ヌルクラインの交点がリペラーであるようなa,bでは，リミットサイクルが存在する。
a=0.08, b=0.6の例を示す。
参考：グラフのコード
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint, solve_ivp
def normalize_vectors(DX, DY):
magnitude = np.sqrt(DX**2 + DY**2)
magnitude[magnitude == 0] = 1 # ゼロ除算回避
return DX / magnitude, DY / magnitude
def plot_vector_field(system, x_range=(-3, 3), y_range=(-3, 3), title="Phase Portrait"):
x_vals = np.linspace(*x_range, 25)
y_vals = np.linspace(*y_range, 25)
X, Y = np.meshgrid(x_vals, y_vals)
DX, DY = system(X, Y)
DX, DY = normalize_vectors(DX, DY)
plt.figure(figsize=(7, 7))
plt.streamplot(X, Y, DX, DY, color='b', arrowsize=1, density=1)
plt.axhline(0, color='black', linewidth=0.8, linestyle='--')
plt.axvline(0, color='black', linewidth=0.8, linestyle='--')
plt.xlim(*x_range)
plt.ylim(*y_range)
plt.xlabel("x")
plt.ylabel("y")
plt.title(title)
plt.grid(True, linestyle='--', linewidth=0.5)
plt.show()
def system1(x, y):
return np.sin(y), x * np.cos(y)
def system2(x, y):
return -x + 4*y, -x - 3*y**2
def system3(x, y):
return 1.00000001*y, -x - 3*y**2
def system_glycolysis(t, z, a=0.08, b=0.6):
x, y = z
dx = -x + a*y + x**2 * y
dy = b - a*y - x**2 * y
return [dx, dy]
def plot_polar_trajectory():
def equations(state, t):
r, theta = state
return [r - r**3, 1]
initial_conditions_list = [[1.5, 0], [0.5, 0], [2.0, 0], [1.5, np.pi]]
t = np.linspace(0, 50, 10000)
plt.figure(figsize=(8, 6))
for initial_state in initial_conditions_list:
solution = odeint(equations, initial_state, t)
x_values = solution[:, 0] * np.cos(solution[:, 1])
y_values = solution[:, 0] * np.sin(solution[:, 1])
plt.plot(x_values, y_values, label=f"Initial r = {initial_state[0]}, θ = {initial_state[1]}")
plt.title("Phase Plot with Different Initial Conditions")
plt.xlabel("x")
plt.ylabel("y")
plt.grid(True)
plt.axis('equal')
plt.legend()
plt.show()
def plot_glycolysis_dynamics():
x_range, y_range = (-0.1, 3), (-0.1, 3)
def vector_field(x, y):
return -x + y + x**2 * y, 2 - y - x**2 * y
plot_vector_field(vector_field, x_range, y_range, "Phase Portrait of Glycolysis Model")
sol = solve_ivp(system_glycolysis, [0, 50], [1, 1], t_eval=np.linspace(0, 50, 1000))
plt.figure(figsize=(7, 5))
plt.plot(sol.t, sol.y[0], label='x(t)')
plt.plot(sol.t, sol.y[1], label='y(t)')
plt.xlabel('Time')
plt.ylabel('Concentration')
plt.title('Time Evolution of x and y')
plt.legend()
plt.grid(True, linestyle='--', linewidth=0.5)
plt.show()
# 各ベクトル場のプロット
plot_vector_field(system1, title="Phase Portrait of dx/dt = sin(y), dy/dt = x cos(y)")
plot_vector_field(system2, title="Phase Portrait of dx/dt = -x + 4y, dy/dt = -x - 3y^2")
plot_vector_field(system3, title="Phase Portrait of dx/dt = 1.00000001*y, dy/dt = -x - 3y^2")
plot_polar_trajectory()
plot_glycolysis_dynamics()本記事のもくじはこちら：
