自励振動と非線形振動

自励振動と非線形振動　Ver. ‘20.9.14.

1　線形振動

1.1　本稿の目的・目標

本稿の目的は、以下の3つである。^※1※2

1. 摩擦による自励振動の発現メカニズムを説明する。
2. 振動系において、線形と非線形の違い、保存系と散逸系（非保存系）の違いを説明する。
3. 非線形振動の近似解を求める方法のひとつとして平均法を説明する。

一体系の非線形振動を考察する。文献ではしばしば運動するベルト上での振動が検討される（図1）が、ここでは後の都合で、摩擦面上で外力により運動する物体（図2では車。外力源と呼ぶこととする。）にバネを介してつながれた質点 \(m\) の振動状態を解析する（図2）。摩擦力は、質量 \(m\) の速度に対し非線形項を持つとする。

図2の系は、非線形（↔線形）かつ散逸系（↔保存系）である。本稿の流れは、線形・保存系をはじめとして、保存系から散逸系での特徴、線形から非線形での特徴を説明し、最後に非線形・散逸系として図2の系を解く。

1.2　線形保存振動＝調和振動

非線形という言葉があるなら線形という言葉が先にあるはずということで、まずは線形振動を考察する。また、まずは保存系を考えます。保存振動系という言葉が一般的か、正しいのかは知りません。とりあえずここでの方言と思ってください。エネルギーが保存している線形微分方程式にしたがう振動系という意味です。

あえて図1に似せて書くと、図3となる。注意してほしいのは、（1）摩擦が0ということ、（2）外力が働いていないこと。縦にぶら下げると図4。図3と図4は少しだけ状態が異なる。図4ではバネの伸び方向に重力が働いているため、振動の中心が自然長より \(mg/k\) だけ長くなっており、その分エネルギーを蓄えた状態になっている。

重力があっても調和振動になることは、次のようにわかる。まず図3の運動方程式は次式である。時間の1回微分を \(\dot{x}\)、2階微分を \(\ddot{x}\) で表す。

$$
m\ddot{X}=-k(X-L_0) \tag{1}
$$

\(x=X-L_0\) と変数変換すると、\(\dot{x}=\dot{X}\)、\(\ddot{x}=\ddot{X}\) なので、

$$
m\ddot{x}=-kx \tag{2}
$$

図4の運動方程式は次式である。下向きを正とする。

$$
m\ddot{X}=-k(X-L_0)+mg \tag{3}
$$

$$
m\ddot{X}=-k\left(X-L_0-\frac{mg}{k}\right) \tag{4}
$$

\(x=X-L_0-mg/k\) と変数変換すると、同様に \(\dot{x}=\dot{X}\)、\(\ddot{x}=\ddot{X}\) なので、

$$
m\ddot{x}=-kx \tag{5}
$$

となり、式(2)と同じ形となっている。ただし、振動の中心はずれている。

これら図3、図4の系は、エネルギーが保存されている。減りも増えもせず一定である。このことを今後の議論に役立つエネルギー積分の方法で議論しておく。式(5)の両辺に \(\dot{x}\) をかけて積分すると、

$$
m\frac{d}{dt}\int \ddot{x}\dot{x}\,dt
=
-k\frac{d}{dt}\int x\dot{x}\,dt \tag{6}
$$

$$
\frac{d}{dt}
\left\{
\frac{1}{2}m\dot{x}^{2}
+
\frac{1}{2}kx^{2}
\right\}
=0 \tag{7}
$$

となり、運動エネルギーとばねの弾性エネルギーの和が、時間の1階微分で0、すなわち全エネルギーは時間変化せず一定値を保つことがわかる。

さて、運動方程式(2)を解いて振動の状態を見てみよう。式(2)の解は、

$$
x=a\sin(\omega_0t)+b\cos(\omega_0t) \tag{8}
$$

$$
x=A_0\sin(\omega_0t+\phi_0) \tag{9}
$$

ここで以下のように定義した。^※3

$$
\omega_0\equiv\sqrt{\frac{k}{m}},
\qquad
A_0\equiv\sqrt{a^2+b^2},
\qquad
\cos\phi_0\equiv\frac{a}{A_0} \tag{10}
$$

\(\omega_0\) は系、すなわち物体の質量 \(m\) とバネ定数 \(k\) によって決まるため、系に固有のという意味を込めて固有振動数と呼ぶ。\(\phi\) は初期位相と呼ぶ。位相は初期条件、すなわち物体の位置と速度に依存するため「固有」という言葉はつけない。

詳細な解き方は教科書を見てもらうとするけど、次のことを覚えておくと理解しやすい。

• 2階微分すると元に戻ってマイナスがつく初等関数は、\(\sin\) と \(\cos\) である。
• 2階微分方程式なので、決まらない定数が2つ、\(a\) と \(b\)、あるいは \(A_0\) と \(\phi_0\) ある。
• 三角関数の公式 \(\sin(\alpha+\beta)=\sin\alpha\cos\beta+\cos\alpha\sin\beta\) を使う。

1.3　座標の設定と基礎方程式

