limit cycles

hbghlyj

2.3.5 Another important example: limit cycles

使用Matlab程式pplane8绘图 蓝实线是orbit, 红虚线是nullcline 红点是equilibrium point(即critical point)

考虑plane autonomous system:
\begin{align*}
\dot{x}&=\left(1-\left(x^2+y^2\right)^{\frac{1}{2}}\right) x-y \\
\dot{y}&=\left(1-\left(x^2+y^2\right)^{\frac{1}{2}}\right) y+x .
\end{align*}
由方程看出，唯一的临界点是 $(0,0)$。线性近似:
Jacobian $\pmatrix{1&-1\\1&1}$, $\lambda_1=1+i$, $\lambda_2=1-i$.
特征向量: $v_1=\pmatrix{\frac1{\sqrt2}\\\frac{-i}{\sqrt2}},v_2=\pmatrix{\frac1{\sqrt2}\\\frac{i}{\sqrt2}}$
所以它是一个 unstable spiral.

在这种情况下，我们可以分析整个非线性系统。我们将转换为极坐标。最简单的方法如下：
令 $x^2+y^2=r^2$ 得
\begin{align*}
& X=x(1-r)-y \\
& Y=y(1-r)+x
\end{align*}
我们计算
\begin{align*}
r \dot{r}&=x \dot{x}+y \dot{y}=x[x(1-r)-y]+y[y(1-r)+x] \\
&=r^2(1-r)
\end{align*}
即
$$
\dot{r}=r(1-r) .
$$
然后从
$$
y=r \sin \theta,
$$
我们发现
$$
\dot{y}=\dot{r} \sin \theta+r \cos \theta \dot{\theta}=y(1-r)+x,
$$
给出了$\dot{\theta}$的表达式. 代入微分方程组变为
\begin{cases}
\dot{\theta}=1 \\
\dot{r}=r(1-r) .
\end{cases}Unlike the system in its previous form, we can solve this. First
$$
\theta=t+\text {const}
$$
and then
$$
\int d t=\int \frac{d r}{r(1-r)}=\int\left(\frac{1}{r}+\frac{1}{1-r}\right) d r
$$
so
$$
\log \frac{r}{|1-r|}=t+\text {const}
$$
i.e.
$$
\frac{r}{1-r}=A e^t
$$
Solve for $r$ and change the constant:
$$
r=\frac{1}{1+B e^{-t}}=\frac{1}{1+\left(\frac{1}{r_0}-1\right) e^{-t}}
$$
where $r(0)=r_0$.
Note that as $t \rightarrow \infty$, $r \rightarrow 1$,
as $t \rightarrow-\infty$ either $r \rightarrow 0$ if $r_0<1$ or $r \rightarrow \infty$ at some finite $t$ if $r_0>1$.