図2について、質量 \(m\) の運動方程式を立てる。

$$
m\ddot{X}
=
-k(X-X_0-L_0)-f(\dot{X}) \tag{11}
$$

以下のように変数変換する。

$$
x=X-X_0-L_0,
\qquad
x_0=X_0 \tag{12}
$$

すると、時間の1階微分について、

$$
\dot{x}=\dot{X}-\dot{X}_0,
\qquad
\dot{x}_0=\dot{X}_0 \tag{13}
$$

時間の2階微分について、

$$
\ddot{x}=\ddot{X}-\ddot{X}_0,
\qquad
\ddot{x}_0=\ddot{X}_0 \tag{14}
$$

となる。まず、\(f(\dot{X})=f(\dot{x}+\dot{x}_0)\) を \(\dot{x}=\dot{x}_0\) のまわりでテイラー展開すると、式(11)は、

$$
m\ddot{X}
=
-k(X-X_0-L_0)
–
\sum_{j=0}^{\infty}
\frac{1}{j!}
f^{(j)}(\dot{x}_0)\dot{x}^{j} \tag{15}
$$

$$
m\ddot{X}
=
-k(X-X_0-L_0)
-f^{(0)}
–
\sum_{j=1}^{\infty}
\frac{1}{j!}
f^{(j)}(\dot{x}_0)\dot{x}^{j} \tag{16}
$$

ここで、

$$
\left.
\frac{d^n}{d\dot{x}^{\,n}}f(\dot{x})
\right|_{\dot{x}=\dot{x}_0}
=
f^{(n)}(\dot{x}_0) \tag{17}
$$

と略した。微分の回数を表す括弧ありの上付き文字と、べき乗を表す括弧なしの上付き文字を混同しないこと。

\(\dot{x}_0=\mathrm{const.}\) の場合を考える。このとき \(\ddot{x}_0=0\) である。式(12)、式(13)、式(14)の変数変換を式(16)に代入すると、

$$
\begin{aligned}
m\ddot{x}
&=
-kx-f^{(0)}
–
\sum_{j=1}^{\infty}
\frac{1}{j!}f^{(j)}(\dot{x}_0)\dot{x}^{j} \\
&=
-k\left(x+\frac{f^{(0)}}{k}\right)
–
\sum_{j=1}^{\infty}
\frac{1}{j!}f^{(j)}(\dot{x}_0)\dot{x}^{j}
\end{aligned}
\tag{18}
$$

ここで、

$$
\tilde{x}=x+\frac{f^{(0)}}{k} \tag{19}
$$

とおくと、\(\dot{\tilde{x}}=\dot{x}\)、\(\ddot{\tilde{x}}=\ddot{x}\) なので、式(18)は、

$$
m\ddot{\tilde{x}}
=
-k\tilde{x}
–
\sum_{j=1}^{\infty}
\frac{1}{j!}
f^{(j)}(\dot{x}_0)
\dot{\tilde{x}}^{\,j} \tag{20}
$$

\(\tilde{x}\) を改めて \(x\) と書き直して、

$$
m\ddot{x}
=
-kx
–
\sum_{j=1}^{\infty}
\frac{1}{j!}
f^{(j)}(\dot{x}_0)\dot{x}^{j} \tag{21}
$$

^※4

1.4　線形減衰振動系

1.4.1　解析的な解

式(21)で第一近似として和の部分を \(j=1\) までで打ち切ると、

$$
m\ddot{x}
+
f^{(1)}(\dot{x}_0)\dot{x}
+
kx
=
0 \tag{22}
$$

となる。\(\dot{X}=0\) の極限で図2の左、すなわち負の向きに摩擦力が働くので、自然長より縮んだ位置が振動の中心となる。重力がかかる場合の図4と比較しよう。

\(f^{(1)}(\dot{x}_0)=\mu\) と置き、全体を \(m\) で割って、

$$
\zeta\equiv\frac{\mu}{2\sqrt{mk}},
\qquad
\omega_0\equiv\sqrt{\frac{k}{m}} \tag{23}
$$

と定義する。^※5

$$
\ddot{x}
+
2\epsilon\zeta\omega_0\dot{x}
+
\omega_0^2x
=
0 \tag{24}
$$

\(\epsilon\) は小さいことを示す印として形式的に付与した。計算の最後に1とする。

解を \(x=\exp(\lambda t)\) と置き、式(24)に代入すると、

$$
\left(
\lambda^2
+
2\epsilon\zeta\omega_0\lambda
+
\omega_0^2
\right)
\exp(\lambda t)
=
0
$$

となるため、括弧内が0のとき式(24)が満たされる。\((\epsilon\zeta)^2-1<0\) のとき、振動しながら減衰する。以下この場合を考察する。

\((\epsilon\zeta)^2-1\geq0\) の場合は、振動せずに減衰する。ここでは省略する。

このとき、

$$
\lambda_1
=
\omega_0
\left(
-\epsilon\zeta
+
i\sqrt{1-(\epsilon\zeta)^2}
\right)
$$

$$
\lambda_2
=
\omega_0
\left(
-\epsilon\zeta
–
i\sqrt{1-(\epsilon\zeta)^2}
\right)
$$

が解となる。ここで \(i=\sqrt{-1}\) は虚数単位である。よって式(24)の一般解は、

$$
x
=
a_1\exp(\lambda_1t)
+
a_2\exp(\lambda_2t) \tag{25}
$$

$$
x
=
\exp(-\epsilon\omega_0\zeta t)
\left\{
a_1
\exp\left(
i\omega_0\sqrt{1-(\epsilon\zeta)^2}\,t
\right)
+
a_2
\exp\left(
-i\omega_0\sqrt{1-(\epsilon\zeta)^2}\,t
\right)
\right\} \tag{26}
$$

$$
\begin{aligned}
x
=
\exp(-\epsilon\omega_0\zeta t)
\biggl\{
&(a_1+a_2)
\cos\left(
\omega_0\sqrt{1-(\epsilon\zeta)^2}\,t
\right)\\
&+
i(a_2-a_1)
\sin\left(
\omega_0\sqrt{1-(\epsilon\zeta)^2}\,t
\right)
\biggr\}
\end{aligned}
\tag{27}
$$

$$
x
=
\exp(-\epsilon\omega_0\zeta t)
\left\{
a\cos\left(
\omega_0\sqrt{1-(\epsilon\zeta)^2}\,t
\right)
+
b\sin\left(
\omega_0\sqrt{1-(\epsilon\zeta)^2}\,t
\right)
\right\} \tag{28}
$$

$$
\begin{aligned}
x
=
\sqrt{a^2+b^2}\,
\exp(-\epsilon\omega_0\zeta t)
\biggl[
&
\frac{a}{\sqrt{a^2+b^2}}
\sin\left(
\omega_0\sqrt{1-(\epsilon\zeta)^2}\,t
\right)\\
&+
\frac{b}{\sqrt{a^2+b^2}}
\cos\left(
\omega_0\sqrt{1-(\epsilon\zeta)^2}\,t
\right)
\biggr]
\end{aligned}
\tag{29}
$$

$$
x
=
A\exp(-\epsilon\omega_0\zeta t)
\sin\left(
\omega_0\sqrt{1-(\epsilon\zeta)^2}\,t+\phi
\right) \tag{30}
$$

\(a_1\)、\(a_2\) は一般に複素数である。また \(x\) は実数なので、\(a\equiv i(a_2-a_1)\)、\(b\equiv a_2+a_1\) が実数となるように選んだ。オイラーの公式 \(\exp(i\theta)=\cos\theta+i\sin\theta\) を使った。

散逸がない場合は減衰しないが、散逸がある場合は \(\exp(-\epsilon\omega_0\zeta t)\) のため減衰していくことがわかる。

減衰時の振動数、すなわち減衰固有振動数は、

$$
\omega_d
\equiv
\omega_0\sqrt{1-(\epsilon\zeta)^2}
$$

であり、周期は \(T=2\pi/\omega_d\) となる。どちらも自由振動の場合と異なることに注意のこと。

1周期で振幅が減衰する割合は、ある時刻 \(t=t_0\) での振幅と、\(T\) 経過した \(t=t_0+T\) での振幅との比を取って、

$$
\frac{
\exp\{-\epsilon\omega_0\zeta(t_0+T)\}
}{
\exp(-\epsilon\omega_0\zeta t_0)
}
=
\exp(-\epsilon\omega_0\zeta T)
=
\exp\left(
-\frac{2\pi\epsilon\zeta}
{\sqrt{1-(\epsilon\zeta)^2}}
\right) \tag{31}
$$

この対数を取ったものを対数減衰率といい、

$$
\delta
\equiv
\frac{2\pi\epsilon\zeta}
{\sqrt{1-(\epsilon\zeta)^2}}
=
\frac{2\pi}
{\sqrt{(1/\epsilon\zeta)^2-1}} \tag{32}
$$

となる。\(\epsilon\zeta\ll1\) のときは近似的に \(2\pi\epsilon\zeta\) となる。^※6

注意してほしいのは、減衰率あるいは対数減衰率は各系での1周期で定義してあることである。異なる系の減衰度合いを比較するときは、同じ経過時間での振幅を比較する必要があるが、系が異なれば一般に周期 \(T\) も異なるので、対数減衰率を比較の指標とすることはできない。

\(t=t_0\) から \(t=t_0+t_1\) まで経過したときを考えると、振幅は初期 \(t_0\) の、

$$
\begin{aligned}
\exp(-\epsilon\omega_0\zeta t_1)
&=
\exp\left[
-\epsilon
\sqrt{\frac{k}{m}}
\frac{\mu}{2\sqrt{mk}}t_1
\right]\\
&=
\exp\left(
-\frac{\epsilon\mu}{2m}t_1
\right)
\end{aligned}
\tag{33}
$$

倍になっている。\(m\) と \(\mu\) に依存するが \(k\) には依存しないことに注意する。質量が大きい方が減衰は小さい。バネ定数には依存しない。^※7 \(\mu\) が大きい方が減衰は大きい。