Now it is clear that the origin is an unstable spiral, and that the trajectories spiral out of it anticlockwise. We can also see that $r=1$ is a closed trajectory and that all other trajectories (except the fixed point at the origin) tend to it; we call such a closed trajectory a limit cycle. It is stable because the other trajectories converge on it. (For an example of an unstable limit cycle we could consider the same system but with $t$ changed to $-t$.

Another system with a limit cycle arises from the Van der Pol equation:
$$
\ddot{x}+\epsilon\left(x^2-1\right) \dot{x}+x=0
$$
where $\epsilon$ is a positive real constant. If $\epsilon=0$ this is the harmonic oscillator again. If $\epsilon \neq 0$ then the usual trick produces a plane autonomous system:
\begin{align*}
\dot{x}&=y\\
\dot{y}&=-\epsilon\left(x^{2}-1\right) y-x
\end{align*}
The only critical point is $(0, 0)$ and it’s an unstable spiral for $\epsilon > 0$ (exercise!).

Claim: Its beyond us to show this, but this system has a unique limit cycle, which is stable. There are some good illustrations for this in e.g. Boyce and DiPrima (pp. 496–500 of the 5th edition).
第八版PDF及Solution

hbghlyj

微分方程 \begin{align*}\frac{dx}{dt}&=p(x,y)\\\frac{dy}{dt}&=q(x,y)\end{align*} 构成的平面自治系统，它的解是由无穷多条互不相交的解曲线$(x,y)=(x(t),y(t))$构成。如果平面上的一点X是该微分方程的一个奇点，那么过X点附近任一点的解曲线在t趋向正无穷或者负无穷时都无限趋近X，或者环绕X周期性地旋转。解曲线趋向X的方式有图1所示的前三种可能，它们对应的X分别称为微分系统的焦点、结点、鞍点，而解曲线环绕X旋转时的X则称为微分系统的中心。

图1：微分系统的4类奇点及其附近解曲线的示意图

(a) 焦点	(b) 结点	(c) 鞍点	(d) 中心

如果平面上的一个闭环T是上述微分方程的一个孤立周期解，那么过T附近任意一点的解曲线在t趋向正无穷或者负无穷时都会无限逼近T。这时T就是微分系统的极限环：如果解曲线都在t趋向正无穷时逼近T，则称该极限环稳定；否则称其不稳定，如图2所示。

图2：稳定和不稳定极限环示意图

(a) 稳定极限环	(b) 不稳定极限环

微分系统的解曲线构成的向量场中会出现极限环,这一现象是由法国大数学家H.Poincaré首先观察发现的。极限环的出现让微分系统解向量场的性态和结构特征陡添神秘色彩。什么时候会出现极限环？会出现多少个？它们会处于怎样的相对位置？

德国大数学家D.Hilbert将这些问题作为第16问题的后半部分列入了他在1900年国际数学家大会上提出的、影响整个二十世纪数学发展的23个问题。Hilbert所提的问题是，当p(x,y)和q(x,y)为n次实系数多项式时，上述平面微分系统极限环的最大个数是多少？它们的相对位置关系如何？

我们来看两个具体例子。首先容易验证原点是三次平面系统 \begin{align*}\frac{dx}{dt}&=-y-x(x^2+y^2)\\\frac{dy}{dt}&=x-y(x^2+y^2)\end{align*} 的唯一孤立奇点，它是稳定焦点，见图4(a)。考虑该系统的扰动系统 \begin{align*}\frac{dx}{dt}&=kx-y-x(x^2+y^2)\\\frac{dy}{dt}&=x+ky-y(x^2+y^2)\end{align*} 当扰动参数k大于0时，扰动微分系统都会有一个稳定的极限环，见图4(b)。

图4：稳定焦点和不稳定焦点附近的极限环

(a) 稳定焦点	(b) 稳定极限环

通过适当扰动焦点型微分系统，通常可以得到具有极限环的扰动系统。这种通过扰动分岔出极限环的经典方法是由Poincaré、Lyapunov、Andronov等数学大家在十九、二十世纪发展起来的。该方法及其奠基理论——微分方程的定性理论，可以用来有效地分析微分方程解的稳定性、周期性、分岔等各种结构性质。结构性质的研究会相对容易如果我们知道解的解析表达式，可是在绝大多数情形微分方程的解都无法用解析表达式给出。所以微分方程的定性理论在解的形式未知时为分析解的性态提供了一条有效的途径。

微分方程的定性分析方法涉及到大量繁琐复杂的推理和演算，往往超出了以纸笔为工具的传统数学家推演能力的极限。很多经典的定性分析结果都未能幸免这样或那样的推演错误。特别值得一提的是，1955年俄国数学家I.G.Petrovskii和E.M.Landis声称证明了平面二次系统最多只有3个极限环。这个结果曾轰动一时，但其中的若干引理很快受到其他学者的质疑。1978年，我国数学家史松龄与陈兰荪、王明淑分别举出了具有4个极限环的实例，从而推翻了俄国数学家的结论。中国数学家令人称奇的这一发现将极限环问题的研究推向高潮，同时也展示了中国学者在这一领域的研究水平。1990年前后，王东明首先将特征列、Gröbner基等符号计算方法用于微分系统的定性分析，很大程度上化解了定性分析中的繁琐推演问题，并取得了若干研究成果，包括给出了一类具有6个极限环的平面三次系统。王东明的工作引发了大量基于符号计算软件和方法分析微分系统性态结构的后续研究。一大批著名学者包括N.G.Lloyd、H.Zoladek、李继彬、李承志、韩茂安等相继对极限环问题展开了深入研究，获得了一系列标志性的成果，从而使该领域的研究热度二十多年持续不减。与此同时，符号计算方法还被用于生物学、控制理论、系统科学等领域中的微分系统和离散动力系统的定性分析。

微分方程是描述非线性现象和规律的最基本数学工具，遍及数学、力学、生物学等众多科学和工程领域，其研究和应用范围都非常广泛。微分方程的定性分析始终会是十分重要的难题。随着符号计算和算法数学的不断发展、计算机计算能力的不断提高、计算机软件的不断改进，微分方程定性分析所涉及的推理和计算问题有望在未来的20年内取得重大突破。正在兴起的超大规模网络并行计算和以构造性数学为基础的计算智能有可能帮助揭开极限环的神秘面纱，为Hilbert极限环问题的最终解决提供系统化、机械化的工具和途径。

hbghlyj

\begin{cases}x'=-y+x(1-x^2-y^2)\\y'=x+y(1-x^2-y^2)\end{cases}
当初始数据为 $x(0)=1,y(0)=0$ 轨迹是单位圆
Asymptote normalized vector field
import graph;

size(300);
real xmin = -2;
real xmax = 2;
real ymin = -2;
real ymax = 2;

real dx = 0.4; // x grid spacing
real dy = 0.4; // y grid spacing

pen vectorPen = black+0.8bp; // pen for the direction vectors

pair field(pair P) {
  real x = P.x;
  real y = P.y;
  return unit((-y + x*(1 - x^2 - y^2), x + y*(1 - x^2 - y^2)))/2;
}

void drawField(real xmin, real xmax, real ymin, real ymax) {
  for (real x = xmin; x <= xmax; x += dx) {
    for (real y = ymin; y <= ymax; y += dy) {
      pair P = (x, y);
      pair F = field(P);
      draw(P -- P + F, vectorPen, Arrow);
    }
  }
}

drawField(xmin, xmax, ymin, ymax);
xaxis(BottomTop, Ticks("%", Size=3));
yaxis(LeftRight, Ticks("%", Size=3));


WolframAlpha
Slope field

Vector field

Integral curves

hbghlyj

LaTeXDraw.com
Phase portrait of Van-Der-Pol oscillator in TikZ
December 12, 2020
TeXLive.net

1. \documentclass[border=0.2cm]{standalone}
3. \usepackage{pgfplots}
4. \pgfplotsset{compat = newest}
6. % Define arrow's style
7. \usetikzlibrary{decorations.markings}
9. % Arrow style
10. \tikzset{decorated arrows/.style={
11.     postaction={
12.         decorate,
13.         decoration={
14.             markings,
15.             mark=between positions 0 and 1 step 15mm with {\arrow[black]{stealth};}
16.             }
17.         },
18.     }
19. }
22. \begin{document}
23. \begin{filecontents}{data.txt}
24. -1.4524 -2.0041
25. -1.5484 -1.8329
26. -1.6354 -1.6464
27. -1.7128 -1.4509
28. -1.7804 -1.2526
29. -1.8381 -1.0571
30. -1.8863 -0.86899
31. -1.9252 -0.69163
32. -1.9556 -0.52729
33. -1.9782 -0.37716
34. -1.9936 -0.24155
35. -2.0026 -0.12011
36. -2.0058 -0.012026
37. -2.004 0.083793
38. -1.9976 0.16859
39. -1.9873 0.24366
40. -1.9734 0.31025
41. -1.9564 0.36955
42. -1.9366 0.42265
43. -1.9142 0.47053
44. -1.8896 0.51406
45. -1.8629 0.55402
46. -1.8342 0.59108
47. -1.8038 0.62583
48. -1.7717 0.65879
49. -1.738 0.69041
50. -1.7027 0.72111
51. -1.6659 0.75124
52. -1.6275 0.78112
53. -1.5877 0.81105
54. -1.5464 0.8413
55. -1.5036 0.87212
56. -1.4592 0.90377
57. -1.4132 0.93647
58. -1.3655 0.97046
59. -1.3161 1.006
60. -1.2649 1.0432
61. -1.2118 1.0825
62. -1.1566 1.1239
63. -1.0994 1.1678
64. -1.0398 1.2144
65. -0.97786 1.264
66. -0.91335 1.3167
67. -0.84613 1.3729
68. -0.77601 1.4326
69. -0.7028 1.4961
70. -0.62633 1.5635
71. -0.54639 1.6348
72. -0.46279 1.7098
73. -0.37535 1.7884
74. -0.2839 1.8701
75. -0.18831 1.954
76. -0.088478 2.0392
77. 0.015613 2.1242
78. 0.1239 2.2068
79. 0.23622 2.2847
80. 0.35224 2.3547
81. 0.4715 2.4132
82. 0.59331 2.4563
83. 0.7168 2.4798
84. 0.8409 2.4798
85. 0.96433 2.4528
86. 1.0857 2.3965
87. 1.2035 2.3099
88. 1.3162 2.1938
89. 1.4224 2.0509
90. 1.5209 1.8853
91. 1.6106 1.7027
92. 1.691 1.5093
93. 1.7615 1.3112
94. 1.8221 1.1144
95. 1.873 0.92364
96. 1.9146 0.7428
97. 1.9475 0.57444
98. 1.9723 0.42004
99. 1.9897 0.28015
100. 2.0005 0.15459
101. 2.0054 0.042661
102. 2.005 -0.05666
103. 2 -0.14459
104. 1.9907 -0.2224
105. 1.9779 -0.29137
106. 1.9617 -0.35271
107. 1.9427 -0.40753
108. 1.9211 -0.45686
109. 1.8971 -0.50159
110. 1.871 -0.54253
111. 1.8429 -0.58038
112. 1.813 -0.61576
113. 1.7814 -0.6492
114. 1.7481 -0.68117
115. 1.7133 -0.71211
116. 1.6769 -0.74237
117. 1.639 -0.77229
118. 1.5997 -0.80218
119. 1.5588 -0.8323
120. 1.5164 -0.86293
121. 1.4725 -0.89431
122. 1.427 -0.92668
123. 1.3798 -0.96026
124. 1.3309 -0.9953
125. 1.2803 -1.032
126. 1.2277 -1.0706
127. 1.1732 -1.1114
128. 1.1165 -1.1546
129. 1.0577 -1.2004
130. 0.99643 -1.249
131. 0.9327 -1.3008
132. 0.8663 -1.3559
133. 0.79705 -1.4146
134. 0.72478 -1.4769
135. 0.64929 -1.5432
136. 0.5704 -1.6133
137. 0.4879 -1.6872
138. 0.40162 -1.7648
139. 0.31137 -1.8456
140. 0.21701 -1.929
141. 0.11844 -2.014
142. 0.015605 -2.0992
143. -0.091456 -2.1828
144. -0.20261 -2.2623
145. -0.31758 -2.335
146. -0.43594 -2.3973
147. -0.55707 -2.4454
148. -0.68017 -2.4751
149. -0.80422 -2.4825
150. -0.92799 -2.4637
151. -1.0501 -2.4163
152. -1.1691 -2.3386
153. -1.2835 -2.2311
154. -1.3917 -2.0957
155. -1.4926 -1.9363
156. \end{filecontents}
157. \begin{tikzpicture}
159. \begin{axis}[
160.     grid=both,
161.     grid style={dashed,red!20},
162.     xmin = -4, xmax = 4,
163.     ymin = -4, ymax = 4,
164.     width = \textwidth,
165.     height = 0.7\textwidth,
166.     xlabel = {$x$},
167.     ylabel = {$y$},
168.     title={Phase Portrait of Van Der Pol Oscillator},
169.     view = {0}{90},
170. ]
172. % Vector Field
173. \addplot3[
174.     quiver = {
175.         u = {y/sqrt(y^2+(0.8*(1-x^2)*y-x)^2)},
176.         v = {(0.8*(1-x^2)*y-x)/sqrt(y^2+(0.8*(1-x^2)*y-x)^2)},
177.         scale arrows = 0.25,
178.     },
179.     -stealth,
180.     domain = -4:4,
181.     domain y = -4:4,
182.     lightgray]
183. {0};
185. % Plot the Limit Cycle      
186. \addplot[ultra thick, red, decorated arrows] file {data.txt};
188. \end{axis}
189. \end{tikzpicture}
190. \end{document}

复制代码