パラメータの異なる系の比較を図6に示す。

減衰振動における振動状態を表すパラメータをまとめておく。

• 振動数（減衰固有振動数）：\(\omega_d=\omega_0\sqrt{1-(\epsilon\zeta)^2}\)。減衰なしと比較して小さくなる。
• 周期：\(T=2\pi/\omega_d\)。減衰なしと比較して大きくなる。
• 対数減衰率：\(\delta=2\pi\epsilon\zeta/\sqrt{1-(\epsilon\zeta)^2}\)。\(\epsilon\zeta\ll1\) のときは近似的に \(2\pi\epsilon\zeta\) となる。

1.4.2　結果の考察、特に振幅変化率について

ここで、振幅変化率について考察する。まずは、先の線形の例で対数減衰率についてもう少し考える。図6を見てみるとわかるが、各周期での振動の変化幅と比較して振幅の大きさの変化幅は小さい。

ここでは、このような振幅の変化が振動周期と比較して「ゆっくり」動く系を考える。でも「ゆっくり」では雰囲気しかわからないので、「ゆっくり」という言葉を定量化しよう。

まず、振幅の大きさの変化が振幅の大きさと同じくらい、代表的な値として減衰前と減衰後の差 \(1-\exp(-1)\sim0.63\) 程度になる時間間隔を \(t_d\) とすると、

$$
t_d
=
\frac{1}{\epsilon\omega_0\zeta}
=
\frac{2m}{\epsilon\mu}
$$

となる。この時間間隔が振動周期と比較して1より十分大きければ、ゆっくりと言えるだろう。結果として、「ゆっくり」を式で書くと、

$$
\frac{T}{t_d}
=
2\pi\epsilon\zeta
\ll1
$$

となる。

もう一点大事な点として、振動数が \(\omega_0\) からずれていることは注意してほしい。ただし \(\epsilon\) が小さいときは、そのずれも小さい。

1.4.3　一般的な系での振幅変化率

ここまで、線形減衰振動系が解析的に解けることを利用して考察してきた。しかし、微分方程式を解析的に解けるケースは、特に非線形な式では限られている。一般的な系での振幅変化率について考えてみる。

1.4では調和振動、すなわち非摂動系に減衰項を加えたが、代わりに弱い摂動項 \(f\) が加わった系を考える。摂動項は物体の速度 \(\dot{x}\) 以外に位置 \(x\) にも依存すると仮定する。

$$
\ddot{x}
+
\omega_0^2x
+
\frac{\epsilon}{m}f(x,\dot{x})
=
0 \tag{34}
$$

\(\epsilon\) は小さいことを表す無次元のパラメータで、実際の物理的なパラメータと考えてもよいし、\(f(x,\dot{x})\) で物理を表せるのなら計算の最後に1としてもよい。^※8

\(\epsilon=0\) のときの解は、調和振動子の式(9)となる。式(34)の一般解は、摂動項 \(\epsilon f(x,\dot{x})\) の形が具体的でないため、これ以上式には書くことができない。

さて、振幅の大きさ、位相の時間変化の大きさを見ておこう。\(A(t+T_0)\)、\(\phi(t+T_0)\) を \(t\) の周りで展開する。

$$
A(t+T_0)
=
A(t)
+
\dot{A}(t)T_0
+
\frac{1}{2}\ddot{A}(t)T_0^2
+\cdots
\sim
A(t) \tag{35}
$$

$$
\phi(t+T_0)
=
\phi(t)
+
\dot{\phi}(t)T_0
+
\frac{1}{2}\ddot{\phi}(t)T_0^2
+\cdots
=
\phi(t)+2\pi \tag{36}
$$

\(T\) ではなく \(T_0\)、2行目の右辺が0でなく \(2\pi\) に注意。

$$
\dot{A}(t)T_0\ll A(t),
\qquad
\frac{1}{2}\ddot{A}(t)T_0^2\ll A(t) \tag{37}
$$

$$
\dot{\phi}(t)T_0\simeq2\pi,
\qquad
\frac{1}{2}\ddot{\phi}(t)T_0^2\ll2\pi \tag{38}
$$

1.5　負勾配を持つ摩擦係数による自励振動

式(16)において \(f^{(1)}\) が負、つまり速度が大きいほど摩擦力が小さくなる場合を考える。まずは \(f^{(j)}=0\) for \(j>1\) とする。すると、わずかでも平衡点からずれると、減衰の逆で振幅が増幅していくことになる。

摩擦係数が負というのは、普通の感覚だと想像しづらいが、例えば濡れた平面上での運動などでは実現される。いわゆるストライベック線図と呼ばれるものを図7に示す。混合潤滑領域では、横軸の速度に対して摩擦係数の傾きが負になっている。

2　近似的な解法＝平均法

2.1　平均法の基礎式

線形振動をひととおり見てきた。非線形振動に続きたいわけだが、解析的な解を求められる場合は限られる。よって近似解の求め方のひとつとして、平均法について見ていくことにする。

もとの運動方程式は先に出てきた式である。

$$
\ddot{x}
+
\omega_0^2x
+
\frac{\epsilon}{m}f(x,\dot{x})
=
0
\qquad\text{再掲} \tag{39}
$$

非摂動系は \(\epsilon=0\) のときで、その解と時間微分は、

$$
x
=
A_0\sin(\omega_0t+\phi_0)
\qquad\text{再掲} \tag{40}
$$

$$
\dot{x}
=
A_0\omega_0\cos(\omega_0t+\phi_0) \tag{41}
$$

平均法では、非摂動系の解とその微分において未定定数とした振幅 \(A_0\)、位相 \(\phi_0\) を時間の関数とする。^※9

$$
x
=
A(t)\sin\{\omega_0t+\phi(t)\}
\equiv
A(t)\sin\psi(t) \tag{42}
$$

$$
\dot{x}
=
A(t)\omega_0\cos\{\omega_0t+\phi(t)\}
\equiv
A(t)\omega_0\cos\psi(t) \tag{43}
$$

見やすさのため、\(\psi(t)=\omega_0t+\phi(t)\) とおいた。

式(42)の時間微分をすべて書き出すと、

$$
\dot{x}
=
A(t)\omega_0\cos\psi(t)
+
A(t)\dot{\phi}(t)\cos\psi(t)
+
\dot{A}(t)\sin\psi(t) \tag{44}
$$

式(44)が式(43)と一致するためには、式(43)の第1項しか要らないから、

$$
\dot{x}_{\mathrm{re}}(t)
\equiv
A(t)\dot{\phi}(t)\cos\psi(t)
+
\dot{A}(t)\sin\psi(t)
=
0 \tag{45}
$$

の必要がある。式(43)を時間微分すると次式を得る。

$$
\ddot{x}
=
\dot{A}(t)\omega_0\cos\psi(t)
–
A(t)\omega_0^2\sin\psi(t)
–
A(t)\omega_0\dot{\phi}(t)\sin\psi(t) \tag{46}
$$

これを式(34)に代入すると、

$$
\dot{A}(t)\omega_0\cos\psi(t)
–
A(t)\omega_0\dot{\phi}(t)\sin\psi(t)
=
-\frac{\epsilon}{m}f(x,\dot{x}) \tag{47}
$$

式(45)と式(47)を組み合わせると、

$$
\dot{A}(t)
=
-\frac{\epsilon}{m\omega_0}
f(x,\dot{x})\cos\psi \tag{48}
$$

$$
\dot{\phi}(t)
=
\frac{\epsilon}{m\omega_0A(t)}
f(x,\dot{x})\sin\psi \tag{49}
$$

ここで \(f(x,\dot{x})\) の \(x\)、\(\dot{x}\) には、\(A(t)\sin\{\omega_0t+\phi(t)\}\)、\(A(t)\omega_0\cos\{\omega_0t+\phi(t)\}\) を用いる。ここまでは近似は行われておらず、厳密な関係式である。この式から厳密解を得るのは難しく、近似的な評価をする。

先に議論したように、\(A(t)\)、\(\phi(t)\) は「ゆっくり」変化する場合を考えている。振幅と位相の時間変化を周期ごとの時間平均により近似する。

$$
A(t_0+T_0)
=
A(t_0)\{1+O(\epsilon)\} \tag{50}
$$

$$
\phi(t_0+T_0)
=
\phi(t_0)\{1+O(\epsilon)\} \tag{51}
$$

とすると、平均を求めるときの積分内の \(A\) を \(A(t_0)\) と近似して悪くないはずである。無視した分の誤差は \(O(\epsilon^2)\) となる。^※10

$$
\begin{aligned}
\dot{\bar{A}}(t_0)
&=
\frac{1}{T_0}
\int_{t_0}^{t_0+T_0}
\dot{A}(t)\,dt\\
&=
\frac{1}{2\pi}
\int_0^{2\pi}
\dot{A}(t)\,d\psi
\end{aligned}
\tag{52}
$$

$$
\dot{\bar{A}}(t_0)
\simeq
-\frac{\epsilon}{2\pi m\omega_0}
\int_0^{2\pi}
f\left(
A(t_0)\sin\psi,
A(t_0)\omega_0\cos\psi
\right)
\cos\psi\,d\psi \tag{53}
$$

$$
\begin{aligned}
\dot{\bar{\phi}}(t_0)
&=
\frac{1}{T_0}
\int_{t_0}^{t_0+T_0}
\dot{\phi}(t)\,dt\\
&=
\frac{1}{2\pi}
\int_0^{2\pi}
\dot{\phi}(t)\,d\psi
\end{aligned}
\tag{54}
$$

$$
\dot{\bar{\phi}}(t_0)
\simeq
\frac{\epsilon}{2\pi m\omega_0A(t_0)}
\int_0^{2\pi}
f\left(
A(t_0)\sin\psi,
A(t_0)\omega_0\cos\psi
\right)
\sin\psi\,d\psi \tag{55}
$$

蛇足として、

$$
\dot{\psi}(t)
=
\dot{\phi}(t)+\omega_0 \tag{56}
$$

である。上線は時間平均されたことを示す。^※11平均法の一般論はいったんここまでとし、続きは後の章で議論する。

2.2　平均法の線形減衰振動系への適用と厳密解との比較

ここでは線形減衰振動系に適用する。摂動項は式(43)より、

$$
\frac{\epsilon}{m}f(x,\dot{x})
=
2\epsilon\zeta\omega_0\dot{x}
=
2\epsilon\zeta\omega_0^2A(t)\cos\psi(t) \tag{57}
$$

となる。

$$
\int_0^{2\pi}\cos^2\psi\,d\psi
=
\pi,
\qquad
\int_0^{2\pi}\sin\psi\cos\psi\,d\psi
=
0 \tag{58}
$$

を使って、式(53)と式(55)は、

$$
\frac{\dot{\bar{A}}(t)}{\bar{A}(t)}
=
-\epsilon\zeta\omega_0 \tag{59}
$$

$$
\dot{\bar{\phi}}(t)
=
0 \tag{60}
$$

となる。これを解いて、

$$
\bar{A}(t)
=
\bar{A}(t_0)
\exp(-\epsilon\zeta\omega_0t) \tag{61}
$$

$$
\bar{\phi}(t)
=
\bar{\phi}(0) \tag{62}
$$

これは、式(30)の振幅の時間依存項に \(O(\epsilon)\) まで合致している。また、周波数、位相はこの次数では摂動の有無で変わらない。

3　非線形振動

3.1　非線形保存振動系

3.1.1　非線形保存振動系の数値解

非摂動系は前と同じ。

$$
\ddot{x}
+
\omega_0^2x
+
\frac{\epsilon}{m}f(x,\dot{x})
=
0
\qquad\text{再掲} \tag{63}
$$

数値的に解いた結果が図9、図10である。図9は時間発展の波形で、単振動の場合と \(x^3\) 項を加えた場合を比較した。周波数がずれているのがわかる。

ちなみに振幅が小さい、例えば初期条件 \(x(0)=0.1\) だと、位相のずれはほとんど見られない。後の式(72)でわかるように、振幅の2乗で効いてくるからである。

図10では位相平面を描いた。\(x\) のみの系と比較して、速度方向が拡大されている。線形振動子、すなわち調和振動子系と比較して、\(x^3\) 項により加速されているからである。

3.1.2　非線形保存振動系の近似解

ここでは、線形非保存系に引き続き、平均法により近似解を求める。非摂動系は \(\epsilon=0\) のときで、その解と時間微分は、

$$
x
=
A_0\sin(\omega_0t+\phi_0)
\qquad\text{再掲} \tag{64}
$$

$$
\dot{x}
=
A_0\omega_0\cos(\omega_0t+\phi_0) \tag{65}
$$

非線形保存振動系の例として、ダッフィング（Duffing）系といわれる系を考える。\(f\) は、

$$
f(x)=kx^3 \tag{66}
$$

である。式(53)、式(55)に \(f\) を代入すると、

$$
\dot{\bar{A}}(t_0)
=
-\frac{\epsilon}{2\pi m\omega_0}
\int_0^{2\pi}
f\left(
A(t_0)\sin\psi
\right)
\cos\psi\,d\psi \tag{67}
$$

$$
\dot{\bar{A}}(t_0)
=
-\frac{\epsilon kA(t_0)^3}{2\pi m\omega_0}
\int_0^{2\pi}
\sin^3\psi\cos\psi\,d\psi
=
0 \tag{68}
$$

$$
\dot{\bar{\phi}}(t_0)
=
\frac{\epsilon}{2\pi m\omega_0A(t_0)}
\int_0^{2\pi}
f\left(
A(t_0)\sin\psi
\right)
\sin\psi\,d\psi \tag{69}
$$

$$
\dot{\bar{\phi}}(t_0)
=
\frac{\epsilon kA(t_0)^2}{2\pi m\omega_0}
\int_0^{2\pi}
\sin^4\psi\,d\psi \tag{70}
$$

$$
\dot{\bar{\phi}}(t_0)
=
\frac{3}{8}
\epsilon A(t_0)^2\omega_0 \tag{71}
$$

よって一般解は、\(A\)、\(\phi\) を任意の定数として、

$$
x(t)
=
A\sin
\left\{
\left(
1+\frac{3\epsilon}{8}A^2
\right)
\omega_0t+\phi
\right\} \tag{72}
$$

初期条件を \(x(0)=a_0\)、\(\dot{x}(0)=0\) とすると、

$$
x(t)
=
a_0\sin
\left\{
\left(
1+\frac{3\epsilon}{8}a_0^2
\right)
\omega_0t+\frac{\pi}{2}
\right\} \tag{73}
$$

となる。結果における注意点として、振動数が \(\omega_0\) から変化しており、振幅に依存している。前者は線形非保存振動系でも見られたが、後者の振幅に依存するという性質は非線形な系での特徴といわれる。

補足：

$$
\int_0^{2\pi}\sin^4\psi\,d\psi
=
\frac{1}{4}
\int_0^{2\pi}
\left(
1-2\cos2\psi+\cos^22\psi
\right)
d\psi \tag{74}
$$

$$
\begin{aligned}
\int_0^{2\pi}\sin^4\psi\,d\psi
&=
\frac{1}{4}
\int_0^{2\pi}
\left[
1-2\cos2\psi
+
\frac{1}{2}(1+\cos4\psi)
\right]
d\psi\\
&=
\frac{3}{8}\,2\pi
\end{aligned}
\tag{75}
$$

3.2　非線形散逸振動系

3.2.1　非線形散逸振動系の数値解

非線形散逸振動系の例として、動摩擦を受ける系を考える。\(f\) は、

$$
f(x,\dot{x})
=
\sum_{j=1}^{n}
f^{(j)}\dot{x}^{\,j} \tag{76}
$$

である。ここでは、\(f^{(j)}\neq0\) for \(j=1,3\)、その他の \(j\) では0のときを考える。まず、数値的に解いてみる。

図11では、初期条件が小さいところから振幅が増幅し、一定の振幅に収束しているように見える。

図12は同じ波形を位相平面上に描いたものである。時間が経過すると同じ軌跡上に収束していく様子がわかる。

図13では、初期条件が大きいところから振幅が小さくなり、一定の振幅に収束しているように見える。

図14では同じ波形を位相平面上に描き、初期条件が小さい場合と比較した。時間が経過すると同じ軌跡上に収束していく様子がわかる。この軌跡はリミットサイクルと呼ばれる。

3.2.2　非線形散逸振動系の近似解

ここで、近似解ではあるが解析的な解を求めるために、平均法を適用する。式(53)、式(55)に \(f\) を代入すると、

$$
\dot{\bar{A}}(t_0)
=
-\frac{\epsilon}{2\pi m\omega_0}
\int_0^{2\pi}
f\left(
A(t_0)\sin\psi,
A(t_0)\omega_0\cos\psi
\right)
\cos\psi\,d\psi \tag{77}
$$

$$
\dot{\bar{A}}^{(j=1)}(t_0) \tag{78}
$$

$$
\dot{\bar{A}}^{(j=1)}(t_0)
\equiv
-\frac{\epsilon}{2\pi m\omega_0}
\int_0^{2\pi}
f^{(1)}\dot{x}\cos\psi\,d\psi \tag{79}
$$

$$
\dot{\bar{A}}^{(j=1)}(t_0)
=
-\frac{\epsilon f^{(1)}}{2m}A_0
=
-\omega_0\nu A_0,
\qquad
\nu
=
\frac{\epsilon f^{(1)}}{2\sqrt{mk}} \tag{80}
$$

$$
\dot{\bar{A}}^{(j=3)}(t_0) \tag{81}
$$

$$
\dot{\bar{A}}^{(j=3)}(t_0)
\equiv
-\frac{\epsilon}{2\pi m\omega_0}
\int_0^{2\pi}
f^{(3)}\dot{x}^{\,3}\cos\psi\,d\psi \tag{82}
$$

$$
\dot{\bar{A}}^{(j=3)}(t_0)
=
-\frac{\epsilon f^{(3)}A_0^3\omega_0^2}{2\pi m}
\int_0^{2\pi}
\cos^4\psi\,d\psi \tag{83}
$$

$$
\dot{\bar{A}}^{(j=3)}(t_0)
=
-\frac{3}{8m}
\epsilon f^{(3)}A_0^3\omega_0^2
=
-\frac{3}{4}
\omega_0^3\nu
\frac{f^{(3)}}{f^{(1)}}A_0^3 \tag{84}
$$

$$
\dot{\bar{\phi}}(t_0)
=
\frac{\epsilon}{2\pi m\omega_0A(t_0)}
\int_0^{2\pi}
f\left(
A(t_0)\sin\psi,
A(t_0)\omega_0\cos\psi
\right)
\sin\psi\,d\psi \tag{85}
$$

$$
\dot{\bar{\phi}}^{(j=1)}(t_0)
\equiv
\frac{\epsilon}{2\pi m\omega_0A(t_0)}
\int_0^{2\pi}
f^{(1)}\dot{x}\sin\psi\,d\psi
=
0 \tag{86}
$$

$$
\dot{\bar{\phi}}^{(j=3)}(t_0)
\equiv
\frac{\epsilon}{2\pi m\omega_0A(t_0)}
\int_0^{2\pi}
f^{(3)}\dot{x}^{\,3}\sin\psi\,d\psi
=
0 \tag{87}
$$

まとめると、

$$
\begin{aligned}
\dot{\bar{A}}(t)
&=
-\frac{\epsilon f^{(1)}A_0}{2m}
\left(
1+
\frac{3}{4}
\frac{f^{(3)}}{f^{(1)}}
\omega_0^2A_0^2
\right)\\
&=
-\omega_0\nu A_0
\left(
1+
\frac{3\epsilon f^{(3)}\omega_0A_0^2}{8m\nu}
\right)
\end{aligned}
\tag{88}
$$

$$
\dot{\bar{\phi}}(t)
=
0 \tag{89}
$$

ダッフィング系の \(x^3\) では周波数がずれるのに、\(\dot{x}^3\) ではこの次数の計算では周波数は \(\omega_0\) のままである。これは、平均法の計算で \(O(\epsilon)\) までしか計算していないからと思われる。^※12

式(88)の解を求める。

$$
\frac{\tilde{\epsilon}}{2}
\equiv
\frac{\epsilon f^{(1)}}{2m} \tag{90}
$$

$$
\tilde{f}
\equiv
\frac{3f^{(3)}}{4f^{(1)}}\omega_0^2 \tag{91}
$$

$$
\tilde{f}A^2
\equiv
-\frac{1}{4}\tilde{A}^{\,2},
\qquad
A
=
\frac{\tilde{A}}{\sqrt{-4\tilde{f}}} \tag{92}
$$

と置くと、式(88)は、

$$
\dot{\tilde{A}}
=
-\frac{\tilde{\epsilon}}{2}
\tilde{A}
\left(
1-\frac{1}{4}\tilde{A}^{\,2}
\right) \tag{93}
$$

$$
\tilde{A}(t)
=
\frac{
\tilde{A}(0)\exp(-\tilde{\epsilon}t/2)
}{
\sqrt{
1+
\frac{1}{4}\tilde{A}(0)^2
\left[
\exp(-\tilde{\epsilon}t)-1
\right]
}
} \tag{94}
$$

となる。文献[2]の168ページ、式(13.2)に対応する。ただし、\(\tilde{\epsilon}\) の前の符号に注意する。\(f^{(1)}\) を正で定義したため、文献[2]と逆符号になっている。

よって、

$$
A(t)
=
\frac{
A(0)
\exp\left(
-\dfrac{\epsilon f^{(1)}}{2m}t
\right)
}{
\sqrt{
1
–
\dfrac{3f^{(3)}\omega_0^2}{4f^{(1)}}
A(0)^2
\left[
\exp\left(
-\dfrac{\epsilon f^{(1)}}{m}t
\right)
-1
\right]
}
} \tag{95}
$$

注意：分母の指数関数の肩は \(1/m\) であり、\(1/(2m)\) ではない。図15に、微分方程式を数値的に求めた解と、この近似解を比較した。ほぼ合致している。\(f^{(1)}\) をもっと大きくすれば、ずれが大きくなることを数値的に確認している。

4　いろんな自励振動

自励振動は、ここで挙げた負の勾配を持つ摩擦係数によるもの以外にも、いろいろなパターンがある。自励振動の例として、文献[7]の158、159ページを参考文献として挙げておく。^※13

参考文献

1. 柴田正和『漸近級数と特異摂動法』森北出版。
2. 安田仁彦『振動工学 応用編』コロナ社。
3. 岩田佳雄・佐伯暢人・小松崎俊彦『機械振動学』数理工学社、新・数理／工学ライブラリ 機械工学5。
4. 井上順吉・末岡淳男『機械力学II―非線形振動論―』理工学社、機械工学基礎講座。
5. スティーヴン・ストロガッツ著、田中久陽・中尾裕也・千葉逸人訳『非線形ダイナミクスとカオス』丸善出版。
6. ホームページ「楽しい物理ノート」力学・振動／対話・非線形振動。
7. 中川憲治・室津義定・岩壺卓三『工業振動学』第2版、森北出版。

脚注

1. 参考文献は私が知る文献の中で、その項目の説明が勧められるものである。文献が多くあるが、特筆すべきものを知らない場合は、特に文献を挙げない。 ↩
2. ちょこちょこ自信のないところがあります。その場合は書いてあります。 ↩
3. \(\omega\) はオメガ、\(\phi\) はファイと読む。 ↩
4. 式(16)で最初から \(x=X-X_0-L_0+f^{(0)}/k\) と置けば、この形になる。 ↩
5. \(\zeta\) は「ツェータ」と読む。 ↩
6. \(a\ll1\) のとき、\(a^2\ll a\ll1\)、\(1/(1-a)=1+a+a^2+a^3+\cdots\) が成り立つ。 ↩
7. 例えば「\(\epsilon\zeta\ll1\) のとき、対数減衰率を \(\epsilon\mu/\sqrt{mk}\) として、バネ定数が大きいと減衰が小さい」とする文献が見られるが、誤っている。 ↩
8. \(\epsilon\) を付けて計算した方が見通しがよい。 ↩
9. 常微分方程式の解法のひとつである定数変化法を利用している。\(A_0\)、\(\phi_0\) がゆっくり変化するという特性はここでは使わなくてよく、この段階では近似はしていない。詳しい説明は文献[6]がよいと思う。 ↩
10. ちょっと怪しい。確認のこと。 ↩
11. 積分内の \(\psi\) は平均化された記号に置き換えなくてよい。位相を置き換えると、例えば式(68)が0にならない可能性がある。 ↩
12. どこかの本に書いてあったけど、確認したい。2次の計算はしていない。 ↩
13. 理解はできていない。2体間の位相差によるエネルギーのやり取りがみそらしい。 ↩